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.