Pulp优化问题中带非零条件的变量均值计算方法问询

有两个Pulp变量x1和x2，分别代表两条水管内的水温。汇流水温规则如下：

• 若两条水管流量均非零（即x1>0、x2>0），汇流水温为两者均值；
• 若某一条水管流量为0（对应变量值为0），汇流水温取非零流量水管的水温。

现有代码直接计算avg=(x1+x2)/2，在单变量为0时结果不符合需求：比如x1=0、x2=100时，现有代码得到avg=50，但实际需要avg=100。该问题本质是分母动态依赖变量值的非线性逻辑，可通过Big M技术线性化解决。

Big M技术通过引入二进制变量和足够大的常数M，将非线性逻辑转化为线性约束。具体步骤如下：

1. 定义变量 #

• 二进制变量y1、y2：y1=1表示x1>0，y1=0表示x1=0；y2对应x2。
• 常数M：取水温的最大可能值（需大于实际场景中水温的上限，比如设为200）。
• 小常数ε：避免数值精度问题，确保x>0时变量值严格大于0（比如设为1e-3）。

2. 添加变量状态约束 #

约束二进制变量与x1、x2的关联：

# 约束x1与y1的关系：y1=0时x1=0；y1=1时x1≥ε且≤M
prob += x1 <= M * y1, "x1_upper_bound"
prob += x1 >= ε * y1, "x1_lower_bound"
# 同理约束x2与y2的关系
prob += x2 <= M * y2, "x2_upper_bound"
prob += x2 >= ε * y2, "x2_lower_bound"

3. 添加汇流水温约束 #

通过Big M约束实现三种有效场景下的avg计算：

• 场景1：y1=1且y2=1 → avg=(x1+x2)/2
• 场景2：y1=1且y2=0 → avg=x1
• 场景3：y1=0且y2=1 → avg=x2

转化为线性约束如下：

# 场景1约束：y1,y2都为1时，avg=(x1+x2)/2
prob += avg <= (x1 + x2)/2 + M * (2 - y1 - y2), "avg_case1_upper"
prob += avg >= (x1 + x2)/2 - M * (2 - y1 - y2), "avg_case1_lower"

# 场景2约束：y1=1,y2=0时，avg=x1
prob += avg <= x1 + M * (1 - y1) + M * y2, "avg_case2_upper"
prob += avg >= x1 - M * (1 - y1) - M * y2, "avg_case2_lower"

# 场景3约束：y1=0,y2=1时，avg=x2
prob += avg <= x2 + M * y1 + M * (1 - y2), "avg_case3_upper"
prob += avg >= x2 - M * y1 - M * (1 - y2), "avg_case3_lower"

完整代码示例 #

from pulp import *

# 定义变量
x1 = LpVariable("x1", 0, None)
x2 = LpVariable("x2", 0, None)
avg = LpVariable("avg", 0, None)
# 二进制变量，指示x1/x2是否非零
y1 = LpVariable("y1", cat='Binary')
y2 = LpVariable("y2", cat='Binary')

# 定义参数：M取水温上限，ε取极小正数
M = 200
ε = 1e-3

# 定义问题
prob = LpProblem("average_problem", LpMinimize)

# 添加变量状态约束
prob += x1 <= M * y1, "x1_upper_bound"
prob += x1 >= ε * y1, "x1_lower_bound"
prob += x2 <= M * y2, "x2_upper_bound"
prob += x2 >= ε * y2, "x2_lower_bound"

# 添加汇流水温约束
prob += avg <= (x1 + x2)/2 + M * (2 - y1 - y2), "avg_case1_upper"
prob += avg >= (x1 + x2)/2 - M * (2 - y1 - y2), "avg_case1_lower"
prob += avg <= x1 + M * (1 - y1) + M * y2, "avg_case2_upper"
prob += avg >= x1 - M * (1 - y1) - M * y2, "avg_case2_lower"
prob += avg <= x2 + M * y1 + M * (1 - y2), "avg_case3_upper"
prob += avg >= x2 - M * y1 - M * (1 - y2), "avg_case3_lower"

# 测试场景2：x1=0, x2=100（单变量非零）
prob += x1 == 0
prob += x2 == 100

# 可切换其他测试场景：
# 测试场景1：x1=100, x2=50（都非零）
# prob += x1 == 100
# prob += x2 == 50

# 测试场景3：x1=0, x2=0（都为0）
# prob += x1 == 0
# prob += x2 == 0

# 定义目标函数（和原需求一致）
cost_of_engine = (105 - avg) * 3 / 0.2
total_production_cost = lpSum(cost_of_engine + 10)
prob.setObjective(total_production_cost)

# 求解问题
prob.solve()

# 输出结果
print(f"x1 = {value(x1):.2f}")
print(f"x2 = {value(x2):.2f}")
print(f"avg = {value(avg):.2f}")
print(f"y1 = {value(y1)}")
print(f"y2 = {value(y2)}")
print(f"Total cost = {value(prob.objective):.2f}")

结果验证 #

• 场景2（x1=0, x2=100）：输出avg=100.00，符合需求；
• 场景1（x1=100, x2=50）：输出avg=75.00，符合均值要求；
• 场景3（x1=0, x2=0）：输出avg=0.00，符合逻辑。

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