🔲 Normal Magic Squares

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.

As a followup to the last post, I created another SCIP example for finding Normal Magic Squares. This is similar to solving a Sudoku problem, except that here the number of binary variables depends on the square size. In the case of Sudoku, each cell has 9 binary variables – one for each potential value it might take. For a normal magic square, there are $n^2$ possible values for each cell, $n^2$ cells, and one variable representing the row, column, and diagonal sums. This makes a total of $n^4$ binary variables and one continuous variables in the model.

However, there are no big-Ms.

I think the neat part of this code is in this section:

# Construct an expression for each cell that is the sum of
# its binary variables with their associated coefficients.
sums = []
for row in matrix:
    sums_row = []
    for cell in row:
        sums_row.append(sum((i + 1) * x for i, x in enumerate(cell)))
    sums.append(sums_row)

It creates sums of the $n^2$ variables for each cell with their appropriate coefficients ($1$ to $n^2$) and stores those expressions to make the subsequent constraint creation simpler.

Another interesting exercise for the reader: Change the code to minimize the sum of each column. How does that impact the solution time?