Pyomo中如何给不等式上下界乘以二进制变量控制约束生效

原代码问题说明 #

你当前写的约束规则逻辑错误：约束规则是在模型构建阶段执行的，此时二进制变量还没有求解结果，无法通过if modelito.p==0这类判断动态跳过约束，该写法完全不符合混合整数规划的约束构建逻辑。

实现原理 #

你需要用大M法实现开关约束，逻辑如下：

• 定义一个足够大的常数M，取值要大于所有决策变量的最大可能取值（你这个场景选1000即可）
• 将上下界约束拆分为两个独立约束，通过二进制变量控制约束是否生效：
  • 当二进制变量为1时，约束为min <= 决策变量 <= max，正常生效
  • 当二进制变量为0时，约束变为0 <= 决策变量 <= M，永远成立，相当于自动失效
• 如果你需要二进制变量为0时决策变量必须取0，直接去掉M项即可，约束写为min*b <= 变量 <= max*b即可。

修改后的代码 #

from pyomo.environ import *

m=4
n=3
# 定义大M常数，大于所有变量的最大可能值
M = 1000

#Boiler minimo
boiler_min= [400, 500, 300, 240]
#boiler máximo
boiler_max= [900,700,600,800]
#boiler cost per ton
boiler_ton= [9,7,8,6]
#boiler fixed cost
boiler_fc= [160,200,190,250]

#turbine minimo
turbine_min= [100,400,300]
#turbine máximo
turbine_max= [700,800,600]
#turbine Kwh per ton of stream
turbine_ton= [7,3,6]
#turbine Processing Cost per Ton
turbine_pc = [9,4,6]

modelito= ConcreteModel()
modelito.m = RangeSet(m)
modelito.n = RangeSet(n)

modelito.b = Var(modelito.m, domain=NonNegativeIntegers)
modelito.t = Var(modelito.n, domain=NonNegativeIntegers)

#Variables binarias
modelito.p = Var(modelito.m, domain=Binary)  #binary boiler        
modelito.q= Var(modelito.n, domain=Binary)   #binary turbine

#declare objective
modelito.cost= Objective(expr = (sum(boiler_ton[i-1]*modelito.b[i] for i in  modelito.m) + sum(turbine_pc[j-1]*modelito.t[j] for j in modelito.n) + sum(boiler_fc[i-1]*modelito.p[i] for i in  modelito.m )), sense=minimize)

# 锅炉约束：p[i]=1时b[i]在上下限范围内，p[i]=0时约束失效
def boiler_low_rule(modelito, i):
    return modelito.b[i] >= boiler_min[i-1] * modelito.p[i]
modelito.boiler_low = Constraint(modelito.m, rule=boiler_low_rule)

def boiler_high_rule(modelito, i):
    return modelito.b[i] <= boiler_max[i-1] * modelito.p[i] + M * (1 - modelito.p[i])
modelito.boiler_high = Constraint(modelito.m, rule=boiler_high_rule)

# 汽轮机约束：q[j]=1时t[j]在上下限范围内，q[j]=0时约束失效
def turbine_low_rule(modelito, j):
    return modelito.t[j] >= turbine_min[j-1] * modelito.q[j]
modelito.turbine_low = Constraint(modelito.n, rule=turbine_low_rule)

def turbine_high_rule(modelito, j):
    return modelito.t[j] <= turbine_max[j-1] * modelito.q[j] + M * (1 - modelito.q[j])
modelito.turbine_high = Constraint(modelito.n, rule=turbine_high_rule)

# 发电量约束
modelito.constraint2 = Constraint(expr=sum(turbine_ton[i-1]*modelito.t[i] for i in modelito.n)>=9000)

modelito.pprint()

可选调整 #

如果要求锅炉/汽轮机关闭时（p[i]/q[j]=0）出力必须为0，直接删除上下限约束中的+ M * (1 - 变量)项即可，此时p[i]=0时约束会自动限定b[i]=0。

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