PYOMO目标函数中结合Param与变量值的方法及报错求助

一、直接解决load_profile[0]未定义的报错 #

报错核心是load_profile这个Param的索引0没有赋值，按以下步骤排查修复：

• 检查Param初始化代码：确保所有时间索引（包括0）都有对应数值。比如时间集是RangeSet(0,23)，初始化时要覆盖所有时间点：

  model.load_profile = pyo.Param(model.Time, initialize={0: 150, 1: 160, ...})  # 手动赋值所有时间点
  

• 外部数据加载校验：如果从CSV/Excel加载需求曲线，确认文件包含索引0的行，且加载逻辑的索引映射正确。
• 临时兜底方案：给Param设置默认值（不推荐长期使用，建议确保所有索引都有有效值）：

  model.load_profile = pyo.Param(model.Time, default=0)
  

二、目标函数中Param与变量结合的正确方式 #

Pyomo线性求解器不支持直接使用Python原生min()/max()，需通过辅助变量+约束线性化非线性表达式，具体实现如下：

1. 定义核心变量 #

# 光伏、风电发电量变量
model.solar_gen = pyo.Var(model.Time, within=pyo.NonNegativeReals)
model.wind_gen = pyo.Var(model.Time, within=pyo.NonNegativeReals)

2. 拆分收益组件并线性化 #

（1）实际收益：Min(总发电量, 需求曲线) * PPA

引入辅助变量actual_delivery，通过约束限制其取值范围：

model.actual_delivery = pyo.Var(model.Time, within=pyo.NonNegativeReals)
# 约束1：实际供电量不超过总发电量
model.constr_actual1 = pyo.Constraint(model.Time, rule=lambda m,t: m.actual_delivery[t] <= m.solar_gen[t]+m.wind_gen[t])
# 约束2：实际供电量不超过需求曲线
model.constr_actual2 = pyo.Constraint(model.Time, rule=lambda m,t: m.actual_delivery[t] <= m.load_profile[t])

实际收益表达式：sum(model.actual_delivery[t] * PPA for t in model.Time)

（2）超额收益：(Max(总发电量, 需求曲线)-需求曲线) * 0.5*PPA

引入辅助变量excess_gen，约束其为非负的超额发电量：

model.excess_gen = pyo.Var(model.Time, within=pyo.NonNegativeReals)
# 约束：超额发电量≥总发电量-需求曲线（小于0时取0）
model.constr_excess = pyo.Constraint(model.Time, rule=lambda m,t: m.excess_gen[t] >= m.solar_gen[t]+m.wind_gen[t]-m.load_profile[t])

超额收益表达式：sum(model.excess_gen[t] * 0.5*PPA for t in model.Time)

（3）短缺惩罚：当总发电量<90%需求曲线时，惩罚(总发电量-0.9*需求曲线)*惩罚率

引入辅助变量shortfall，约束其为非负的短缺量：

model.shortfall = pyo.Var(model.Time, within=pyo.NonNegativeReals)
# 约束：短缺量≥0.9*需求曲线-总发电量（大于0时取0）
model.constr_shortfall = pyo.Constraint(model.Time, rule=lambda m,t: m.shortfall[t] >= 0.9*m.load_profile[t]-(m.solar_gen[t]+m.wind_gen[t]))

短缺惩罚表达式：sum(model.shortfall[t] * penalty_rate for t in model.Time)

3. 构建目标函数 #

def obj_rule(m):
    actual_rev = sum(m.actual_delivery[t] * PPA for t in m.Time)
    excess_rev = sum(m.excess_gen[t] * 0.5*PPA for t in m.Time)
    penalty_cost = sum(m.shortfall[t] * penalty_rate for t in m.Time)
    return actual_rev + excess_rev - penalty_cost

model.obj = pyo.Objective(rule=obj_rule, sense=pyo.maximize)

三、完整可运行代码片段 #

import pyomo.environ as pyo

# 全局参数定义
PPA = 60
DFR = 0.9
penalty_rate = 120
time_range = range(24)

# 创建模型
model = pyo.ConcreteModel()
model.Time = pyo.RangeSet(0,23)

# 初始化需求曲线（示例数据）
load_data = {t: 120 + t*3 for t in time_range}
model.load_profile = pyo.Param(model.Time, initialize=load_data)

# 变量定义
model.solar_gen = pyo.Var(model.Time, within=pyo.NonNegativeReals)
model.wind_gen = pyo.Var(model.Time, within=pyo.NonNegativeReals)
model.actual_delivery = pyo.Var(model.Time, within=pyo.NonNegativeReals)
model.excess_gen = pyo.Var(model.Time, within=pyo.NonNegativeReals)
model.shortfall = pyo.Var(model.Time, within=pyo.NonNegativeReals)

# 约束定义
model.constr_actual1 = pyo.Constraint(model.Time, rule=lambda m,t: m.actual_delivery[t] <= m.solar_gen[t]+m.wind_gen[t])
model.constr_actual2 = pyo.Constraint(model.Time, rule=lambda m,t: m.actual_delivery[t] <= m.load_profile[t])
model.constr_excess = pyo.Constraint(model.Time, rule=lambda m,t: m.excess_gen[t] >= m.solar_gen[t]+m.wind_gen[t]-m.load_profile[t])
model.constr_shortfall = pyo.Constraint(model.Time, rule=lambda m,t: m.shortfall[t] >= 0.9*m.load_profile[t]-(m.solar_gen[t]+m.wind_gen[t]))

# 目标函数
model.obj = pyo.Objective(rule=lambda m: sum(m.actual_delivery[t]*PPA + m.excess_gen[t]*0.5*PPA - m.shortfall[t]*penalty_rate for t in m.Time), sense=pyo.maximize)

# 求解
solver = pyo.SolverFactory('glpk')
result = solver.solve(model)
print(result)

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