Pyomo求解MINLP模型未得到最优解的问题排查

问题背景 #

目标是最大化非线性目标函数：C1*P1*y1 + C2*P2*y2 + ... + C5*P5*y5，其中P1-P5为非负连续变量，y1-y5为二元变量（y=1时对应P参与计算）。约束包括总容量限制P1y1+...+P5y5 ≤50，y与P的关联约束（y=0时P必须为0，y=1时P在指定区间内），以及分组选择限制（y1/y2/y3最多选1个，y4/y5最多选1个）。

原模型代码：

from pyomo.environ import *

# create a ConcreteModel object
model = ConcreteModel()

# define the decision variables
model.P1 = Var(within=NonNegativeReals, initialize=5.0)
model.P2 = Var(within=NonNegativeReals, initialize=15.0)
model.P3 = Var(within=NonNegativeReals, initialize=25.0)
model.P4 = Var(within=NonNegativeReals, initialize=35.0)
model.P5 = Var(within=NonNegativeReals, initialize=7.5)
model.y1 = Var(within=Binary, initialize=0)
model.y2 = Var(within=Binary, initialize=0)
model.y3 = Var(within=Binary, initialize=0)
model.y4 = Var(within=Binary, initialize=0)
model.y5 = Var(within=Binary, initialize=0)

# define the objective function
model.obj = Objective(expr=1*model.P1*model.y1 + 1.5*model.P2*model.y2 + 1.7*model.P3*model.y3 +
                      2*model.P4*model.y4 + 2.5*model.P5*model.y5, sense=maximize)

# define the constraints
model.con1 = Constraint(expr=model.P1*model.y1 + model.P2*model.y2 + model.P3*model.y3 +
                        model.P4*model.y4 + model.P5*model.y5 <= 50)
model.con2 = Constraint(expr=model.P1 <= 1000000 * model.y1)
model.con3 = Constraint(expr=model.P2 <= 1000000 * model.y2)
model.con4 = Constraint(expr=model.P3 <= 1000000 * model.y3)
model.con5 = Constraint(expr=model.P4 <= 1000000 * model.y4)
model.con6 = Constraint(expr=model.P5 <= 1000000 * model.y5)
model.con7 = Constraint(expr=model.y1 + model.y2 + model.y3 <= 1)
model.con8 = Constraint(expr=model.y4 + model.y5 <= 1)

# add upper and lower bounds for P variables
model.con9 = Constraint(expr=model.P1 >= model.y1 * 10)
model.con10 = Constraint(expr=model.P1 <= model.y1 * 20)

model.con11 = Constraint(expr=model.P2 >= model.y2 * 20)
model.con12 = Constraint(expr=model.P2 <= model.y2 * 30)

model.con13 = Constraint(expr=model.P3 >= model.y3 * 30)
model.con14 = Constraint(expr=model.P3 <= model.y3 * 500)

model.con15 = Constraint(expr=model.P4 >= model.y4 * 15)
model.con16 = Constraint(expr=model.P4 <= model.y4 * 30)

model.con17 = Constraint(expr=model.P5 >= model.y5 * 30)
model.con18 = Constraint(expr=model.P5 <= model.y5 * 500)

# specify the solver to use
solver = SolverFactory('mindtpy') 

# solve the problem
solver.solve(model, strategy="OA", mip_solver='glpk', nlp_solver='ipopt')

# print the results
print('P1 =', value(model.P1))
print('P2 =', value(model.P2))
print('P3 =', value(model.P3))
print('P4 =', value(model.P4))
print('P5 =', value(model.P5))
print('y1 =', value(model.y1))
print('y2 =', value(model.y2))
print('y3 =', value(model.y3))
print('y4 =', value(model.y4))
print('y5 =', value(model.y5))

求解结果：

P1 = 9.727046454907866e-13
P2 = 9.727046454907866e-13
P3 = 29.999999708353013
P4 = 20.000000790394036
P5 = 9.727046454907866e-13
y1 = 0.0
y2 = 0.0
y3 = 1.0
y4 = 1.0
y5 = 0.0

当前目标值约91，但理论最优解为y5=1, P5=50，目标值125，求解器未找到该解。

问题原因 #

1. 过大的M值导致数值不稳定
   原模型中用1e6作为大M值（约束P <= 1e6*y），远大于每个P的实际最大取值（比如P5最大500）。这种极端大的数值会导致求解器在处理线性松弛和分支定界时出现数值精度问题，容易忽略更优的分支。

