This work covers the variational-based modeling of magneto-mechanical hysteresis effects in hard magnetic magneto-active polymers (MAPs). We discuss basic ingredients of the constitutive theory within the concept of generalized standard materials that necessitates suitable definitions of (i) the total energy density function and (ii) the dissipation potential. A key feature of the developed energy-based constitutive model is its modular structure that allows to associate individual components of the total energy density with the magnetic field strength in vacuum, the magnetization of the material, and specific components of the total stress tensor. The utilization of spectral invariants allows the formulation of a compact constitutive model that is able to accurately capture the highly nonlinear material behavior of hard magnetic MAPs with stochastic microstructures. The developed constitutive functions are subsequently embedded in incremental variational principles. This is supplemented by the discussion of corresponding conforming finite element methods. The performance of the developed variational-based models is demonstrated by solving some application-oriented initial boundary value problems.