In this work, we address the simulation and testing of MEMS micromirrors with hardening and softening behaviour excited with patches of piezoelectric materials. The forces exerted by the piezoelectric patches are modelled by means of the theory of ferroelectrics developed by Landau–Devonshire and are based on the experimentally measured polarisation hysteresis loops. The large rotations experienced by the mirrors also induce geometrical nonlinearities in the formulation up to cubic order. The solution of the proposed model is performed by discretising the device geometry using the Finite Element Method, and the resulting large system of coupled differential equations is solved by means of the Harmonic Balance Method. Numerical results were validated with experimental data collected on the devices.