Gekko求解非线性规划需预定义整数变量值才能找到解

问题背景 #

我用Gekko求解一个模拟资金在交易所流转的非线性规划问题，每个交易所启用需支付固定费用。测试网络是1→2→3（箭头代表使用的交易所），用二进制变量Q标记交易所是否启用（1为启用，0为禁用）。由于资金必须从1流向3，两个交易所必须启用，但Gekko无法自动识别这一约束，直接求解会提示无解。只有手动将二进制变量q121、q231初始值设为1时，才能得到正确解：

• i121=10，q121=1
• i231=5，o231≈3.3333，q231=1
• 目标函数值≈-2.3333

但复杂网络中需要Gekko自主判断是否启用交易所，因此需要找到无需手动初始化的解决方案。

原代码如下：

from gekko import GEKKO

m = GEKKO(remote=False)
m.options.SOLVER = 1

# Constants
vol = m.Const(10, 'Vol')
x = m.Const(10, 'x')
y = m.Const(10, 'y')
k = m.Const(100, 'k')
g = m.Const(0.5, 'gas_fee')

# Inputs and Outputs
i121, i231, o231 = [m.Var(lb=0) for i in range(3)]

# Activation
q121, q231 = [m.Var(lb=0, ub=1, integer=True) for i in range(2)]

# Constraints
m.Equation((y - (k/(x+i121))) == i231)
m.Equation((y - (k/(x+i231))) == o231)

m.Equation(i121 == i121*q121)
m.Equation(i231 == i231*q231)

m.Equation(i121 == vol)

# Objective Function
m.Minimize(g*q121 + g*q231 - o231)
m.solve(disp=True)

手动初始化二进制变量的代码：

q121, q231 = [m.Var(value = 1, lb=0, ub=1, integer=True) for i in range(2)]

问题分析 #

当前约束i121 == i121*q121和i231 == i231*q231存在逻辑漏洞：

• 当q=0时，i必须为0；但当i>0时，q必须为1的逻辑没有被严格约束。Gekko的混合整数非线性规划求解器可能陷入局部最优，或无法识别隐含的q=1约束。
• 测试场景中i121=vol=10>0，因此q121必须为1；而i231由i121计算得出也大于0，因此q231必须为1，但求解器无法自动推导这层关系。

解决方案 #

方法1：添加显式逻辑约束（大M法） #

引入足够大的常数M（大M法），强制当流量大于0时对应的二进制变量必须为1：

• 对于每个流量变量i，添加约束i <= M * q，其中M是该流量变量的最大可能取值。当q=0时，i必须为0；当q=1时，i可取值范围不受限（只要符合其他约束）。

修改后的代码 #

from gekko import GEKKO

m = GEKKO(remote=False)
m.options.SOLVER = 1

# Constants
vol = m.Const(10, 'Vol')
x = m.Const(10, 'x')
y = m.Const(10, 'y')
k = m.Const(100, 'k')
g = m.Const(0.5, 'gas_fee')
# 大M常数：取对应流量变量的最大可能值
M1 = vol.value  # i121的最大值是vol=10
M2 = y.value - k.value/(x.value + vol.value)  # i231的最大值是5

# Inputs and Outputs
i121, i231, o231 = [m.Var(lb=0) for i in range(3)]

# Activation
q121, q231 = [m.Var(lb=0, ub=1, integer=True) for i in range(2)]

# Constraints
m.Equation((y - (k/(x+i121))) == i231)
m.Equation((y - (k/(x+i231))) == o231)

# 改进的逻辑约束：强制流量存在时q必须为1
m.Equation(i121 <= M1 * q121)
m.Equation(i231 <= M2 * q231)

m.Equation(i121 == vol)

# Objective Function
m.Minimize(g*q121 + g*q231 - o231)
m.solve(disp=True)

# 打印结果
print(f"i121: {i121.value[0]}, q121: {q121.value[0]}")
print(f"i231: {i231.value[0]}, o231: {o231.value[0]}, q231: {q231.value[0]}")
print(f"目标函数值: {m.options.objfcnval}")

说明 #

• 大M值的选择需合理：既不能过大导致数值不稳定，也不能过小无法覆盖变量的最大可能取值。复杂网络中可通过动态计算每个流量变量的理论最大值来设置M。
• 修改后的约束能让Gekko明确识别流量与二进制变量的绑定关系，无需手动初始化即可找到正确解。

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