Pyomo中变量比较与二进制变量的条件建模：Big-M方法详解_技术教程

本文深入探讨了在pyomo优化模型中如何正确处理涉及决策变量的条件逻辑。当尝试直接在python的`if`语句中使用pyomo变量进行比较时，会导致`typeerror`。针对此问题，本文详细介绍了如何利用big-m方法和辅助二进制变量，将复杂的条件逻辑转化为线性约束，从而在pyomo中实现变量间的有效比较与状态控制，并提供了具体的代码示例与原理分析。

在构建Pyomo优化模型时，经常需要根据模型中的决策变量（Var）的状态来触发或控制其他变量或约束。例如，当某个生产量达到一定阈值时，一个二进制变量需要被激活。然而，直接在Python的if语句中对Pyomo的Var对象进行比较，会导致TypeError: Relational expression used in an unexpected Boolean context。这是因为Pyomo变量是符号表达式，它们的值在模型求解之前是未知的，不能像普通Python数值一样直接进行布尔判断。

理解问题根源

考虑以下Pyomo约束定义中的条件语句：

def gen3_on_off(model, m):
    # 这行代码会导致错误，因为model.gen1_use[m]和model.gen2_use[m]是Pyomo变量
    if model.gen1_use[m] + model.gen2_use[m] <= 0.90 * model.load_profile[m]:
        return model.gen3_status[m] == 1
    else:
        return model.gen3_status[m] == 0

当Pyomo尝试构建gen3_on_off约束时，它会遍历model.m_index中的每个索引m。在每次迭代中，Python解释器会尝试评估if语句中的条件表达式。然而，model.gen1_use[m]和model.gen2_use[m]是pyomo.environ.Var类型的对象，它们代表了优化问题中的未知决策变量。对这些变量执行

解决方案：Big-M方法

解决这类问题的标准方法是使用Big-M（大M）方法，将条件逻辑转化为一组线性不等式。这通常需要引入一个辅助的二进制变量来表示条件的真假。

假设我们希望实现这样的逻辑： 如果 X >= T，则 Y = 1；否则，Y = 0。 其中 X 是一个或多个Pyomo变量的组合，T 是一个阈值，Y 是一个二进制变量（model.gen3_status[m]）。

我们可以通过以下两个Big-M约束来实现这一条件逻辑：

约束1：当 Y=1 时，强制 X >= T + epsilonX >= T + epsilon - BigM * (1 - Y)
- 如果 Y = 1，则 X >= T + epsilon。这确保了只有当 X 严格大于 T 时，Y 才能为1。
- 如果 Y = 0，则 X >= T + epsilon - BigM。由于 BigM 是一个足够大的正数，这个约束在 Y=0 时几乎总是满足的，不会限制 X 的下限。
约束2：当 Y=0 时，强制 X X

参数说明：

BigM：一个足够大的正数。它必须足够大，以至于当条件不激活时，它不会错误地限制变量的范围。然而，过大的BigM可能导致数值不稳定性，因此应选择一个合理的值，通常是问题中变量可能的最大范围。
epsilon (eps)：一个非常小的正数，用于处理严格不等式（> 或 = T + epsilon 确保 X 必须严格大于 T 才能满足条件。如果条件是 X >= T（包含等于），则可以省略 epsilon。

Pyomo中的实现

将上述Big-M方法应用于原问题，假设我们希望当 model.gen1_use[m] + model.gen2_use[m] 大于等于 0.30 * model.load_profile[m] 时，model.gen3_status[m] 为1，否则为0。

首先定义 eps 和 bigm 常量：

# 一些常量
eps = 1e-3  # 用于创建间隙，处理严格不等式
bigm = 1e3  # 选择一个足够大但不过大的值，避免数值不稳定性

然后，将原有的 gen3_on_off 约束替换为两个新的Big-M约束：

# 创建2个Big-M约束
def gen3_on_off_rule1(model, m):
    # 逻辑：如果 gen1_use + gen2_use >= 0.30 * load_profile + eps，则 gen3_status 倾向于 1
    return model.gen1_use[m] + model.gen2_use[m] >= 0.30 * model.load_profile[m] + eps - bigm * (1 - model.gen3_status[m])

def gen3_on_off_rule2(model, m):
    # 逻辑：如果 gen1_use + gen2_use <= 0.30 * load_profile，则 gen3_status 倾向于 0
    return model.gen1_use[m] + model.gen2_use[m] <= 0.30 * model.load_profile[m] + bigm * model.gen3_status[m]

model.gen3_on_off_bigm1 = Constraint(model.m_index, rule=gen3_on_off_rule1)
model.gen3_on_off_bigm2 = Constraint(model.m_index, rule=gen3_on_off_rule2)

代码整合示例：

from pyomo.environ import *
import numpy as np
import pandas as pd

# 创建一个具体的模型
model = ConcreteModel()
idx = 20
np.random.seed(idx)
model.m_index = Set(initialize=list(range(idx)))

model.load_profile = Param(model.m_index, initialize=dict(zip(model.m_index, np.random.randint(10, 350, idx))))
model.solar_profile = Param(model.m_index, initialize=dict(zip(model.m_index, np.random.uniform(0, 0.6, idx))))

