Pyomo optimisation not working (gas plant dispatch)

background

I am trying to write an pyomo script to optimally dispatch a gas plant based on perfect foresight of electricity prices. I believe I am 90% of the way there, just a few issues.

Problem

My script works, but the solver is never dispatching the plant, even where it should be, in the example provided below, manually I can calculate at least $8131 of potential profit.

I suspect the reason for my zero results is due to how I've written the constraints, of which there are 2;

Gas Plant takes 10 minutes to boot up from a cold start
Once warmed up, the gas plant has a min load it must operate at/above.
Gas Plant can only consume 9000 GJ of gas in a single day

Specifically on further testing, I think it is the 'gas_volume_used' constraint which is causing the issue.

Help Requested

Could someone please have a look at my code and see what I am missing in the constraint equations?

Code

import pandas as pd
from pyomo.environ import *
from pyomo.core import *


# ---------------------
# Main Functions
# ---------------------
def model_to_df(model, first_period, last_period):

    # Need to increase the first & last period by 1 because of pyomo indexing
    # and add 1 to the value of last model period because of the range
    periods = range(model.T[first_period + 1], model.T[last_period + 1] + 1)
    rrp = [model.P.extract_values()[None][i] for i in periods]
    Gascold = [value(model.Gascold[i]) for i in periods]
    Gaswarm = [value(model.Gaswarm[i]) for i in periods]
    Gashot = [value(model.Gashot[i]) for i in periods]
    GasUnit_mw = [value(model.GasUnit_mw[i]) for i in periods]
    GasUse = [value(model.GasUse[i]) for i in periods]

    df_dict = dict(
        period=periods,
        rrp=rrp,
        Gascold=Gascold,
        Gaswarm=Gaswarm,
        Gashot=Gashot,
        GasUnit_mw=GasUnit_mw,
        GasUse=GasUse
    )

    df = pd.DataFrame(df_dict)

    return df


def optimize_year(df, gunits, gnet, gprice, vom, cold_gas_p_unit, hr, max_gas, sys_use,
                  first_model_period, last_model_period):

    """
    Optimize behavior of a gas plant over a time period
    Assume perfect foresight of electricity prices.

    Parameters
    ----------
    df : dataframe
        dataframe with columns of hourly LBMP and the hour of the year
    gunits : int
        number of gas turbines
    gnet : double
        size in MW's of each gas turbine
    gprice : double
        the price of gas $/GJ
    vom  : double
        the variable cost of gas $/MWh dispatched
    cold_gas_p_unit  : double
        the amount of gas in GJ used in 5 minutes while booting up from cold start
    hr  : double
        heat rate, the amount of gas consumed per MWh produced (GJ/MWh)
    sys_use  : double
        auxillary gas usage, percentage
    max_gas  : double
        the maximum amount of gas (GJ) available to use in a period
    first_model_period : int, optional
        Set the first period of the year to be considered in the optimization
    last_model_period : int, optional
        Set the last period of the year to be considered in the optimization

    Returns
    -------
    dataframe
        Gas plant operational stats and time stamp
    """

    # Filter the data
    df = df.loc[first_model_period:last_model_period, :]

    model = ConcreteModel()

    # Define model parameters
    model.T = Set(initialize=df.period.tolist(), doc='periods of year', ordered=True)
    model.P = Param(initialize=df.rrp.tolist(), doc='price for each period', within=Any)

    model.MaxGas = Param(initialize=(gnet/12) * gunits * hr * (1+sys_use), doc='Max gas usage (GJ) in interval', within=Any)
    model.MaxVol = Param(initialize=max_gas, doc='Max gas usage (GJ) in 24 hour period', within=Any)
    model.MaxMW = Param(initialize=(gnet/12) * gunits, doc='Max MWs', within=Any)

    model.Gascold = Var(model.T, doc='binary', within=Binary, initialize=1)
    model.Gaswarm = Var(model.T, doc='binary', within=Binary, initialize=0)
    model.Gashot = Var(model.T, doc='binary', within=Binary, initialize=0)
    model.GasUse = Var(model.T, bounds=(0, model.MaxGas))
    model.GasUnit_mw = Var(model.T, bounds=(0, model.MaxMW))

    # *********************
    # GAS CONSTRAINTS
    # *********************
    "Ensure that all three Gas State flags sum to 1"
    def set_on(model, t):
        return model.Gascold[t] + model.Gaswarm[t] + model.Gashot[t] == 1

    model.limit_flags = Constraint(model.T, rule=set_on)

    "Set t0 to GAS COLD = 1"
    def gas_cold(model, t):
        if t == model.T.first():
            return model.Gascold[t] == 1
        else:
            return model.Gascold[t] >= 0

    model.gas_cold_check = Constraint(model.T, rule=gas_cold)

    "Gas can only warm up from cold"
    def gas_warm(model, t):
        if t < 2:
            return model.Gaswarm[t] == 0
        elif (model.Gaswarm[t-2] + model.Gaswarm[t-1] + model.Gascold[t-1]) == 1:
            return model.Gaswarm[t] >= 0
        else:
            return model.Gaswarm[t] == 0

    model.gas_warm_check = Constraint(model.T, rule=gas_warm)

    "gas can only get hot from warm"
    def gas_hot(model, t):
        if t <= 2:
            return model.Gashot[t] == 0
        elif value(model.Gaswarm[t-2] + model.Gaswarm[t-1] + model.Gashot[t-2] + model.Gashot[t-1]) < 2:
            return model.Gashot[t] == 0
        else:
            return model.Gashot[t] >= 0

    model.gas_hot_check = Constraint(model.T, rule=gas_hot)

    "Gas hot min load"
    def gas_hot_min_load(model, t):
        if model.Gashot[t] == 1:
            return model.GasUnit_mw[t] >= 30 / 12
        else:
            return model.GasUnit_mw[t] == 0

    model.gas_hot_min_load_check = Constraint(model.T, rule=gas_hot_min_load)

    "Gas volume used"
    def gas_volume_used(model, t):
        'how much gas is being used?'
        return model.GasUse[t] == (model.GasUnit_mw[t] * hr * (1+sys_use)) + (model.Gaswarm[t] * gunits * cold_gas_p_unit * (1+sys_use))

    model.gas_vol = Constraint(model.T, rule=gas_volume_used)

    "only 9000 GJ of gas can be used in a single day"
    def daily_max_gas(model, t):
        min_t = model.T.first()
        max_t = model.T.last()

        # Check all t until the last 24 hours
        if t <= max_t:
            return sum(model.GasUse[i] for i in range(min_t, max_t+1)) <= model.MaxVol
        else:
            return Constraint.Skip

    model.max_gas_vol = Constraint(model.T, rule=daily_max_gas)

    # Define the income, expenses, and profit
    gas_income = sum(df.loc[t, 'rrp'] * model.GasUnit_mw[t] for t in model.T)
    gas_expenses = sum(model.GasUse[t] * gprice + model.GasUnit_mw[t] * vom for t in model.T)
    profit = gas_income - gas_expenses
    model.objective = Objective(expr=profit, sense=maximize)

    # Solve the model
    solver = SolverFactory('glpk')
    solver.solve(model)

    results_df = model_to_df(model, first_period=first_model_period, last_period=last_model_period)
    results_df['local_time'] = df.loc[:, 'local_time']

    return results_df


