*Note: This post was updated to work with Python 3 and PySCIPOpt. The original version used Python 2 and python-zibopt. It has also been edited for clarity.*

Back in October of 2011, I started toying with a model for finding magic squares using SCIP. This is a fun modeling exercise and a challenging problem. First one constructs a square matrix of integer-valued variables.

```
from pyscipopt import Model
# [...snip...]
m = Model()
matrix = []
for i in range(size):
row = [m.addVar(vtype="I", lb=1) for _ in range(size)]
for x in row:
m.addCons(x <= M)
matrix.append(row)
```

Then one adds the following constraints:

- All variables β₯ 1.
- All rows, columns, and the diagonal sum to the same value.
- All variables take different values.

The first two constraints are trivial to implement, and relatively easy for the solver. What I do is add a single extra variable then set it equal to the sums of each row, column, and the diagonal.

```
sum_val = m.addVar(vtype="M")
for i in range(size):
m.addCons(sum(matrix[i]) == sum_val)
m.addCons(sum(matrix[j][i] for j in range(size)) == sum_val)
m.addCons(sum(matrix[i][i] for i in range(size)) == sum_val)
```

It’s the third that messes things up. You can think of this as saying, for every possible pair of integer-valued variables $x$ and $y$:

$$ x \ge y + 1 \quad \text{or} \quad x \le y - 1 $$

Why is this hard? Because we can’t add both constraints to the model. That would make it infeasible. Instead, we add write them in such a way that exactly one will be active for any any given solution. This requires, for each pair of variables, an additional binary variable $z$ and a (possibly big) constant $M$. Thus we reformulate the above as:

$$ x \ge (y + 1) - M z \ x \le (y - 1) + M (1-z) \ z \in {0,1} $$

In code this looks like:

```
from itertools import chain
all_vars = list(chain(*matrix))
for i, x in enumerate(all_vars):
for y in all_vars[i+1:]:
z = m.addVar(vtype="B")
m.addCons(x >= y + 1 - M*z)
m.addCons(x <= y - 1 + M*(1-z))
```

However, here be dragons. We may not know how big (or small) to make $M$. Generally we want it as small as possible to make the LP relaxation of our integer programming model tighter. Different values of $M$ have unpredictable effects on solution time.

Which brings us to an interesting idea:

SCIP now supports bilinear constraints. This means that I can make $M$ a variable in the above model.

```
import sys
try:
M = int(sys.argv[2])
except IndexError:
M = m.addVar(vtype="M", lb=size * size)
else:
assert M >= size * size
```

The magic square model linked to in this post provides both options. The first command line argument it requires is the matrix size. The second one, $M$, is optional. If not given, it leaves $M$ up to the solver.

An interesting exercise for the reader: Change the code to search for a *minimal* magic square, which minimizes either the value of $M$ or the sums of the columns, rows, and diagonal.