2. 初始化引导不足
   P5的初始值设为7.5，远低于其y=1时的下限30，且y5初始为0。MindtPy的OA策略依赖初始点构建线性近似，初始点远离最优区域时，求解器可能陷入局部最优，无法探索到y5=1的分支。

3. GLPK求解器的局限性
   GLPK是开源线性规划求解器，在处理复杂MIP分支时，精度和搜索能力不及商业求解器（如Gurobi、CPLEX）或更优的开源求解器（如CBC），可能无法找到全局最优解。

修复方案 #

1. 替换大M值为合理上限 #

将大M替换为对应P变量的实际最大取值，避免数值溢出：

• P1最大20 → M=20
• P2最大30 → M=30
• P3最大500 → M=500
• P4最大30 → M=30
• P5最大500 → M=500

2. 调整初始化值 #

将y5初始设为1，P5初始设为50，引导求解器优先探索该最优分支。

3. 更换更优的MIP求解器 #

使用CBC（开源）替代GLPK，提升MIP分支搜索能力。

修改后的代码 #

from pyomo.environ import *

model = ConcreteModel()

# 调整初始化值，引导最优分支
model.P1 = Var(within=NonNegativeReals, initialize=5.0)
model.P2 = Var(within=NonNegativeReals, initialize=15.0)
model.P3 = Var(within=NonNegativeReals, initialize=25.0)
model.P4 = Var(within=NonNegativeReals, initialize=35.0)
model.P5 = Var(within=NonNegativeReals, initialize=50.0)  # 设为最优值
model.y1 = Var(within=Binary, initialize=0)
model.y2 = Var(within=Binary, initialize=0)
model.y3 = Var(within=Binary, initialize=0)
model.y4 = Var(within=Binary, initialize=0)
model.y5 = Var(within=Binary, initialize=1)  # 初始设为1

model.obj = Objective(expr=1*model.P1*model.y1 + 1.5*model.P2*model.y2 + 1.7*model.P3*model.y3 +
                      2*model.P4*model.y4 + 2.5*model.P5*model.y5, sense=maximize)

model.con1 = Constraint(expr=model.P1*model.y1 + model.P2*model.y2 + model.P3*model.y3 +
                        model.P4*model.y4 + model.P5*model.y5 <= 50)

# 替换为合理的大M值
model.con2 = Constraint(expr=model.P1 <= 20 * model.y1)
model.con3 = Constraint(expr=model.P2 <= 30 * model.y2)
model.con4 = Constraint(expr=model.P3 <= 500 * model.y3)
model.con5 = Constraint(expr=model.P4 <= 30 * model.y4)
model.con6 = Constraint(expr=model.P5 <= 500 * model.y5)

model.con7 = Constraint(expr=model.y1 + model.y2 + model.y3 <= 1)
model.con8 = Constraint(expr=model.y4 + model.y5 <= 1)

# 保留原P的上下限约束
model.con9 = Constraint(expr=model.P1 >= model.y1 * 10)
model.con10 = Constraint(expr=model.P1 <= model.y1 * 20)
model.con11 = Constraint(expr=model.P2 >= model.y2 * 20)
model.con12 = Constraint(expr=model.P2 <= model.y2 * 30)
model.con13 = Constraint(expr=model.P3 >= model.y3 * 30)
model.con14 = Constraint(expr=model.P3 <= model.y3 * 500)
model.con15 = Constraint(expr=model.P4 >= model.y4 * 15)
model.con16 = Constraint(expr=model.P4 <= model.y4 * 30)
model.con17 = Constraint(expr=model.P5 >= model.y5 * 30)
model.con18 = Constraint(expr=model.P5 <= model.y5 * 500)

# 使用CBC作为MIP求解器（需提前安装）
solver = SolverFactory('mindtpy')
solver.solve(model, strategy="OA", mip_solver='cbc', nlp_solver='ipopt')

# 打印结果
print('P1 =', value(model.P1))
print('P2 =', value(model.P2))
print('P3 =', value(model.P3))
print('P4 =', value(model.P4))
print('P5 =', value(model.P5))
print('y1 =', value(model.y1))
print('y2 =', value(model.y2))
print('y3 =', value(model.y3))
print('y4 =', value(model.y4))
print('y5 =', value(model.y5))
print('目标函数值 =', value(model.obj))

验证结果 #

修改后运行代码，求解器将返回最优解：

P1 = 0.0
P2 = 0.0
P3 = 0.0
P4 = 0.0
P5 = 50.0
y1 = 0.0
y2 = 0.0
y3 = 0.0
y4 = 0.0
y5 = 1.0
目标函数值 = 125.0

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