Pyomo MINLP解算器没有选择最优结果

我使用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。

怎么了?

当您编写模型时，在编辑之后，约束2-6是多余的。当y为零时，其他后续约束将P变量箝位为零。

因此，P变量在该模型中是半连续的，这意味着它们要么为零，要么在一个范围内。如果这是您的愿望，则可以通过从目标函数中删除y变量来简单地使模型线性化，因为它们也是多余的。如果y为零，那么相关的P也为零，留给您的目标函数是sum(P * weight)。

最后，我在你的代码中没有看到你在检查解决方案对象以确保状态是"OPTIMAL"。这是一个巨大的错误。除非解算器声称是最优的，否则结果是乱码。

编辑:添加代码。

从上次迭代…

您仍然在约束1中无缘无故地使用了多个变量。
不需要在这个小模型上初始化变量，特别是当初始化无效时
打印模型并检查，结果与所示相同
如果这个变大，我会索引P和y，而不是单独手动阻塞它们。
如果你正在制作线性模型，请使用线性求解器…

from pyomo.environ import *
# create a ConcreteModel object
model = ConcreteModel()
# define the decision variables
model.P1 = Var(within=NonNegativeReals)
model.P2 = Var(within=NonNegativeReals)
model.P3 = Var(within=NonNegativeReals)
model.P4 = Var(within=NonNegativeReals)
model.P5 = Var(within=NonNegativeReals)
model.y1 = Var(within=Binary)
model.y2 = Var(within=Binary)
model.y3 = Var(within=Binary)
model.y4 = Var(within=Binary)
model.y5 = Var(within=Binary)
# define the objective function
model.obj = Objective(expr=1*model.P1 + 1.5*model.P2 + 1.7*model.P3 +
2*model.P4 + 2.5*model.P5, sense=maximize)
# define the constraints
model.con1 = Constraint(expr=model.P1 + model.P2 + model.P3 +
model.P4 + model.P5 <= 50)
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)

# QA the model:
model.pprint()
# specify the solver to use
solver = SolverFactory('cbc') 
# solve the problem
result = solver.solve(model) # strategy="OA", mip_solver='glpk', nlp_solver='ipopt', tee=True, nlp_solver_tee=True, mip_solver_tee=True)
# review the solution and solve status:
print(result)

# 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))

约束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可以取其边界内的任何值，但是当yi == 0.

时，它们不强制Pi为零。

编辑# 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会更方便。

编辑:添加代码。

相关内容

最新更新

热门标签：