Yttria stabilized zirconia (YSZ) ceramics have been used for various engineering applications including structural ceramics, biomedical materials, and thermal barrier coatings. The versatile and excellent properties of YSZ stem from its unique microstructure consisting of monoclinic, tetragonal, and cubic phases, whose stability depends on yttria concentration and temperature. However, there are no empirical interatomic potentials (EIPs) that can reproduce the structures and energies of ZrO₂ and YSZ polymorphs, limiting the atomic-scale investigation of lattice defect structures and their interactions that affect the YSZ microstructure and properties. Here, using a genetic algorithm and ab initio training datasets, we have optimized EIPs to sufficiently reproduce the structures and stability of ZrO₂ and YSZ polymorphs, as well as the properties of the tetragonal and cubic phases at finite temperature. The potentials have also been applied to the search for a tetragonal grain boundary structure, showing that the obtained grain boundary structure is consistent with that obtained by ab initio calculations. The developed EIPs will aid in revealing the microstructure-property relationships in YSZ by performing large-scale and systematic calculations, which are practically difficult to perform with ab initio and machine-learning-potential calculations.