This paper proposes an effective numerical method for shape control of bi-directional functionally graded plates (2D-FGPs) with piezoelectric layers. Isogeometric analysis (IGA) based on non-uniform rational B-splines (NURBS) related to third-order shear deformation theory (TSDT) is employed for the static analysis of the 2D-FGPs with piezoelectric layers. The B-spline basis functions are utilized to represent the distribution of the ceramic volume fractions, where the control points placed along the plane corresponding to the ceramic volume fraction and the applied voltages are taken as the design variables. In addition, an improved moth flame optimization algorithm is utilized to solve the optimization problem of minimizing the static shape error, which effectively balances the exploratory and exploitative capabilities of the algorithm. Various numerical examples of square, skew, and dart-shaped 2D-FGPs are analyzed to validate the proposed method and demonstrated the superior mechanical performance of 2D-FGPs over 1D-FGPs.