I'm trying to use the integer programming option to find a magic square. The algorithm finds a solution if I remove the requirement the entries are unique. I can add up to six requirements for entries not being equal but fails to find a solution when I add the seventh. As I know there is a solution to this 3x3 case, is there a way to change the solver settings, or reframe the constraints to get a solution?
from gekko import GEKKO
m=GEKKO(remote=False)
m.options.SOLVER=1
n=3 #size of magic square
S=15 #sum of each row, col, diag
M = m.Array(m.Var,(n,n),lb=1,ub=9,integer=True)
#constraint that no two entries are the same: (hard-coded to find where it fails)
m.Equation(abs(M[0,0]-M[0,1])>=1)
m.Equation(abs(M[0,0]-M[0,2])>=1)
m.Equation(abs(M[0,0]-M[1,0])>=1)
m.Equation(abs(M[0,0]-M[1,1])>=1)
m.Equation(abs(M[0,0]-M[1,2])>=1)
m.Equation(abs(M[0,0]-M[2,0])>=1)
m.Equation(abs(M[0,0]-M[2,1])>=1)
#m.Equation(abs(M[0,0]-M[2,2])>=1)
#m.Equation(abs(M[0,1]-M[0,2])>=1)
#m.Equation(abs(M[0,1]-M[1,0])>=1)
#m.Equation(abs(M[0,1]-M[1,1])>=1)
#m.Equation(abs(M[0,1]-M[1,2])>=1)
#m.Equation(abs(M[0,1]-M[2,0])>=1)
#m.Equation(abs(M[0,1]-M[2,1])>=1)
#m.Equation(abs(M[0,1]-M[2,2])>=1)
#m.Equation(abs(M[0,2]-M[1,0])>=1)
#m.Equation(abs(M[0,2]-M[1,1])>=1)
#m.Equation(abs(M[0,2]-M[1,2])>=1)
#m.Equation(abs(M[0,2]-M[2,0])>=1)
#m.Equation(abs(M[0,2]-M[2,1])>=1)
#m.Equation(abs(M[0,2]-M[2,2])>=1)
#m.Equation(abs(M[1,0]-M[1,1])>=1)
#m.Equation(abs(M[1,0]-M[1,2])>=1)
#m.Equation(abs(M[1,0]-M[2,0])>=1)
#m.Equation(abs(M[1,0]-M[2,1])>=1)
#m.Equation(abs(M[1,0]-M[2,2])>=1)
#m.Equation(abs(M[1,1]-M[1,2])>=1)
#m.Equation(abs(M[1,1]-M[2,0])>=1)
#m.Equation(abs(M[1,1]-M[2,1])>=1)
#m.Equation(abs(M[1,1]-M[2,2])>=1)
#m.Equation(abs(M[1,2]-M[2,0])>=1)
#m.Equation(abs(M[1,2]-M[2,1])>=1)
#m.Equation(abs(M[1,2]-M[2,2])>=1)
#m.Equation(abs(M[2,0]-M[2,1])>=1)
#m.Equation(abs(M[2,0]-M[2,2])>=1)
#m.Equation(abs(M[2,1]-M[2,2])>=1)
#Each col sums to S
for i in range(n):
m.Equation(m.sum([M[i,j] for j in range(n)])==S)
#Each row sum to S
for j in range(n):
m.Equation(m.sum([M[i,j] for i in range(n)])==S)
#Diagonals sum to S
m.Equation(m.sum([M[i,i] for i in range(n)])==S)
m.Equation(m.sum([M[i,n-1-i] for i in range(n)])==S)
m.solve()
print(M)
I've tried rewriting the constraints as abs(M[0,0] - M[0,1])>0. I've also changed to squares. It seems the issue is I've hit an upper bound on the number of constraints--it doesn't mention that in the output, but wondering if that is the issue.
Reinderien is correct that it is often easier to translate to a LP problem, even for the MINLP solver (APOPT) in GEKKO. If you do want to use nonlinear models, use m.abs3() or m.abs2() instead of abs() that doesn't work well with gradient-based solvers. Here is an MILP version of the program in Gekko.
from gekko import GEKKO
n = 3 # size of the magic square
S = n * (n*n + 1) // 2
# Create a Gekko model
m = GEKKO(remote=False)
# Create matrix of decision variables
selectors = m.Array(m.Var,(n,n,n*n),lb=0,ub=1,integer=True)
# Sum of each row, column must equal S
for i in range(n):
m.Equation(m.sum([selectors[i][j][k]*(k+1) for j in range(n) for k in range(n*n)]) == S) # Rows
m.Equation(m.sum([selectors[j][i][k]*(k+1) for j in range(n) for k in range(n*n)]) == S) # Columns
m.Equation(m.sum([selectors[i][i][k]*(k+1) for i in range(n) for k in range(n*n)]) == S) # Diagonal
m.Equation(m.sum([selectors[i][n-i-1][k]*(k+1) for i in range(n) for k in range(n*n)]) == S) # Anti-diagonal
# Each number must be used exactly once
for k in range(n*n):
m.Equation(sum(selectors[i][j][k] for i in range(n) for j in range(n)) == 1)
# Each number can only be used once for each entry
# Avoid solutions for 3x3 such as:
# [4, 2, 9]
# [7+3, 5, 0]
# [1, 8, 6]
for i in range(n):
for j in range(n):
m.Equation(sum(selectors[i][j][k] for k in range(n*n)) == 1)
# Solve the model
m.options.SOLVER=1
m.solve(disp=True)
# Extract and print the solution
solution = [[sum((k+1) * selectors[i][j][k].value[0] for k in range(n*n)) for j in range(n)] for i in range(n)]
for row in solution:
print(row)
It produces the following solutions:
3x3 Magic Square
----------------------------------------------
Steady State Optimization with APOPT Solver
----------------------------------------------
Iter: 1 I: 0 Tm: 0.01 NLPi: 2 Dpth: 0 Lvs: 3 Obj: 0.00E+00 Gap: NaN
Iter: 2 I: -1 Tm: 0.00 NLPi: 1 Dpth: 1 Lvs: 2 Obj: 0.00E+00 Gap: NaN
Iter: 3 I: 0 Tm: 0.01 NLPi: 2 Dpth: 1 Lvs: 4 Obj: 0.00E+00 Gap: NaN
Iter: 4 I: -1 Tm: 0.00 NLPi: 1 Dpth: 2 Lvs: 3 Obj: 0.00E+00 Gap: NaN
--Integer Solution: 0.00E+00 Lowest Leaf: 0.00E+00 Gap: 0.00E+00
Iter: 5 I: 0 Tm: 0.01 NLPi: 2 Dpth: 2 Lvs: 3 Obj: 0.00E+00 Gap: 0.00E+00
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 0.05400000000000005 sec
Objective : 0.
Successful solution
---------------------------------------------------
[8.0, 1.0, 6.0]
[3.0, 5.0, 7.0]
[4.0, 9.0, 2.0]
4x4 Magic Square
Iter: 179 I: 0 Tm: 0.01 NLPi: 2 Dpth: 11 Lvs: 8 Obj: 0.00E+00 Gap: 0.00E+00
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 29.3463 sec
Objective : 0.
Successful solution
---------------------------------------------------
[11.0, 13.0, 6.0, 4.0]
[8.0, 2.0, 9.0, 15.0]
[1.0, 7.0, 16.0, 10.0]
[14.0, 12.0, 3.0, 5.0]
The 3x3 solution is very fast, but the 4x4 (and likely larger problems) will take much longer to run with the APOPT solver. A dedicated MILP solver may work better for larger problems with number of binary variables:
For the 4x4 problem, I tried without using m.sum() and used the Python sum() function with faster results (29 sec vs 113 sec).
from gekko import GEKKO
n = 3 # size of the magic square
S = n * (n*n + 1) // 2
# Create a Gekko model
m = GEKKO(remote=False)
# Create matrix of decision variables
selectors = m.Array(m.Var,(n,n,n*n),lb=0,ub=1,integer=True)
# Sum of each row, column must equal S
for i in range(n):
m.Equation(sum(selectors[i][j][k]*(k+1) for j in range(n) for k in range(n*n)) == S) # Rows
m.Equation(sum(selectors[j][i][k]*(k+1) for j in range(n) for k in range(n*n)) == S) # Columns
m.Equation(sum(selectors[i][i][k]*(k+1) for i in range(n) for k in range(n*n)) == S) # Diagonal
m.Equation(sum(selectors[i][n-i-1][k]*(k+1) for i in range(n) for k in range(n*n)) == S) # Anti-diagonal
# Each number must be used exactly once
for k in range(n*n):
m.Equation(sum(selectors[i][j][k] for i in range(n) for j in range(n)) == 1)
for i in range(n):
for j in range(n):
m.Equation(sum(selectors[i][j][k] for k in range(n*n)) == 1)
# Solve the model
m.options.SOLVER=1
m.solve(disp=True)
# Extract and print the solution
solution = [[sum((k+1) * selectors[i][j][k].value[0] for k in range(n*n)) for j in range(n)] for i in range(n)]
for row in solution:
print(row)
The m.sum() is used for larger problems if a model equation exceeds 15,000 characters.