# 参数：发电机容量
model.gen1_cap = Param(initialize=100)
model.gen2_cap = Param(initialize=100)
model.gen3_cap = Param(initialize=100)
model.backup_cap = Param(initialize=300)

# 参数：发电机运行成本
model.gen1_cost = Param(initialize=10)
model.gen2_cost = Param(initialize=10)
model.gen3_cost = Param(initialize=10)
model.backup_cost = Param(initialize=-3)

# 变量：发电机输出能量
model.gen1_use = Var(model.m_index, domain=NonNegativeReals)
model.gen2_use = Var(model.m_index, domain=NonNegativeReals)
model.gen3_use = Var(model.m_index, domain=NonNegativeReals)
model.backup_use = Var(model.m_index, domain=NonNegativeReals)

# 二进制变量：发电机3的状态
model.gen3_status = Var(model.m_index, domain=Binary)


# 目标函数：最大化总成本
def production_cost(model):
    total_cost = sum(
        model.gen1_use[m] * model.gen1_cost +
        model.gen2_use[m] * model.gen2_cost +
        model.gen3_use[m] * model.gen3_cost * model.gen3_status[m] + # gen3成本与状态挂钩
        model.backup_use[m] * model.backup_cost
        for m in model.m_index)
    return total_cost

model.obj = Objective(rule=production_cost, sense=maximize)

# Big-M 方法实现条件逻辑
# 一些常量
eps = 1e-3  # 用于创建间隙
bigm = 1e3  # 大M值

# 创建2个Big-M约束来控制gen3_status
# 条件：model.gen1_use[m] + model.gen2_use[m] >= 0.30 * model.load_profile[m]
# 结果：model.gen3_status[m] = 1 (如果条件满足) 或 0 (如果条件不满足)
def gen3_on_off_bigm_rule1(model, m):
    # 如果 gen3_status[m] = 1，则 L.H.S >= R.H.S + eps
    # 如果 gen3_status[m] = 0，则 L.H.S >= R.H.S + eps - bigm (几乎总是满足)
    return model.gen1_use[m] + model.gen2_use[m] >= 0.30 * model.load_profile[m] + eps - bigm * (1 - model.gen3_status[m])

def gen3_on_off_bigm_rule2(model, m):
    # 如果 gen3_status[m] = 0，则 L.H.S <= R.H.S
    # 如果 gen3_status[m] = 1，则 L.H.S <= R.H.S + bigm (几乎总是满足)
    return model.gen1_use[m] + model.gen2_use[m] <= 0.30 * model.load_profile[m] + bigm * model.gen3_status[m]

model.gen3_on_off_bigm1 = Constraint(model.m_index, rule=gen3_on_off_bigm_rule1)
model.gen3_on_off_bigm2 = Constraint(model.m_index, rule=gen3_on_off_bigm_rule2)


def energy_balance(model, m):
    eq = model.gen1_use[m] + model.gen2_use[m] + model.gen3_use[m] + model.backup_use[m] >= model.load_profile[m]
    return eq

model.energy_balance = Constraint(model.m_index, rule=energy_balance)

def gen1_max(model, m):
    eq = model.gen1_use[m] <= model.gen1_cap
    return eq
model.gen1_max = Constraint(model.m_index, rule=gen1_max)

def gen2_max(model, m):
    eq = model.gen2_use[m] <= model.gen2_cap
    return eq
model.gen2_max = Constraint(model.m_index, rule=gen2_max)

def gen3_max(model, m):
    eq = model.gen3_use[m] <= model.gen3_cap * model.solar_profile[m] * model.gen3_status[m] # gen3容量也与状态挂钩
    return eq
model.gen3_max = Constraint(model.m_index, rule=gen3_max)

def backup_max(model, m):
    eq = model.backup_use[m] <= model.backup_cap
    return eq
model.backup_max = Constraint(model.m_index, rule=backup_max)


# 求解器配置
Solver = SolverFactory('gurobi') # 确保Gurobi已安装并配置
Solver.options['LogFile'] = "gurobiLog"
print('\nConnecting to Gurobi Server...')
results = Solver.solve(model)

if (results.solver.status == SolverStatus.ok):
    if (results.solver.termination_condition == TerminationCondition.optimal):
        print("\n\n***Optimal solution found***")
        print('obj returned:', round(value(model.obj), 2))
    else:
        print("\n\n***No optimal solution found***")
        if (results.solver.termination_condition == TerminationCondition.infeasible):
            print("Infeasible solution")
            exit()
else:
    print("\n\n***Solver terminated abnormally***")
    exit()

# 结果提取与输出
load = []
gen1 = []
gen2 = []
gen3 = []
solar_profile = []
backup = []
gen3_status = []

for i in range(idx):
    load.append(value(model.load_profile[i]))
    gen1.append(value(model.gen1_use[i]))
    gen2.append(value(model.gen2_use[i]))
    gen3.append(value(model.gen3_use[i]))
    solar_profile.append(value(model.solar_profile[i]))
    backup.append(value(model.backup_use[i]))
    gen3_status.append(value(model.gen3_status[i]))

df = pd.DataFrame({
    'Load profile': load,
    'Gen 1': gen1,
    'Gen 2': gen2,
    'Gen 3': gen3,
    'Solar profile': solar_profile,
    'Backup': backup,
    'Gen 3 Status': gen3_status,
})
df.to_excel('binary.xlsx')
print(df)