基于PYOMO的收益最大化容量扩展模型构建及约束问题排查

问题描述 #

我用Pyomo构建了一个容量扩展模型，包含两台主发电机gen1、gen2和一台备用机组lost_load。核心逻辑为：gen1、gen2运行满足负荷曲线获取收益，同时它们的容量存在负向固定成本。

原模型代码如下：

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


model = ConcreteModel()
np.random.seed(24)
load_profile = np.random.randint(90, 120, 24)
model.m_index = Set(initialize=list(range(len(load_profile))))

model.grid = Var(model.m_index, domain=NonNegativeReals)

# Gen1变量
model.gen1_cap = Var(domain=NonNegativeReals)
model.gen1_use = Var(model.m_index, domain=NonNegativeReals)

# Gen2变量
model.gen2_cap = Var(domain=NonNegativeReals)
model.gen2_use = Var(model.m_index, domain=NonNegativeReals)

# 负荷曲线参数
model.load_profile = Param(model.m_index, initialize=dict(zip(model.m_index, load_profile)))

model.lost_load = Var(model.m_index, domain=NonNegativeReals)

# 目标函数
def revenue(model):
    total_revenue = sum(
        model.gen1_use[m] * 5.2 +
        model.gen2_use[m] * 6.1 +
        model.lost_load[m] * -100
        for m in model.m_index)

    total_fixed_cost = model.gen1_cap * -45 + model.gen2_cap * -50

    total_cost = total_revenue + total_fixed_cost
    return total_cost


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

# 能量平衡约束（原设置为<=时出错）
def energy_balance1(model, m):
    return model.grid[m] <= model.gen1_use[m] + model.gen2_use[m] + model.lost_load[m]


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


def grid_limit(model, m):
    return model.grid[m] == model.load_profile[m]


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


def max_gen1(model, m):
    eq = model.gen1_use[m] <= model.gen1_cap
    return eq


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


def max_gen2(model, m):
    eq = model.gen2_use[m] <= model.gen2_cap
    return eq


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


Solver = SolverFactory('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()

grid_use = []
gen1 = []
gen2 = []
lost_load = []
load = []

for i in range(len(load_profile)):
    grid_use.append(value(model.grid[i]))
    gen1.append(value(model.gen1_use[i]))
    gen2.append(value(model.gen2_use[i]))
    lost_load.append(value(model.lost_load[i]))
    load.append(value(model.load_profile[i]))

print('gen1 capacity: ', value(model.gen1_cap))
print('gen2 capacity: ', value(model.gen2_cap))

pd.DataFrame({
    'Grid': grid_use,
    'Gen1': gen1,
    'Gen2': gen2,
    'Shortfall': lost_load,
    'Load': load
}).to_excel('capacity expansion.xlsx')

当能量平衡约束设为等式时，模型正常运行，得到的最优容量为：

gen1 capacity:  9MW
gen2 capacity:  108MW

但我需要将能量平衡约束设为<=，目的是让gen1和gen2全天24时段100%满负荷运行——毕竟目标函数中gen1_use和gen2_use对应正向收益，理论上满发能最大化收益。但修改后模型报不可行错误，且当前结果不符合我期望的调度效果（期望两台机组全程满发，容量匹配负荷峰值）。

问题原因 #

1. 约束逻辑与问题性质冲突：将能量平衡约束改为grid[m] <= gen1_use[m]+gen2_use[m]+lost_load[m]后，模型无多余电量消纳途径，且目标函数中gen1_cap、gen2_cap的总收益系数为正（gen1每MW年收益：5.2*24-45=79.8，gen2：6.1*24-50=96.4），导致模型可无限增大容量提升收益，属于无界问题，而非不可行，你看到的提示可能是对求解器输出的误判。

2. 未明确满发约束：仅靠目标函数的正向收益无法强制机组满发，模型会自动权衡容量固定成本与运行收益，选择最优组合而非强制满发。

解决方案 #

要实现机组全天满发的需求，需添加强制满发约束，同时调整能量平衡约束逻辑：

1. 添加强制满发约束 #

在模型中加入以下代码，确保gen1和gen2在所有时段都运行在最大容量：

# 强制gen1全天满负荷运行
def gen1_full_load(model, m):
    return model.gen1_use[m] == model.gen1_cap
model.gen1_full_load = Constraint(model.m_index, rule=gen1_full_load)

# 强制gen2全天满负荷运行
def gen2_full_load(model, m):
    return model.gen2_use[m] == model.gen2_cap
model.gen2_full_load = Constraint(model.m_index, rule=gen2_full_load)

2. 调整能量平衡约束 #

保持能量平衡约束为等式（多余电量无消纳途径，发电量+失负荷必须等于负荷）：

def energy_balance1(model, m):
    return model.grid[m] == model.gen1_use[m] + model.gen2_use[m] + model.lost_load[m]
model.energy_balance1 = Constraint(model.m_index, rule=energy_balance1)

3. 可选：限制总容量匹配负荷峰值 #

若希望机组总容量刚好覆盖最大负荷（避免不必要的容量成本），可添加以下约束：

max_load = max(load_profile)
model.total_capacity_constraint = Constraint(expr=model.gen1_cap + model.gen2_cap == max_load)

修改后核心约束代码 #

# 强制满发约束
def gen1_full_load(model, m):
    return model.gen1_use[m] == model.gen1_cap
model.gen1_full_load = Constraint(model.m_index, rule=gen1_full_load)

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

# 能量平衡约束设为等式
def energy_balance1(model, m):
    return model.grid[m] == model.gen1_use[m] + model.gen2_use[m] + model.lost_load[m]
model.energy_balance1 = Constraint(model.m_index, rule=energy_balance1)

# 可选：总容量匹配负荷峰值
max_load = max(load_profile)
model.total_capacity_constraint = Constraint(expr=model.gen1_cap + model.gen2_cap == max_load)

效果说明 #

修改后，模型会强制gen1和gen2全天满发，同时由于lost_load的高惩罚，会优先让总容量覆盖最大负荷以避免失负荷。若添加了总容量约束，模型会根据单位容量收益（gen2更高）分配两台机组的容量，最大化整体收益。

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