Pyomo中含max/min的Pavg与H_ad约束编写问题求助

错误代码与异常信息 #

尝试直接用Python内置max()函数编写Pavg的约束时触发异常：

# 错误的Pavg约束实现
def Avg_pressure(model, I, J, T, S, r):
    return model.Pavg[I, J, T, S] == max(model.A1p[r] * model.Pin[I, J, T, S] + model.A2p[r] * model.Pout[I, J, T, S] + model.Bp[r] for r in model.R1_2)
model.constraint52 = Constraint(model.Pipe, model.T, model.S, model.R1_2, rule=Avg_pressure)

异常提示：

PyomoException: Cannot convert non-constant Pyomo expression (0.54886306Pin[2,3,1,1] + 0.44727615Pout[2,3,1,1] + 3.02961114e-05 < 0.62043087Pin[2,3,1,1] + 0.31263309Pout[2,3,1,1] + 0.00039507888) to bool.
This error is usually caused by using a Var, unit, or mutable Param in a
Boolean context such as an "if" statement, or when checking container
membership or equality. For example,
m.x = Var()
if m.x >= 1:
pass
and
m.y = Var()
if m.y in [m.x, m.y]:
pass
would both cause this exception.

问题原因 #

Pyomo的变量/表达式是符号化的，无法直接用Python内置的max()/min()函数进行比较（这类函数会尝试将表达式转换为布尔值判断大小）。必须将max/min逻辑转化为一组等价的数学约束。

正确的约束实现方式 #

1. Pavg = max(...) 约束实现 #

要实现Pavg[I,J,T,S] = max_{r∈R1_2} (A1p[r]·Pin[I,J,T,S] + A2p[r]·Pout[I,J,T,S] + Bp[r])，需拆分为以下约束：

步骤1：定义二进制选择变量

用于标记哪个表达式取到最大值：

model.z_p = Var(model.Pipe, model.T, model.S, model.R1_2, domain=Binary)

步骤2：Pavg大于等于所有表达式（已实现部分）

def Avg_pressure_lower(model, I, J, T, S, r):
    return model.Pavg[I, J, T, S] >= model.A1p[r] * model.Pin[I, J, T, S] + model.A2p[r] * model.Pout[I, J, T, S] + model.Bp[r]
model.constraint51 = Constraint(model.Pipe, model.T, model.S, model.R1_2, rule=Avg_pressure_lower)

步骤3：Pavg小于等于选中的表达式（引入大M常数）

M_p需根据实际数据设置为足够大的常数（大于表达式的最大可能差值）：

M_p = 1e6  # 根据你的模型数据调整

def Avg_pressure_upper(model, I, J, T, S, r):
    return model.Pavg[I, J, T, S] <= model.A1p[r] * model.Pin[I, J, T, S] + model.A2p[r] * model.Pout[I, J, T, S] + model.Bp[r] + M_p * (1 - model.z_p[I, J, T, S, r])
model.constraint52 = Constraint(model.Pipe, model.T, model.S, model.R1_2, rule=Avg_pressure_upper)

步骤4：确保仅一个表达式被选中

def Avg_pressure_z_sum(model, I, J, T, S):
    return sum(model.z_p[I, J, T, S, r] for r in model.R1_2) == 1
model.constraint53 = Constraint(model.Pipe, model.T, model.S, rule=Avg_pressure_z_sum)

2. h_ad = min(...) 约束实现 #

同理，实现h_ad[I,J,T,S] = min_{r∈R1_2} (Ah[I,J,r]·Pout[I,J,T,S] + Bh[I,J,r])：

步骤1：定义二进制选择变量

model.z_h = Var(model.C, model.T, model.S, model.R1_2, domain=Binary)

步骤2：h_ad小于等于所有表达式（已实现部分）

def H_ad_upper(model, I, J, T, S, r):
    return model.h_ad[I, J, T, S] <= model.Ah[I, J, r] * model.Pout[I, J, T, S] + model.Bh[I, J, r]
model.constraint55 = Constraint(model.C, model.T, model.S, model.R1_2, rule=H_ad_upper)

步骤3：h_ad大于等于选中的表达式（引入大M常数）

M_h = 1e6  # 根据实际数据调整

def H_ad_lower(model, I, J, T, S, r):
    return model.h_ad[I, J, T, S] >= model.Ah[I, J, r] * model.Pout[I, J, T, S] + model.Bh[I, J, r] - M_h * (1 - model.z_h[I, J, T, S, r])
model.constraint56 = Constraint(model.C, model.T, model.S, model.R1_2, rule=H_ad_lower)

步骤4：确保仅一个表达式被选中

def H_ad_z_sum(model, I, J, T, S):
    return sum(model.z_h[I, J, T, S, r] for r in model.R1_2) == 1
model.constraint57 = Constraint(model.C, model.T, model.S, rule=H_ad_z_sum)

简化方案（无需二进制变量） #

如果Pavg是目标函数的最大化对象，只需保留Pavg >= 所有表达式的约束，然后在目标函数中最大化Pavg，求解器会自动让Pavg等于所有表达式的最大值。同理，如果h_ad是目标函数的最小化对象，只需保留h_ad <= 所有表达式的约束，然后最小化h_ad即可得到等于min的结果。

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