Gekko求解整数规划幻方问题：添加第7个约束后无解
2026-7-3
问题描述
尝试用Gekko的整数规划功能求解3×3幻方：
- 移除元素唯一性约束时，算法能正常找到解；
- 添加最多6个元素两两不等的约束时，程序可正常运行；
- 但添加第7个不等约束后，求解器无法找到可行解。
已知3×3幻方存在有效解，已尝试将约束改写为abs(M[0,0]-M[0,1])>0或平方形式，怀疑是约束数量或形式导致的问题，询问是否可通过调整求解器设置或重构约束解决。
原测试代码
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)
解决方案
1. 重构唯一性约束（推荐）
原代码用两两绝对值约束实现元素唯一性，会引入大量非线性约束，导致求解器分支定界过程中难以搜索到可行解。可以用以下两种更高效的方式替代：
方式1：利用排列的和与平方和约束
3×3幻方的元素是1-9的排列，因此满足：
- 所有元素的和为
1+2+...+9=45 - 所有元素的平方和为
1²+2²+...+9²=285
仅需2个约束即可替代原有的36个两两不等约束，且均为多项式约束，求解器更容易处理。
修改后的约束代码：
# 唯一性约束：元素是1-9的排列 m.Equation(m.sum([M[i][j] for i in range(n) for j in range(n)]) == 45) m.Equation(m.sum([M[i][j]**2 for i in range(n) for j in range(n)]) == 285)
方式2：线性化两两不等约束
对于整数变量a和b，a≠b可通过引入二进制变量线性化，避免非线性的绝对值约束。示例代码如下：
# 遍历所有两两不同的元素对 for i in range(n): for j in range(n): for p in range(i, n): for q in range(j+1, n) if i == p else range(n): if (i,j) != (p,q): # 引入二进制变量，控制两个方向的不等约束 b = m.Var(lb=0, ub=1, integer=True) # M[i][j] >= M[p][q] + 1 或者 M[p][q] >= M[i][j] + 1 m.Equation(M[i][j] - M[p][q] >= 1 - 10*(1 - b)) m.Equation(M[p][q] - M[i][j] >= 1 - 10*b)
这里10是足够大的常数（大于M中元素的最大差值8），确保约束生效。
2. 调整求解器参数
如果坚持使用原约束形式，可以调整APOPT求解器的参数，提升搜索能力：
# 设置求解器为APOPT（整数规划专用） m.options.SOLVER = 1 # 增加最大迭代次数 m.options.MAX_ITER = 10000 # 调整分支定界的节点数上限 m.options.NODES = 1000 # 放宽收敛容差 m.options.TOL = 1e-4 # 开启详细输出，便于调试 m.options.DIAGLEVEL = 2
3. 更换求解器
尝试使用Gekko的BPOPT求解器（适用于混合整数非线性规划），修改求解器设置：
m.options.SOLVER = 2
内容的提问来源于stack exchange，提问作者Bob Hesse
相关产品推荐
相关产品推荐
