Background :
I am trying to write an pyomo script to optimally dispatch a ccgt plant based on perfect foresight of electricity prices. x1, x2 and x5 are binary variables and each of them are corresponding to different modes of ccgt where x1 is full base load mode, x2 is part load mode and x5 is idle mode.
Problem :
The script functions properly and precisely as intended in the absence of the min_downtime limitation.
Help requested :
Could someone please have a look at my code and see what I am missing in the min_downtime constraint? The ccgt needs to be idle for minimum 6 hours if X5 selected 1.
Code :
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)
Take a look at below. I modified a handful of things. Can't test as no data is provided. If this doesn't work as intended and you are stuck, edit your question and include about 20 rows of data in a table.
The logic implemented is:
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)