Pyomo MINLP 求解器没有选择最佳结果

Question

我正在使用 Pyomo 来解决 MINLP 问题。它似乎运行没有问题，但它选择的解决方案似乎并不是最佳的。

问题的要点是我正在尝试优化形式为 (C1 * P1 * y1) + (C2 * P2 * y2)...(C5 * P5 * y5)= Z 的函数，其中 C1。 ..C5 是常量，P1...P5 和 y1...y5 是决策变量。 P 变量是连续的，y 变量是二元的，其中 y 变量用作 P 变量的选择器（即，如果 y5 为 1，则 P5 被分解为目标等）。我们还得到一个约束，即 P1y1 + ... P5y5 <= 50.

的总和

代码如下：

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

请注意，求解器已打开 y3 和 y4，并给出了 P3 和 P4 值。如果我们进行数学计算，目标函数将等于 91 左右。但是，最优解应该是 y5 = 1 和 P5 = 50，以产生 125 的目标函数值。

这里有什么问题吗？

Answer 1

约束条件

con2

、

con3

、

con4

、

con5

和

con6

未按预期工作。尝试将这些约束替换为：

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)

这确保连续变量（

P1

、

P2

、

P3

、

P4

、

P5

）在相应的二进制变量（

y1

、

y2

、

y3）有效地“关闭”时

,

y4

,

y5

) 等于零。在您的原始约束（

con2

到

con6

）中，您将连续变量（

Pi

）乘以二元变量（

yi

），这并不能有效地强制执行所需的行为。原始约束仅确保如果

yi == 1

，则

Pi

可以在其范围内取任何值，但它们不会强制

Pi

在

yi == 0

时为零。

---

编辑 #1： 您可以向模型添加额外的约束，以帮助引导求解器找到更好的解决方案。您提到最优解应该有

y4 == 1

和

P4 == 50

，因此我们可以设置一个额外的约束条件，强制至少一个

y

变量等于

1

：

model.con11 = Constraint(expr=model.y1 + model.y2 + model.y3 + model.y4 + model.y5 >= 1)

---

编辑#2：您可以增加求解器的最大迭代次数，这将使其有更多时间找到更好的解决方案：

solver.solve(model, strategy='OA', mip_solver='glpk', nlp_solver='ipopt', max_iter=5000)

另一个值得考虑的选择是创建一个额外的约束来防止当前的次优解决方案。例如，如果当前解是

y3 == 1

和

y4 == 1

与

P3 == 29.999

和

P4 == 20

，您可以添加一个约束来禁止这种组合并强制求解器搜索其他解：

model.con12 = Constraint(expr=model.y3 + model.y4 <= 1)

最后，如果仍然不成功，也许看看其他求解器如

bonmin

、

couenne

或

baron

是否会更方便。

Answer 2

在编写模型时，编辑后，约束 2-6 是多余的。当

P

为零时，其他后来的约束将

y

变量限制为零。

因此

P

变量在此模型中是半连续的，这意味着它们要么为零，要么在一个范围内。如果这是您的愿望，可以通过从目标函数中删除

y

变量来简单地使模型线性化，因为它们也是多余的。如果

y

为零，则相关的

P

也为零，从而为您留下一个目标函数，即

sum(P * weight)

.

最后，我没有在您的代码中看到您检查解决方案对象以确保状态为“最佳”的位置。这是一个巨大的错误。除非求解器声称最优，否则结果是乱码。