基于BigM约束线性化Pyomo模型的问题与修正方案

我已经找到一个可行解法，替代了之前的model.bigM_cons约束，现在模型能实现合理线性化并返回可行解，作为新手希望得到大家的反馈。

我新增了变量model.new来表示model.auction*model.volume的乘积，然后为这个新变量添加三个约束，并将其用于目标函数中：

新BigM约束实现 #

model.m = Param(initialize = 100)
model.new = Var(model.WEEK_PROD_flat, within = NonNegativeIntegers)

# 当auction为0时强制new为0；auction为1时new介于0和M之间
def c1_rule(model,w,p):
    return model.new[w,p] <= model.auction[w,p]*model.m
model.c1 = Constraint(model.WEEK_PROD_flat, rule = c1_rule)

# 强制new始终不超过volume
def c2_rule(model,w,p):
    return(model.new[w,p] <= model.volume[w,p])
model.c2 = Constraint(model.WEEK_PROD_flat, rule = c2_rule)

# 当auction为1时强制new等于volume
def c3_rule(model,w,p):
    return model.new[w,p] >= model.volume[w,p]-(1-model.auction[w,p])*model.m
model.c3 = Constraint(model.WEEK_PROD_flat, rule = c3_rule)

原问题描述 #

我正在构建一个优化模型，包含两个以稀疏(周,产品)索引集model.WEEK_PROD_flat为索引的决策变量：

• model.volume：非负整数变量，代表待售产品数量；
• model.auction：二进制变量，标记特定(周,产品)组合是否开展销售。

模型原本存在非线性问题，因为多个约束和目标函数中包含model.auction * model.volume的乘积，model.auction的作用是控制model.volume是否启用。

我知道需要用BigM约束来线性化这个乘积，于是尝试构建约束，要求：当model.auction为0时，model.volume必须为0；仅当model.auction为1时，model.volume可以赋值。

我最初的约束实现如下：

model.bigM = Param(initialize = 100) # 定义BigM值

def bigM_rule(model,w,p):
    return model.volume[w,p] <= model.auction[w,p] * model.bigM

model.bigM_cons = Constraint(model.WEEK_PROD_flat, rule = bigM_rule)  

这个约束在简单的最小工作示例(MWE)中能正常运行，但加入完整模型后返回了不可行解。需要说明的是，完整模型使用非线性约束和目标函数时能得到可行解（尽管不一定是全局最优），只有应用这个BigM约束时才会不可行，所以我怀疑约束公式有问题，希望大家能指出可能的错误。

完整最小工作示例(MWE) #

weekly_products = {
    1: ['Q24'],
    2: ['Q24', 'J24'],
    3: ['Q24', 'J24','F24'],
    4: ['J24', 'F24'],
    5: ['F24']
}

product_weeks = {'Q24': [1, 2, 3],
                 'J24': [2, 3, 4],
                 'F24': [3, 4, 5]}

prices = {(1, 'Q24'):43.42,
      (2, 'Q24'):43.73,
      (2, 'J24'):24.89,
      (3, 'Q24'):44.03,
      (3, 'J24'):25.54,
      (3, 'F24'):43.10,
      (4, 'J24'):26.15,
      (4, 'F24'):43.45,
      (5, 'F24'):43.77}

from pyomo.environ import *
model = ConcreteModel()

# 定义集合
model.WEEKS = Set(initialize = [1,2,3,4,5])
model.PRODS = Set(initialize = ['Q24','J24','F24'])
model.WEEK_PROD = Set(model.WEEKS, initialize=weekly_products)
model.WEEK_PROD_flat = Set(initialize=[(w, p) for w in model.WEEKS for p in model.WEEK_PROD[w]])
model.PROD_WEEK = Set(model.PRODS, initialize = product_weeks)

# 定义变量
model.volume = Var(model.WEEK_PROD_flat, within = NonNegativeIntegers, bounds = (0,60))
model.auction = Var(model.WEEK_PROD_flat, within = Binary)

# 定义参数
model.price = Param(model.WEEK_PROD_flat, initialize = prices)
model.weekMax = Param(initialize = 1)
model.prodMax = Param(initialize = 3)
model.bigM = Param(initialize = 100)

# 定义约束
def weekMax_rule(model,i):
    return sum(model.auction[i,j] for j in model.WEEK_PROD[i]) <= model.weekMax
model.weekMax_const = Constraint(model.WEEKS, rule = weekMax_rule)

def prodMax_rule(model,j):
    return sum(model.auction[i,j] for i in model.PROD_WEEK[j]) <=model.prodMax
model.prodMax_const = Constraint(model.PRODS, rule = prodMax_rule)

def bigM_rule(model,w,p):
    return model.volume[w,p] <= model.auction[w,p] * model.bigM
model.bigM_cons= Constraint(model.WEEK_PROD_flat, rule = bigM_rule)

# 目标函数
def objective_rule(model):
    return sum(model.volume[w,p] * model.price[w,p] for p in model.PRODS for w in model.PROD_WEEK[p])
model.maximiseRev = Objective(rule = objective_rule, sense = maximize)

optimizer = SolverFactory('scip')
results = optimizer.solve(model)
model.display()

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