Complex ramping on plant discharge

This example shows the types of complex ramping constraints that can be applied to the plant discharge in SHOP under the SHOP_COMPLEX_RAMPING license. Similar constraints can also be applied to reservoir level, as shown in this example. To understand the underlying concepts of rolling and period ramping better, the reader is referred to the complex ramping section in the brief ramping tutorial.

from pyshop import ShopSession
import pandas as pd
import plotly.graph_objs as go
import numpy as np
from complex_ramping import build_model, run_model

Baseline run without ramping constraints

The system used in the example is a single reservoir above a plant with one generator:

shop = build_model()
shop.model.build_connection_tree()

The reservoir starts half full, and has a constant inflow of 10 m3/s incoming during the 24 hour optimization horizon. The market price is above the constant water value in all hours, but the price is also linearly increasing in each hour.

First, we run the model without any complex ramping constraints. The resulting plant discharge is shown below, and we see that SHOP will wait a few hours before producing at best point.

run_model(shop)

plant = shop.model.plant["Plant1"]
discharge_no_ramp = plant.discharge.get().asfreq('1min').ffill()

fig = go.Figure()
fig.add_trace(go.Scatter(x=discharge_no_ramp.index, y=discharge_no_ramp.values, name="Plant discharge"))

fig.update_layout(title="Optimization without ramping constraints", xaxis_title="Time", yaxis_title="m3/s")

fig.show()

Adding complex ramping constraints

Now we will add rolling and period ramping constraints to the plant discharge to limit how much the average, minimum, and maximum discharge can change over certain time windows. First, we create a new SHOP session.

shop = build_model()
plant = shop.model.plant["Plant1"]

res = shop.get_time_resolution()
starttime = res["starttime"]
endtime = res["endtime"]

The plant discharge will be limited downwards by a rolling ramping constraint with a window size of 2 hours (120 min). This is a limit rolling ramping constraint in the downward direction, which means that the plant discharge in every time step is constrained by the maximal discharge value in the rolling window after subtracting the ramping limit.

The limit_discharge_rolling_ramping_down attribute is of type XYT, which makes it possible to change the ramping limit inside the optimization horizon. Note that each XY attribute (in the form of a pd.Series object in the code below) has a reference time stamp. The time stamp specifies from when the ramping constraint is valid. In this example, the ramping limit is 20 m3/s for the first 6 hours of the optimization. It is then set to ‘NaN’, which deactivates the rolling ramping constraint for the remainder of the optimization horizon. Note that only a single rolling ramping constraint with a window size of 2 hours is defined in this example, which means that each XY table has a single point. The reservoir level example shows how to model multiple rolling ramping constraints with different rolling time windows.

rolling_window_sizes = [120]
limit_rolling_ramping_down_1 = pd.Series([20], index = rolling_window_sizes, name=starttime)
limit_rolling_ramping_down_2 = pd.Series([np.nan], index = rolling_window_sizes, name=starttime + pd.Timedelta(hours=6))

plant.limit_discharge_rolling_ramping_down.set([limit_rolling_ramping_down_1, limit_rolling_ramping_down_2])

To limit the discharge in the upward direction, a period ramping constraint with 6 hour (360 min) long periods is added. Unlike the rolling ramping constraint, the period ramping constraint is constraining the average period discharge. This means that the average discharge in each 6 hour period is limited upward by the average discharge in the previous period.

The average_discharge_period_ramping_up attribute is used to add this constraint, and the ramping limit changes from 20 m3/s to 10 m3/s 12 hours into the optimization. Note that changing the ramping limit for period ramping constraints should only be done when a new period begins. Changing it in the middle of a period makes little sense, and an average ramping limit value will be used in those cases.

period_sizes = [360]
average_period_ramping_up_1 = pd.Series([20], index = period_sizes, name=starttime)
average_period_ramping_up_2 = pd.Series([10], index = period_sizes, name=starttime + pd.Timedelta(hours=12))

plant.average_discharge_period_ramping_up.set([average_period_ramping_up_1, average_period_ramping_up_2])

Since the optimization horizon is 24 hours and the ramping periods are 6 hours long, there will be no fractional ramping periods if the first period start at the start time of the optimization. It is not strictly necessary to define an offset for the period ramping when the offset from the start time is zero, but a warning will be printed to the log file. Therefore, the SY attribute average_discharge_period_ramping_up_offset is explicitly defined below, and the offset will be zero since the periods start at the optimization start time:

offset = pd.Series(period_sizes, index = [starttime])
plant.average_discharge_period_ramping_up_offset.set(offset)