# *********************
# RUN SCRIPT
# *********************
# Import file
df = pd.read_csv('raw_prices.csv')

# New Index
df['period'] = df.reset_index().index

# get current Year
first_model_period = 0
last_model_period = df['period'].iloc[-1]

# Gas Inputs
gunits = 3 # No units
gnet = 58.5 # MW/unit
srmc = 115.49 # $/MWh
vom = 6.18 # $/MWh
gprice = 10.206 # $/GJ
cold_gas_p_unit = 15.75 # GJ/unit/5mins
hr = 10.5 # GJ/MWh
sys_use = 0.02 # %
max_gas = 9000 # GJ

# Run Model
results_df = optimize_year(df=df, gunits=gunits, gnet=gnet, gprice=gprice, vom=vom, cold_gas_p_unit=cold_gas_p_unit,
                           hr=hr, max_gas=max_gas, sys_use=sys_use,
                           first_model_period=first_model_period, last_model_period=last_model_period)

# ---------------------
# Write the file out and remove the index
# ---------------------
print('ran successfully')
results_df.to_csv('optimised_output.csv', index=False)

Example Data

local_time      rrp
3/07/2022 16:40 64.99
3/07/2022 16:45 73.10744
3/07/2022 16:50 74.21746
3/07/2022 16:55 95.90964
3/07/2022 17:00 89.60728
3/07/2022 17:05 90.37018
3/07/2022 17:10 92.52947
3/07/2022 17:15 100.1449
3/07/2022 17:20 103.70199
3/07/2022 17:25 290.2848
3/07/2022 17:30 115.40232
3/07/2022 17:35 121.00009
3/07/2022 17:40 294.62543
3/07/2022 17:45 121
3/07/2022 17:50 294.69446
3/07/2022 17:55 124.90664
3/07/2022 18:00 117.21929
3/07/2022 18:05 115.60193
3/07/2022 18:10 121
3/07/2022 18:15 148.20186
3/07/2022 18:20 123.11763
3/07/2022 18:25 120.99996
3/07/2022 18:30 121
3/07/2022 18:35 121
3/07/2022 18:40 121
3/07/2022 18:45 120.77726
3/07/2022 18:50 105.44435
3/07/2022 18:55 104.22928
3/07/2022 19:00 102.26646
3/07/2022 19:05 99.5896
3/07/2022 19:10 95.5881
3/07/2022 19:15 88
3/07/2022 19:20 88
3/07/2022 19:25 88
3/07/2022 19:30 96.97537
3/07/2022 19:35 88
3/07/2022 19:40 86.39322
3/07/2022 19:45 88
3/07/2022 19:50 85.97703
3/07/2022 19:55 86
3/07/2022 20:00 87.61047
3/07/2022 20:05 88
3/07/2022 20:10 88
3/07/2022 20:15 85.34908
3/07/2022 20:20 85.1
3/07/2022 20:25 71.01181
3/07/2022 20:30 84.88422
3/07/2022 20:35 85.39669
3/07/2022 20:40 86.02542
3/07/2022 20:45 84.87425
3/07/2022 20:50 70.0652
3/07/2022 20:55 70.65106
3/07/2022 21:00 70.53063
3/07/2022 21:05 68.27
3/07/2022 21:10 68.86159
3/07/2022 21:15 68.32018
3/07/2022 21:20 67.30473
3/07/2022 21:25 66.78292
3/07/2022 21:30 66.02195
3/07/2022 21:35 66.40755

