I want to implement a minimisation on a 2D spin Ising model with 30x30 grid. The spin variables is 0,1 and the objective is to minimize the sum of products of spins. For simplicity, I only include NN pair interaction and a "triangular" interaction term. The formulation I give is as followed. And I wish to ask for suggestions whether there would be some stronger formulation that could make computation faster. Thank you :D


the original objective is:


$$\sum\limits\_{i = 1}^{29} {\sum\limits\_{j = 1}^{30} {{J\_1}{s\_{i,j}}{s\_{i + 1,j}}} } + \sum\limits\_{i = 1}^{30} {\sum\limits\_{j = 1}^{29} {{J\_1}{s\_{i,j}}{s\_{i,j + 1}}} } + \sum\limits\_{i = 1}^{29} {\sum\limits\_{j = 1}^{29} {{J\_3}{s\_{i,j}}{s\_{i + 1,j}}{s\_{i,j + 1}}} } $$


$${s\_{i,j}} \in \{ 0,1\} $$


we linearize it so that:


[\begin{align}
 & \min : \\ 
 & \sum\limits\_{i=1}^{29}{\sum\limits\_{j=1}^{30}{{{J}\_{1}}{{s}\_{i,j,1}}}}+\sum\limits\_{i=1}^{30}{\sum\limits\_{j=1}^{29}{{{J}\_{1}}{{s}\_{i,j,2}}}}+\sum\limits\_{i=1}^{29}{\sum\limits\_{j=1}^{29}{{{J}\_{3}}{{s}\_{i,j,3}}}} \\ 
 & \text{subject to:} \\ 
 & \forall i=1...29\,j=1...30\,\ k=1...4 \\ 
 & \left[ \begin{matrix}
 {{s}\_{i,j}} \\
 {{s}\_{i+1,j}} \\
 {{s}\_{i,j,1}} \\
\end{matrix} \right]={{y}\_{i,j,1,1}}\left[ \begin{matrix}
 0 \\
 0 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,1,2}}\left[ \begin{matrix}
 1 \\
 0 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,1,3}}\left[ \begin{matrix}
 0 \\
 1 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,1,4}}\left[ \begin{matrix}
 1 \\
 1 \\
 1 \\
\end{matrix} \right] \\ 
 & \sum\limits\_{k=1}^{4}{{{y}\_{i,j,1,k}}=1} \\ 
 & 0\le y\le 1 \\ 
 & {{y}\_{i,j,1,k}}\text{ is integral} \\ 
 & \\ 
 & \forall i=1...30\,j=1...29\,\,k=1...4 \\ 
 & \left[ \begin{matrix}
 {{s}\_{i,j}} \\
 {{s}\_{i,j+1}} \\
 {{s}\_{i,j,1}} \\
\end{matrix} \right]={{y}\_{i,j,2,1}}\left[ \begin{matrix}
 0 \\
 0 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,2,2}}\left[ \begin{matrix}
 1 \\
 0 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,2,3}}\left[ \begin{matrix}
 0 \\
 1 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,2,4}}\left[ \begin{matrix}
 1 \\
 1 \\
 1 \\
\end{matrix} \right] \\ 
 & \sum\limits\_{k=1}^{4}{{{y}\_{i,j,2,k}}=1} \\ 
 & 0\le {{y}\_{i,j,2,k}}\le 1 \\ 
 & {{y}\_{i,j,2,k}}\text{ is integral} \\ 
 & \\ 
 & \forall i=1...29\,j=1...29\,\,k=1...8 \\ 
 & \left[ \begin{matrix}
 {{s}\_{i,j}} \\
 {{s}\_{i+1,j}} \\
 {{s}\_{i,j+1}} \\
 {{s}\_{i,j,3}} \\
\end{matrix} \right]={{y}\_{i,j,3,1}}\left[ \begin{matrix}
 0 \\
 0 \\
 0 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,3,2}}\left[ \begin{matrix}
 1 \\
 0 \\
 0 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,3,3}}\left[ \begin{matrix}
 0 \\
 1 \\
 0 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,3,4}}\left[ \begin{matrix}
 1 \\
 1 \\
 0 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,3,5}}\left[ \begin{matrix}
 0 \\
 0 \\
 1 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,3,6}}\left[ \begin{matrix}
 1 \\
 0 \\
 1 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,3,7}}\left[ \begin{matrix}
 0 \\
 1 \\
 1 \\
 0 \\
\end{matrix} \right]+{{y}\_{i,j,3,8}}\left[ \begin{matrix}
 1 \\
 1 \\
 1 \\
 1 \\
\end{matrix} \right] \\ 
 & \sum\limits\_{k=1}^{8}{{{y}\_{i,j,3,k}}=1} \\ 
 & 0\le {{y}\_{i,j,3,k}}\le 1 \\ 
 & {{y}\_{i,j,3,k}}\text{ is integral} \\ 
 & \\ 
 & \forall i=1...30\,j=1...30 \\ 
 & 0\le {{s}\_{i,j}}\le 1 \\ 
 & {{s}\_{i,j}}\quad \text{integral} \\ 
\end{align}]



Mixed integer programming formulation for Ising model
This problem is studied in Condensed Matter. You may wish to look at Projected Entangled Pair States - the 2D generalization of famous Density Matrix Renormalization Group / Matrix Product States ansatz. The recent development is reported in <http://arxiv.org/abs/1405.3259> but do consider references therein as well.


