背景：
我正在尝试编写一个 pyomo 脚本，以根据对电价的完美预见来优化调度 ccgt 电厂。 x1、x2、x5是二元变量，分别对应ccgt的不同模式，其中x1是全基本负载模式，x2是部分负载模式，x5是空闲模式。
问题：
在没有 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)
请看下面。我修改了一些东西。无法测试，因为没有提供数据。如果这不能按预期工作并且您被卡住了，请编辑您的问题并在表中包含大约 20 行数据。
实现的逻辑是：
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)
# NOT NEEDED. If you know X5, then you know IDLE vs. RUNNING so this is duplicate
# Information...
# 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)
# NOT NEEDED ?
# 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)
# re-aligned to capture the first period where mode == 5
def stop_flag(model, hour):
return model.shutdown[hour] >= model.X5[hour] - model.X5[hour-1]
model.stop_flag = Constraint(hours[1:], rule=stop_flag)
# now that we have the start of shutdown flagged, we just need to suppress
# modes 1 and 2 for the next 5 time periods
suppress_periods = 5
def suppress_running_modes(model, hour):
following_periods = {hour_prime for hour_prime in hours
if hour_prime > hour
and hour_prime <= hour + suppress_periods}
return sum(model.X1[h] + model.X2[h] for h in following_periods) <= suppress_periods*(1-model.shutdown[hour]))
model.suppress_start = Constraint(hours, rule=suppress_running_modes)
# need to capture a startup event. Note: this doesnt "charge" for startup in 1st period... model could
# be running at start with no charge....
def startup_capture(model, hour):
return model.startup[hour] >= model.X5[hour - 1] - model.X5[hour]
model.startup_capture = Constraint(hours[1:], rule=startup_capture)
# 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)