Pyomo中MILP模型max(0,x1-x2)约束的线性化求解求助

要实现y[t] = max(0, x1[t] - x2[t])的线性化约束，需引入二进制辅助变量区分x1[t]-x2[t]的正负状态，具体方案如下：

1. 添加二进制辅助变量 #

定义二进制变量z[t]（每个时间点t对应一个，z[t] ∈ {0,1}），用于标识状态：

• z[t] = 1：x1[t] - x2[t] > 0，此时要求y[t] = x1[t] - x2[t]
• z[t] = 0：x1[t] - x2[t] ≤ 0，此时要求y[t] = 0

2. 确定大M值 #

M需取足够大的正数，确保M大于x1[t]-x2[t]的最大可能值。由于x1、x2为非负变量，可设M为x1的已知上限（若未知则根据场景估算合理值，避免数值溢出）。

3. 构建线性约束 #

在Pyomo中添加以下5组约束，覆盖所有场景：

约束1：限制y的上限（对应z=1时的边界） #

def constraint1(m, t):
    return m.y[t] <= m.x1[t] - m.x2[t] + M*(1 - m.z[t])
model.constraint1 = pyo.Constraint(model.timeindex, rule=constraint1)

• 当z[t]=1时，约束简化为y[t] ≤ x1[t]-x2[t]；当z[t]=0时，右边为x1[t]-x2[t]+M（因x1-x2≤0，右边≥0），非负的y[t]自动满足约束。

约束2：限制y的下限（对应z=1时的等式） #

def constraint2(m, t):
    return m.y[t] >= m.x1[t] - m.x2[t] - M*(1 - m.z[t])
model.constraint2 = pyo.Constraint(model.timeindex, rule=constraint2)

• 当z[t]=1时，约束简化为y[t] ≥ x1[t]-x2[t]，结合约束1可得y[t] = x1[t]-x2[t]；当z[t]=0时，右边≤-M<0，非负的y[t]自动满足约束。

约束3：强制z=0时y=0 #

def constraint3(m, t):
    return m.y[t] <= M * m.z[t]
model.constraint3 = pyo.Constraint(model.timeindex, rule=constraint3)

• 当z[t]=0时，y[t] ≤0，结合y[t]≥0的变量定义，强制y[t]=0；当z[t]=1时，y[t]≤M（符合x1-x2≤M的设定）。

约束4：强制x1-x2>0时z=1 #

def constraint4(m, t):
    return m.x1[t] - m.x2[t] <= M * m.z[t]
model.constraint4 = pyo.Constraint(model.timeindex, rule=constraint4)

• 当x1[t]-x2[t]>0时，M*z[t]必须大于正数，因此z[t]只能取1；当z[t]=0时，强制x1[t]-x2[t]≤0。

约束5：补充边界条件（避免极端情况） #

def constraint5(m, t):
    return m.x1[t] - m.x2[t] >= -M*(1 - m.z[t])
model.constraint5 = pyo.Constraint(model.timeindex, rule=constraint5)

• 当z[t]=1时，x1[t]-x2[t]≥-M（因x1-x2>0自动满足）；当z[t]=0时，x1[t]-x2[t]≥-M（符合大M的边界设定）。

完整代码示例 #

# 假设model已定义，timeindex为时间索引
# 定义二进制变量z
model.z = pyo.Var(model.timeindex, domain=pyo.Binary)

# 设置大M，替换为你的问题中x1的实际上限
U1 = 1000
M = U1

# 添加所有约束
def constraint1(m, t):
    return m.y[t] <= m.x1[t] - m.x2[t] + M*(1 - m.z[t])
model.constraint1 = pyo.Constraint(model.timeindex, rule=constraint1)

def constraint2(m, t):
    return m.y[t] >= m.x1[t] - m.x2[t] - M*(1 - m.z[t])
model.constraint2 = pyo.Constraint(model.timeindex, rule=constraint2)

def constraint3(m, t):
    return m.y[t] <= M * m.z[t]
model.constraint3 = pyo.Constraint(model.timeindex, rule=constraint3)

def constraint4(m, t):
    return m.x1[t] - m.x2[t] <= M * m.z[t]
model.constraint4 = pyo.Constraint(model.timeindex, rule=constraint4)

def constraint5(m, t):
    return m.x1[t] - m.x2[t] >= -M*(1 - m.z[t])
model.constraint5 = pyo.Constraint(model.timeindex, rule=constraint5)

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