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.