添加min_downtime约束后Pyomo CCGT调度优化模型失效求排查

背景 #

我正在编写基于电价完美预见的Pyomo脚本，用于CCGT机组的最优调度。X1、X2、X5为二进制变量，分别对应CCGT的全基载模式、部分负载模式和停机模式。

问题 #

未加入min_downtime约束时，脚本运行完全符合预期。

请求协助 #

请帮忙检查我的代码中min_downtime约束存在的问题，需求为：当X5取值为1时，CCGT需至少保持停机状态6小时。

代码 #

from pyomo.environ import *
from pyomo.opt import SolverFactory
import pandas as pd
import matplotlib.pyplot as plt

model = ConcreteModel()

inputs_df = pd.read_excel('input.xlsx')

electricity_prices = pd.DataFrame(inputs_df['PTF'])
ancillary_service_prices = pd.DataFrame(inputs_df['SFC'])
gas_prices = pd.DataFrame(inputs_df['Gasp'])
msp_prices = pd.DataFrame(inputs_df['AUF'])
model.max_gas_quantity = Param(initialize=150000000)  # or any other value

hours = list(electricity_prices.index)

model.X1 = Var(hours, within=Binary)
model.X2 = Var(hours, within=Binary)
model.X5 = Var(hours, within=Binary)
model.startup = Var(hours, within=NonNegativeIntegers)
model.shutdown = Var(hours, within=NonNegativeIntegers)
model.plant_running = Var(hours, within=Binary)
model.plant_idling = Var(hours, within=Binary)


# generation
def generation(model, hour):
    return model.X1[hour]*770 + model.X2[hour]*550 + model.X5[hour]*0 + model.startup[hour]*225 + model.shutdown[hour]*70
model.generation = Expression(hours, rule=generation)

# ancillary services
def ancillary_service(model, hour):
    return model.X1[hour]*0 + model.X2[hour]*220 + model.X5[hour]*0
model.ancillary_service = Expression(hours, rule=ancillary_service)

# gas consumption
def gas_consumption(model, hour):
    return model.X1[hour]*143990 + model.X2[hour]*108350 + model.X5[hour]*0 + model.startup[hour]*58500 + model.shutdown[hour]*17500
model.gas_consumption = Expression(hours, rule=gas_consumption)

# EOH consumption
def eoh_consumption(model, hour):
    return model.X1[hour]*2 + model.X2[hour]*2 + model.X5[hour]*0 + model.startup[hour]*20
model.eoh_consumption = Expression(hours, rule=eoh_consumption)

# Revenues
def electricity_revenue(model, hour):
    return model.generation[hour] * electricity_prices.loc[hour, "PTF"]
model.electricity_revenue = Expression(hours, rule=electricity_revenue)

def ancillary_service_revenue(model, hour):
    return model.ancillary_service[hour] * ancillary_service_prices.loc[hour, "SFC"]
model.ancillary_service_revenue = Expression(hours, rule=ancillary_service_revenue)

# Costs
def gas_cost(model, hour):
    return model.gas_consumption[hour] * gas_prices.loc[hour, "Gasp"]
model.gas_cost = Expression(hours, rule=gas_cost)

def eoh_cost(model, hour):
    return model.eoh_consumption[hour] * 7500
model.eoh_cost = Expression(hours, rule=eoh_cost)

def tso_var_cost(model, hour):
    return model.generation[hour] * 93.25
model.tso_var_cost = Expression(hours, rule=tso_var_cost)

def msp_cost(model, hour):
    if electricity_prices.loc[hour, "PTF"] > msp_prices.loc[hour, "AUF"] :
        return model.generation[hour] * (electricity_prices.loc[hour, "PTF"] - msp_prices.loc[hour, "AUF"])
    else:
        return 0
model.msp_cost = Expression(hours, rule=msp_cost)

# Constraints
def one_mode_per_hour(model, hour):
    return model.X1[hour] + model.X2[hour] + model.X5[hour] == 1
model.one_mode_per_hour = Constraint(hours, rule=one_mode_per_hour)

def running(model, hour):
    return model.plant_running[hour] >= 1 - model.X5[hour]
model.running = Constraint(hours, rule=running)