Solution

Well, I went a little geek on this. Got hooked, kinda interesting problem.

So, I made a bunch of changes and left some of your code in this example. I also chopped down a handful of the cost variables and made them rather simple as I was getting a little lost in the sauce and so that I was (mostly) convinced things were working, so the units/conversions/costs are a bit nonsensical now, but should be easily recovered.

Hopefully there are a couple concepts in here that you can use as you work through this. A few notes...

Needed a binary variable to indicate that the plant was started, and another to keep track of whether it was "running" or not in a particular period, these were linked with a constraint
Added a little trickery with the time windows to make a rolling evaluation period for total gas use
Added a minimum use for the plant to run or else once it was "started" it could arbitrarily run with 0 gas when not profitable, now a minimum-run or off decision is forced

Plot shows pretty convincing evidence that it is running as hoped, it starts up, runs at max blast when price is high, and adheres to rolling gas limit, then shuts down and does it again.

Code

import pandas as pd
import numpy as np
from pyomo.environ import *
#from pyomo.core import *   # not needed
import matplotlib.pyplot as plt

# Import file
df = pd.read_csv('raw_prices_full.csv')

# New Index
df['period'] = df.reset_index().index

print(df.head())
print(df.info())

# get current Year
first_model_period = 0
last_model_period = df['period'].iloc[-1]

# Gas Inputs
gunits = 2 # No units
gnet = 20 # MW/unit
#srmc = 115.49 # $/MWh
#vom = 6.18 # $/MWh
gprice = 10 # $/GJ
startup_gas_rate = 8 # GJ/unit/5mins
min_running = 5 # GJ/unit/5mins   <--- assumed minimal output...  or plant could "run" at zero
hr = 10 # GJ/MWh
sys_use = 0.10 # % overhead in GJ of gas
max_gas = 100 # GJ/30 minutes (6 periods), cut down for testing

# Filter the data
df = df.loc[first_model_period:last_model_period, :]

model = ConcreteModel()

# Define model parameters
model.T = Set(initialize=df.period.tolist(), doc='periods of year', ordered=True)
# P needs to be indexed by T, and should be init by dictionary
model.P = Param(model.T, initialize=df.rrp.to_dict(), doc='price for each period', within=NonNegativeReals)  
model.MaxGas = Param(initialize=max_gas, doc='Max gas usage (GJ) in 24 hour period', within=NonNegativeReals)

# variables
model.turn_on = Var(model.T, domain=Binary)         # decision to turn the plant on in period t, making it available at t+2
model.plant_running = Var(model.T, domain=Binary)   # plant is running (able to produce heat) in period t.  init is not necessary
model.total_gas = Var(model.T, bounds=(0, gunits*gnet / 12 * hr * (1+sys_use)))  # the total gas flow in any period [GJ/5min]
model.gas_for_power = Var(model.T, domain=NonNegativeReals )  # the gas flow for power conversion at any period [GJ/5min]
model.profit = Var(model.T)  # non-essential variable, just used for bookeeping

# *********************
# GAS CONSTRAINTS
# *********************

