中等规模MINLP问题下Gekko求解比值目标0-1优化异常求助

问题成因拆解 #

小维度场景返回全0解的核心原因 #

• 缺少分母非零约束：当x为全0向量时，分子sum(Ax)和分母Bx同时为0，构成数值计算上的0/0不定式，求解器在连续松弛模式下很容易锚定到这个无定义但数值可近似的点，输出全0的“可行解”。
• 错误的求解器配置：minlp_as_nlp 1的作用是把混合整数非线性规划完全松弛为连续非线性规划求解，不触发整数分支定界逻辑，相当于把0-1变量放开到[0,1]连续区间，求解器自然可以滑到边界上的全0点。
• 从附的求解日志可以看到，首次迭代的最优目标对应比值2.17，但最终输出全0解、目标值趋近于0，本质就是松弛求解器没有强制分母为正，直接落到了平凡解上。

大维度场景触发Max Equation Length错误的核心原因 #

• 长表达式拼接：df_num.mul(x).sum(1).sum()这类写法会把矩阵逐元素相乘的所有加项拼接成单条APM方程，当列数、行数上升到一定规模后，拼接后的字符串长度会突破APM单条方程15000字符的硬上限。
• 没有利用矩阵的可预计算特性：把所有逐行求和的计算逻辑全部丢给Gekko生成方程，哪怕用了Intermediate分层，也只是把两层求和压成两个长表达式，没有从根源上减少单条方程的长度。
• 从调整后的求解日志看，拆分逐行方程的写法虽然绕开了超长报错，但会把模型规模从2个中间变量放大到6万多变量，仅R=400、C=76的规模就跑了1500多秒，对于真实场景的7e6行矩阵，这种写法完全不可行。

现有方案的不合理之处 #

• 没有利用线性分式规划的结构特性：直接把除法形式的比值作为目标函数，引入不必要的非线性，还要额外处理分母为0的数值问题。
• 没有利用0-1矩阵的稀疏特性：把逐行逐列的乘加项全部生成Gekko方程，对于百万级行的矩阵，会生成数亿个乘加表达式，求解器根本无法加载。
• 求解器配置错位：用松弛求解整数问题，相当于放弃了整数变量的分支逻辑，自然拿不到非平凡的整数解。

针对性调整方案 #

1. 模型转形，从根源消除非线性和分母为0的问题 #

你的目标是最大化f(x) = (c^T x)/(d^T x)，其中c是A矩阵按列求和得到的常数向量（c_j = sum_i A[i,j]），d是B矩阵按列求和得到的常数向量（d_j = sum_i B[i,j]），x是0-1决策向量，这是典型的线性分式规划问题。
用Charnes-Cooper变换可以把这个问题完全转成线性整数规划，没有除法运算：

• 令t = 1/(d^T x)，y_j = t * x_j
• 变换后目标转化为最大化c^T y
• 配套约束：d^T y =1（强制分母非零），y_j >=0，同时结合x的0-1特性，用大M法把y和x的绑定关系写成线性约束。
• 转形后的问题没有除法运算，从根源上排除了x全0的平凡解，求解难度大幅下降。

2. 预计算矩阵和，解决长方程报错问题 #

所有列和计算都在Python端用numpy完成，和Gekko模型完全解耦：

# 不管R是20还是7e6，只需要计算每列的和即可
c = num.sum(axis=0)
d = den.sum(axis=0)

这一步把原来和行数绑定的长表达式，转成长度仅和列数成正比的短表达式，哪怕R是7e6，c和d也只是长度为500的常数向量，不会出现方程超长的问题。

3. 修正求解器配置 #

• 移除minlp_as_nlp 1参数，转形后的问题是纯线性整数规划，APOPT可以直接处理，不需要做NLP松弛。
• 保留x的整数约束，给t变量设置最小正值下限1e-6，彻底阻断x全0的可行域。

4. 适配大规模场景的性能优化 #

• 不需要把原始矩阵传入Gekko模型，预计算完c和d之后，整个Gekko模型的变量规模仅为C+2，也就是500多个变量，哪怕列数到500，求解时间也不会超过1秒。
• 不要生成逐行的中间变量，否则会把模型规模放大R倍，对于7e6行的规模，直接生成几亿个变量，完全无法求解。

修正后的可运行参考代码 #

import numpy as np
from gekko import GEKKO

C = 43
R = 20
np.random.seed(42)
num = np.random.randint(2, size=(R, C))
den = np.random.randint(2, size=(R, C))

# 预计算列和，与矩阵行数解耦
c = num.sum(axis=0)
d = den.sum(axis=0)

m = GEKKO(remote=False)
# 0-1决策变量
x = m.Array(m.Var, C, lb=0, ub=1, integer=True)
# Charnes-Cooper变换引入的辅助变量
t = m.Var(lb=1e-6, ub=1)
y = m.Array(m.Var, C, lb=0)

# 变换后的线性约束
m.Equation(m.sum(d * y) == 1)
for j in range(C):
    # 绑定y和x的关系
    m.Equation(y[j] >= t - 1 * (1 - x[j]))
    m.Equation(y[j] <= 1 * x[j])

# 线性目标：最大化分子项
m.Obj(-m.sum(c * y))

m.options.SOLVER = 1
m.solve(disp=True)

# 还原最优解
x_val = np.array([xi.value[0] for xi in x]).round().astype(int)
print("最优0-1向量：", x_val)
print("最大比值：", c @ x_val / (d @ x_val + 1e-8))

内容的提问来源于stack exchange，提问作者Rafael Silva