def idling(model, hour):
    return model.plant_idling[hour] >= model.X5[hour]
model.idling = Constraint(hours, rule=idling)

def start_flag(model, hour):
    if hour > len(hours)-2:
        return Constraint.Skip
    else:
        return model.startup[hour] >= model.X5[hour] - model.X5[hour+1]
model.start_flag = Constraint(hours, rule=start_flag)

def stop_flag(model, hour):
    if hour > len(hours)-2:
        return Constraint.Skip
    else:
        return model.shutdown[hour+1] >= model.X5[hour+1] - model.X5[hour]
model.stop_flag = Constraint(hours, rule=stop_flag)


window_size = 6
def rolling_min_downtime(model, hour):
    preceding_periods = {hour_prime for hour_prime in hours if hour - window_size <= hour_prime < hour}
    return 6 - sum(model.X5[hour_prime] for hour_prime in preceding_periods) <= 1
eval_periods = {hour for hour in hours if hour >= window_size}
model.rolling_min_downtime = Constraint(eval_periods, rule=rolling_min_downtime)


def total_gas(model):
    return sum(model.gas_consumption[hour] for hour in hours) <= model.max_gas_quantity
model.total_gas = Constraint(rule=total_gas)

# Objective Function
def objective(model):
    return sum(model.electricity_revenue[hour] for hour in hours) + sum(model.ancillary_service_revenue[hour] for hour in hours) \
- sum(model.gas_cost[hour] for hour in hours) - sum(model.eoh_cost[hour] for hour in hours) - sum(model.tso_var_cost[hour] for hour in hours) - sum(model.msp_cost[hour] for hour in hours)

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


solver = SolverFactory('glpk')
result = solver.solve(model)

问题分析与修正 #

当前的rolling_min_downtime约束逻辑完全错误，没有正确实现“停机后至少保持6小时停机状态”的需求。当前约束6 - sum(X5) <=1等价于sum(X5) >=5，意味着过去6小时里至少有5小时是停机状态，和需求完全不符。

正确的最小停机时间约束应该从停机触发后的连续时段入手：当机组从运行状态转为停机状态（即X5[hour] =1且X5[hour-1]=0）时，接下来的5个小时（加上当前小时共6小时）必须保持X5=1。

修正后的约束实现 #

替换原有的rolling_min_downtime约束，改用以下代码：

window_size = 6

# 最小停机时间约束：停机后必须连续保持6小时停机状态
def min_downtime_constraint(model, hour):
    # 从第5小时开始生效（对应window_size=6）
    if hour < window_size -1:
        return Constraint.Skip
    # 若在hour - window_size +1时刻触发停机（从运行转停机），则当前hour必须保持停机
    return model.X5[hour] >= model.X5[hour - window_size +1] - model.X5[hour - window_size]
model.min_downtime = Constraint(hours[window_size-1:], rule=min_downtime_constraint)

# 处理初始状态：如果第0小时是停机，强制前6小时全部保持停机
def initial_downtime_constraint(model):
    if hours[0] == 0 and len(hours) >= window_size:
        return sum(1 - model.X5[h] for h in hours[:window_size]) <= 0
    return Constraint.Skip
model.initial_downtime = Constraint(rule=initial_downtime_constraint)

约束逻辑说明 #

1. 主约束：对于每个小时hour（从第5小时开始），检查是否在hour-5时刻发生了停机触发（即X5[hour-5] - X5[hour-6] =1）。如果是，则强制X5[hour] =1，确保从触发停机后连续6小时都是停机状态。
2. 初始状态约束：如果调度的第一个小时机组已经处于停机状态，强制要求前6小时全部保持停机，避免初始停机状态不满足最小时长要求。

其他优化建议 #

• startup和shutdown变量定义为NonNegativeIntegers不合理，启停事件是二进制的（要么0要么1），应改为within=Binary，既符合实际逻辑，也能减少求解器计算量。
• plant_running和plant_idling变量可简化：因为X1+X2+X5=1，所以plant_running[hour] = X1[hour]+X2[hour]，plant_idling[hour] = X5[hour]，可直接用表达式替代，避免冗余变量和不等式约束。

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