For the description of the hysteresis behavior of ferroelectrics/ferroelastics under multiaxial combined electrical and/or mechanical loading, a thermodynamically consistent microstructural model of a ferroelectroelastic material is proposed taking into account the presence and evolution of polar point defects. The model also takes into account multiphase composition, anisotropy of properties, domain structure, and dissipative motion of domain walls. The linear theory of the charged point defects evolution is proposed based on the free energy of defects in the quadratic form of the polarization vector and strain tensor of defects. The dependence of the hysteresis loop shift (due internal field bias) on parameters of the free energy of defects is shown. Comparison of computation results with experimental curves of dielectric, mechanical, and electromechanical hysteresis for polycrystalline piezoelectric PZT PIC-151, BaTiO3, single-crystal PMN-PZT and KTS doped with acceptor additives, showed a good agreement.