# logic for the plant_running control variable
def running(model, t):
    # plant cannot be running in first 2 periods
    if t < 2:
        return model.plant_running[t] == 0
    else:
        # plant can be running if it was running in last period or started 2 periods ago.
        return model.plant_running[t] <= model.plant_running[t-1] + model.turn_on[t-2]
model.running_constraint = Constraint(model.T, rule=running)

# charge the warmup gas after a start decision.  Note the conditions in this constraint are all
# mutually exclusive, so there shouldn't be any double counting
# this will constrain the minimum gas flow to either the startup flow for 2 periods, 
# the minimum running flow, or zero if neither condition is met
def gas_mins(model, t):
    # figure out the time periods to inspect to see if a start event occurred...
    if t==0:
        possible_starts = model.turn_on[t]
    else:
        possible_starts = model.turn_on[t] + model.turn_on[t-1]
    return model.total_gas[t] >= \
            startup_gas_rate * gunits * (1 + sys_use) * possible_starts + \
            min_running * gunits * (1 + sys_use) * model.plant_running[t]
model.gas_mins = Constraint(model.T, rule=gas_mins)

# constrain the gas used for power to zero if not running, or the period max
def max_flow(model, t):
    return model.gas_for_power[t] <= model.plant_running[t] * gunits * gnet / 12 * hr
model.running_max_gas = Constraint(model.T, rule=max_flow)

# add the overhead to running gas to capture total gas by when running
def tot_gas_for_conversion(model, t):
    return model.total_gas[t] >=  model.gas_for_power[t] * (1 + sys_use)
model.tot_gas_for_conversion = Constraint(model.T, rule=tot_gas_for_conversion)


# I don't believe the constraint below is needed.  The model will try to maximize
# gas use when profit is available, and we have already established a minimum

# def gas_volume_used(model, t):
#     'how much gas is being used?'
#     'during cold start, a small amount of gas is used'
#     'during operation, gas is consumed based on the Heat Rate'
#     if model.plant_running[t] == 1 and model.GasUnit_mw[t] == 0:
#         return model.gas_use[t] == startup_gas_rate * gunits * (1+sys_use)
#     else:
#         return model.gas_use[t] == model.GasUnit_mw[t] * hr * (1+sys_use)

# model.gas_vol = Constraint(model.T, rule=gas_volume_used)

# this isn't doing what you think, you are not controlling the range of summation to 24 time periods
# def daily_max_gas(model, t):
#     "only 9000 GJ of gas can be used in a single day"
#     min_t = model.T.first()
#     max_t = model.T.last()

#     # Check all t until the last 24 hours
#     if t <= max_t:
#         return sum(model.gas_use[i] for i in range(min_t, max_t+1)) <= model.MaxGas
#     else:
#         return Constraint.Skip

# model.max_gas_vol = Constraint(model.T, rule=daily_max_gas)

# constraint to limit consumption to max per 24hr rolling period.  Note, this is *rolling* constraint,
# if a "daily" constraint tied to particular time of day is needed, more work will need to be done
# on the index to identify the indices of interest
window_size = 6  # for testing on small data, normally would be: 24hrs * 12intervals/hr
def rolling_max(model, t):
    preceding_periods = {t_prime for t_prime in model.T if t - window_size <= t_prime < t}
    return sum(model.total_gas[t_prime] for t_prime in preceding_periods) <= model.MaxGas
eval_periods = {t for t in model.T if t >= window_size}
model.rolling_max = Constraint(eval_periods, rule=rolling_max)

# Define the income, expenses, and profit
gas_income = sum(df.loc[t, 'rrp'] * model.gas_for_power[t] / hr  for t in model.T)
gas_expenses = sum(model.total_gas[t] * gprice for t in model.T)  # removed the "vom" computation
profit = gas_income - gas_expenses
model.objective = Objective(expr=profit, sense=maximize)

# Solve the model
solver = SolverFactory('glpk')
results = solver.solve(model)
print(results)
# model.display()
# model.pprint()

cols = []
cols.append(pd.Series(model.P.extract_values(), name='energy price'))
cols.append(pd.Series(model.total_gas.extract_values(), name='total gas'))
cols.append(pd.Series(model.gas_for_power.extract_values(), name='converted gas'))
df_results = pd.concat(cols, axis=1)

fig, ax1 = plt.subplots()
width = 0.25
ax1.bar(np.array(df.index)-width, df_results['energy price'], width, color='g', label='price of power')
ax1.set_ylabel('price')
ax2 = ax1.twinx()
ax2.set_ylabel('gas')
ax2.bar(df.index, df_results['total gas'], width, color='b', label='tot gas')
ax2.bar(np.array(df.index)+width, df_results['converted gas'], width, color='xkcd:orange', label='converted gas')
fig.legend()
fig.tight_layout()
plt.show()