Historical discharge values are needed to add both rolling and period ramping constraints at the beginning of the optimization. The period ramping constraint requires the most historical values since the rolling ramping window is only 2 hours long. There is no fractional ramping period at the beginning of the optimization since the specified offset is zero, so it is sufficient to include 6 hours of values to the historical_discharge attribute. However, there is no harm in giving in more historical values than strictly required. 24 hours of historical discharge is randomly generated below:

np.random.seed(5)

#Generate random on/off status of generator
on_off = np.random.randint(0,2,size=24)

#Generate random discharge between 25 and 100
discharge_level = 25 + 75*np.random.rand(24)

discharge = on_off * discharge_level
t = [starttime - pd.Timedelta(hours=24-i) for i in range(24)]

historical_discharge = pd.Series(discharge, index=t)

plant.historical_discharge.set(historical_discharge)

historical_discharge = historical_discharge.asfreq('1min').ffill()

fig = go.Figure()
fig.add_trace(go.Scatter(x=historical_discharge.index, y=historical_discharge.values, name="Historical discharge"))
fig.update_layout(title="Historical plant discharge values", xaxis_title="Time", yaxis_title="m3/s")

fig.show()

Comparing results

Now we can optimize and compare the optimization results to the initial SHOP run without ramping constraints.

run_model(shop)

As seen in the figure below, the plant discharge can’t ramp up from zero to 90 m3/s in the way the unconstrained solution suggests. Note that the historical discharge values have been appended to the optimized discharge results with ramping.

discharge = plant.discharge.get()

total_discharge = pd.concat((historical_discharge,discharge))
total_discharge = total_discharge.asfreq('1min').ffill()

fig = go.Figure()
fig.add_trace(go.Scatter(x=total_discharge.index, y=total_discharge.values, name="With ramping"))
fig.add_trace(go.Scatter(x=discharge_no_ramp.index, y=discharge_no_ramp.values, name="Without ramping"))

fig.update_layout(title="Plant discharge comparison", xaxis_title="Time", yaxis_title="m3/s")

fig.show()

To validate that the complex ramping constraints are not violated, we will post-calculate the rolling and period discharge values and their implied ramping bounds.

Period ramping validation

The figure below shows the plant discharge together with the average discharge in all ramping periods, including the purely historical ramping period before the optimization starts. The upper ramping limit for the average discharge in a period is calculated by adding 20 m3/s (or 10 m3/s after 12 hours) to the average of the previous period. The vertical dotted lines show the beginning of each fixed ramping period.

#Relevant discharge for period ramping
total_discharge_period = total_discharge[total_discharge.index >= starttime - pd.Timedelta(hours = 6)]

historic = total_discharge_period[total_discharge_period.index <= starttime]
optimized = total_discharge_period[total_discharge_period.index >= starttime]

#Plot discharge for historic and optimization periods
fig = go.Figure()
fig.add_trace(go.Scatter(x=historic.index, y=historic.values, name="Historical plant discharge"))
fig.add_trace(go.Scatter(x=optimized.index, y=optimized.values, name="Optimized plant discharge"))

t_mean = []
mean_period_dis = []

t_bound = []
mean_upper_bound = []

for i in range(5):
    #The start and end time of each fixed ramping period
    period_start = starttime + pd.Timedelta(hours = 6*(i - 1))
    period_end = period_start + pd.Timedelta(hours = 6)
    
    #Get the discharge in the current ramping period, calculate the average
    period_dis = total_discharge_period[total_discharge_period.index >= period_start]
    period_dis = period_dis[period_dis.index <= period_end]
    average_dis = period_dis.mean()
    
    t_mean += [period_start, period_end]    
    mean_period_dis += [average_dis, average_dis]  
    
    #Calculate the upper ramping limit for the next period based on the average in this period
    if period_end < endtime:
        t_bound += [period_end, period_end + pd.Timedelta(hours = 6)]
        
        limit = 20
        if period_end >= starttime + pd.Timedelta(hours = 12):
            limit = 10
        mean_upper_bound += [average_dis + limit, average_dis + limit]            
    
    
    #Plot vertical dashed lines to indicate the change between two periods
    fig.add_trace(go.Scatter(x=[period_start, period_start], y=[0, 95],mode="lines",showlegend=False, line = dict(color = 'rgb(0, 0, 0 )', dash = 'dot')))

#Plot the average discharge in each period
fig.add_trace(go.Scatter(x=t_mean, y=mean_period_dis, name="Average period discharge"))
#Plot the upper limit for the average discharge in each period
fig.add_trace(go.Scatter(x=t_bound, y=mean_upper_bound, name="Upper average discharge limit",line = dict(color = 'rgb(0, 0, 0)',dash = 'dash')))

fig.update_layout(title="Visualizing the upward period ramping constraint", xaxis_title="Time", yaxis_title="m3/s")

fig.show()

