Roughly speaking, the system has stable configurations where the free magnetic energy $W_{\rm mag}$ assumes a strict local minimum.
That local minimum of $W_{\rm mag}$ is assumed where the combined permeance of the magnet and the iron plate from pole to pole has a strict local maximum.
Symmetry, is not always a valid argument in this regard. E.g., assume a circular iron ring with a hole that is larger than the magnet. Then the central position of the magnet (i.e., magnet in hole) will not be a local maximum of combined permeance. It will be a local minimum. That means the central configuration will be an equilibrium, but an unstable one.
Pityingly, this description is rather rough since it works with concentrated parameters like permeance and the problem actually has distributed parameters.
$\def\rmHS{{\rm HS}}$ $\def\rmFe{{\rm Fe}}$ $\def\rmm{{\rm m}}$ $\def\vr{\vec{r}}$ $\def\vH{\vec{H}}$ $\def\vB{\vec{B}}$ $\def\vA{\vec{A}}$ $\def\nR{{\mathbb{R}}}$ $\def\ph{{\varphi}}$ $\def\grad{\mathrm{grad}}$ $\def\div{\mathrm{div}}$ A simple example demonstrates how it works with concentrated parameters. We discuss a horseshoe magnet standing on an iron bar (see the following Figure).
We assume that the field within the horseshoe magnet is almost homogeneous in longitudinal direction and denote the longitudinal field strength with $H_\rmHS$ and the flux density with $B_\rmHS$ (that works quite good if the thickness of the arc segment is small enough).
The field in the magnet can be described by $$ B_\rmHS = \mu_0 \left(H_\rmHS + M_\rmHS\right) $$ with the magnetization $M_\rmHS$ of the horseshoe magnet.
Thereby, we assume that the changes in magnetization displacement on the iron bar are so small that $M_\rmHS$ can be approximated as a linear-affine function of the field strength: $$ M_\rmHS = M_{\rmHS0} + \chi H_\rmHS $$ We define the flux $\Phi:=A_\rmHS B_\rmHS$ leaving the horseshoe at $x_2$ and the magnetic voltage drop as the integral of $H_\rmHS$ from $x_1$ to $x_2$ over the mean path inside the horseshoe $V_\rmm:= l_\rmHS H_\rmHS$. Substituting this into the last two formulae gives $$ \Phi = A_\rmHS B_\rmHS = A_\rmHS\mu_0\left((1+\chi)\frac{V_\rmm}{l_\rmHS} + M_{\rmHS0}\right) = G_\rmHS V_\rmm + \Phi_{\rmHS0} $$ The magnet acts at its poles like a parallel connection of a permeance $G_\rmHS=\frac{(1+\chi)\mu_0A}{l_\rmHS}$ with a flux source $\Phi_0=\mu_0AM_{\rmHS0}$.
On the other side the flux goes into the iron bar (we neglect the stray flux here). The permeance $G_\rmFe(x)$ of the iron bar depends on the position of the magnet on the bar.
With the magnetic voltage drop $-V_\rmm$ from $x_2$ to $x_1$ on the iron permeance we get the equation $$ G_\rmHS V_\rmm + \Phi_0 = \Phi = G_\rmFe(x)\cdot(-V_\rmm) $$ for the flux balance or equivalently $$ V_\rmm = \frac{\Phi_0}{G_\rmHS+G_\rmFe(x)}. $$ The stored free magnetic energy is $$ W_\rmm = \frac12\int_{\nR^n} \vB\cdot\vH d V. $$ We neglect the stray flux and its energy. Therefore, the magnetic energy can be decomposed into the part $W_\rmHS$ stored in the horseshoe magnet and the part $W_\rmFe$ stored in the iron.
In the magnet we get $$ W_\rmHS = \frac12 A_\rmHS l_\rmHS B_\rmHS H_\rmHS = \frac12 V_\rmm \Phi $$ For the iron material we get similarly $$ W_\rmFe = \frac12 V_\rmm\Phi $$ and the overall energy can be calculated in the following way:
$W_\rmm = \frac12 V_\rmm\Phi + \frac12 V_\rmm\Phi$
$\phantom{W_\rmm}= \frac12 V_\rmm(G_\rmHS V_\rmm + \Phi_0) + \frac12 V_\rmm G_\rmFe(x) V_\rmm$
$\phantom{W_\rmm}= \frac12\left(G_\rmHS+G_\rmFe(x)\right)V_\rmm^2 + \frac12 V_\rmm\Phi_0$
and with $V_\rmm = \frac{\Phi_0}{G_\rmHS+G_\rmFe(x)}$: $$ W_\rmm = \frac{\Phi_0^2}{G_\rmHS+G_\rmFe(x)}. $$ And with this formula the magnetic energy has a strong local minimum at a strong local maximum of $G_\rmFe(x)$.
There remains to discuss the magnetic energy within the iron block.
The field within the iron block can be calculated with the help of the magneto-static potential $\ph_\rmm$ with $H_\rmFe = -\grad\ph_\rmm$ defined through the equation $$ 0=\div(\vB)=-\div(\mu\grad\ph_\rmm) $$ and the boundary conditions: $\ph_\rmm(\vr_1)=0$ for points $\vr_1$ on the pole 1, $\ph_\rmm(\vr_2)=V_\rmm$ for points $\vr_2$ on pole 2 and for points $\vr_2$ on pole 2 and $\grad\ph(\vr) \cdot d\vA=0$ at all the other boundary points $\vr$ whereby $d\vA$ is the outwards directed area element. The stored energy is $$ W_\rmFe = \int_{V_\rmFe} \vH\cdot\vB dV = \int_{V_\rmFe} (\grad\ph_\rmm)\cdot(\mu\grad\ph_\rmm) dV $$ With the differentiation rule $\div(\ph_\rmm\mu\grad\ph_\rmm) = (\grad\ph_\rmm)\cdot(\mu\grad\ph_\rmm) + \ph_\rmm \div(\mu\grad\ph_\rmm)$ and Gauß's theorem for the divergence operator one obtains $$ W_\rmFe = \oint_{\partial V_\rmFe} \ph_\rmm \mu (\grad\ph_\rmm) \cdot d\vA - \int_{V_\rmFe} \ph_\rmm \underbrace{\div(\mu\grad\ph_\rmm)}_{=0}d V $$ Thereby, the volume integral is zero because of $\div(\mu\grad\ph_\rmm)=0$.
The surface integral has only a nonzero component at pole 2 since at pole 1 $\ph_\rmm=0$ and on the complement surface $\grad\ph_\rmm \cdot d \vA=0$. Furthermore, on pole 2 the magnetic potential $\ph_\rmm=V_\rmm$ is assumed to be constant. Thus, we have $$ W_\rmFe = (-V_\rmm) \int_{A_2} \mu (-\grad\ph_\rmm)\cdot d\vA = (-V_\rmm) \int_{A_2} \vB \cdot d\vA = - V_\rmm \Phi. $$