The average in the first purely historical period adds an upper limit for the average discharge in the first ramping period within the optimization interval. Note that the discharge in individual hours can be above the plotted upper limit, while the average in each period stays within the bound. We see that the upward average period ramping constraint is only binding in two of the ramping periods.

Rolling ramping validation

The downward limit rolling ramping constraint is visualized in the figure below. The plant discharge is shown together with the post-calculated rolling maximal discharge value with a 2 hour window. The lower discharge limit is calculated by subtracting 20 m3/s from the rolling max for the first 6 hours, while it is set to zero afterwards since the constraint is deactivated.

#Relevant plant discharge for 2 hour rolling ramping
total_dis_rolling = total_discharge[total_discharge.index >= starttime - pd.Timedelta(minutes=120)]

historic = total_dis_rolling[total_dis_rolling.index <= starttime]
optimized = total_dis_rolling[total_dis_rolling.index >= starttime]

#Plot discharge for historic and optimization periods
fig = go.Figure()
fig.add_trace(go.Scatter(x=historic.index, y=historic.values, name="Historical plant discharge"))
fig.add_trace(go.Scatter(x=optimized.index, y=optimized.values, name="Optimized plant discharge"))

#Upsample the discharge time series to 1 second frequency and calculate the rolling average
total_dis_rolling = total_dis_rolling.asfreq('1s').ffill()
rolling_max = total_dis_rolling.rolling(pd.Timedelta(minutes=120)).max()
rolling_max = rolling_max[rolling_max.index >= starttime]

#Constraint is deactivated for last 12 hours
lower_limit = rolling_max - 20
lower_limit[lower_limit.index >= starttime + pd.Timedelta(hours = 6)] = 0

#Plot rolling 6 hour maximum and lower discharge due to the rolling max
fig.add_trace(go.Scatter(x=rolling_max.index, y=rolling_max.values, name="Maximal rolling discharge, 2h window"))        
fig.add_trace(go.Scatter(x=lower_limit.index, y=lower_limit.values, name="Lower discharge limit", line = dict(color = 'rgb(0, 0, 0)', dash = 'dash')))        

fig.update_layout(title="Visualizing the downward rolling ramping constraint with a 2 hour window",xaxis_title="Time", yaxis_title="masl")

fig.show()

We see that the rolling ramping constraint is binding in the first hours of the optimization where the discharge is at the lower bound. As soon as the constraint is deactivated, the plant discharge drops rapidly in one hour to accommodate the upward period ramping constraint while operating at higher efficiency for longer.

complex_ramping.py

from pyshop import ShopSession
import pandas as pd

def build_model() -> ShopSession:

    shop = ShopSession()
    starttime=pd.Timestamp("2023-01-01 00:00:00")
    endtime=pd.Timestamp("2023-01-02 00:00:00")
    shop.set_time_resolution(starttime=starttime, endtime=endtime, timeunit="minute", timeresolution=pd.Series([60],index=[starttime]))

    rsv1=shop.model.reservoir.add_object("Reservoir1")
    rsv1.hrl.set(100)
    rsv1.lrl.set(90)
    rsv1.max_vol.set(12)
    rsv1.vol_head.set(pd.Series([90,100,101],index=[0,12,14]))

    plant1=shop.model.plant.add_object("Plant1")
    plant1.main_loss.set([0.0002])
    plant1.penstock_loss.set([0.0001])
    plant1.outlet_line.set(0)

    p1g1=shop.model.generator.add_object("P1G1")
    p1g1.connect_to(plant1)
    p1g1.penstock.set(1)
    p1g1.p_min.set(10)
    p1g1.p_nom.set(100)
    p1g1.p_max.set(100)
    p1g1.startcost.set(500)
    p1g1.gen_eff_curve.set(pd.Series([95,98], index=[0,100]))
    p1g1.turb_eff_curves.set([pd.Series([80,95,90],index=[25,90,100],name=90),pd.Series([82,98,92],index=[25,90,100],name=100)])

    rsv1.connect_to(plant1)
   
    rsv1.start_head.set(95)
    rsv1.inflow.set(10)
    rsv1.energy_value_input.set(40)

    da=shop.model.market.add_object('Day_Ahead')
    da.sale_price.set(pd.Series([40 + i for i in range(24)],index=[starttime+pd.Timedelta(hours=i) for i in range(24)]))
    da.buy_price.set(da.sale_price.get()+0.1)
    da.max_sale.set(9999)
    da.max_buy.set(9999)

    gs = shop.model.global_settings.global_settings
    gs.plant_discharge_ramp_penalty_cost.set(1000000)

    return shop
 
def run_model(shop:ShopSession) -> None:
    
    shop.set_universal_mip("on",[])
    shop.start_sim([],3)
    shop.set_code("incremental",[])
    shop.start_sim([],5)