Complex ramping on reservoir level

This example shows the types of complex ramping constraints that can be applied to the reservoir level in SHOP under the SHOP_COMPLEX_RAMPING license. Similar constraints can also be applied to plant discharge, 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 reservoir level is shown below, and we see that SHOP will wait a few hours before producing heavily and almost emptying the reservoir.

run_model(shop)

rsv = shop.model.reservoir["Reservoir1"]
level_no_ramp = rsv.head.get()
lrl = rsv.lrl.get()

fig = go.Figure()
fig.add_trace(go.Scatter(x=level_no_ramp.index, y=level_no_ramp.values, name="Reservoir level"))
fig.add_trace(go.Scatter(x=level_no_ramp.index, y=[lrl]*len(level_no_ramp),name="LRL",line = dict(color = 'rgb(0, 0, 0)', dash = 'dash')))

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

fig.show()

Adding complex ramping constraints

Now we will add rolling and period ramping constraints to the reservoir level to limit how much the reservoir level can change over certain time windows. First, we create a new SHOP session:

shop = build_model()
rsv = shop.model.reservoir["Reservoir1"]

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

The reservoir level will be constrained upwards by two rolling ramping constraints with different rolling window lengths. One has a 6 hour (360 min) rolling window and the other has a 12 hour (720 min) rolling time window. The ramping limits are 4 cm and 5 cm, respectively, and this limit is calculated upward from the average level in the rolling time windows. The appropriate attribute to use for these constraints is therefore average_level_rolling_ramping_up. This is an XYT attribute, which makes it possible to specify several rolling ramping constraints with different window sizes. An XY table (represented by a pd.Series object in the code below) is defined where the window sizes in minutes are the x-values and the corresponding ramping limits are the y-values. Since this is an XYT attribute the XY table is valid from the specified time stamp, which is the start time of the optimization in this case. The ramping limits don’t change during the optimization horizon in this example, so a single XY table is used in the XYT. The plant discharge example shows how the XYT can be defined when the ramping limits change.

rolling_window_sizes = [360, 720]
ramping_limits = [0.04,0.05]
average_rolling_ramping_up = pd.Series(ramping_limits, index = rolling_window_sizes, name=starttime)

rsv.average_level_rolling_ramping_up.set([average_rolling_ramping_up])

For downwards ramping, a period ramping constraint with 10 hour (600 min) long periods is added. This periodic ramping constraint limits the change from the highest reservoir level in one period to the lowest reservoir level in the next period to 75 cm. Since the ramping is applied between extreme values, from max to min in downward ramping and from min to max in upward ramping, the constraint is called limit period ramping. The limit_level_period_ramping_down attribute is defined in a similar way to the rolling ramping constraints. It is possible to specify multiple period constraints with different period lengths, but a single period constraint is used in this example.

period_sizes = [600]
ramping_limits = [0.75]
limit_period_ramping_down = pd.Series(ramping_limits, index = period_sizes, name=starttime)

rsv.limit_level_period_ramping_down.set([limit_period_ramping_down])

The default behaviour in SHOP is to assume that the first fixed ramping period starts at the beginning of the optimization horizon. Since the ramping periods are 10 hours long and the optimization horizon is 24 hours, the four last hours of the optimization will therefore be part of an incomplete ramping period. No constraints will be built from the last full ramping period to the final ramping period if it is too short. However, it is also possible to specify an offset to signal that the first full ramping period should start after the start time of the optimization. This is done below with the SY attribute limit_level_period_ramping_down_offset so that the 4 first hours are fractional instead of the last four:

offset = pd.Series(period_sizes, index = [starttime + pd.Timedelta(hours=4)])
rsv.limit_level_period_ramping_down_offset.set(offset)

Historical reservoir level values should be supplied to SHOP when running with complex ramping constraints since both rolling and periodic constraints look backwards in time. Historical values will therefore allow SHOP to add proper ramping constraints for the first hours of the optimization, and can be defined using the TXY attribute historical_level.

For simplicity, hourly values for the first day before the optimization start time is used in this example. To cover all rolling ramping constraints in this example, it is sufficient to define historical reservoir levels that begin 12 hours before optimization start. Historical values must be defined for 16 hours before the optimization start time to properly initialize the periodic ramping constraints, 6 hours to complete the initial fractional period and another 10 hours to add a full and purely historical ramping period before that.

np.random.seed(0)

#We generate 24 historic hourly reservoir values by going on a random walk backwards from the start level 95 masl
#Each value is within 5 cm of the previous one
prev_level = 95
level = []
for i in range(24):
    l = prev_level + 0.05*(1-2*np.random.rand())
    level.append(l)
    prev_level = l
    
#Flip values so that the oldest comes first
level = level[::-1] 
t = [starttime - pd.Timedelta(hours=24-i) for i in range(24)]
historical_level = pd.Series(level, index=t)

rsv.historical_level.set(historical_level)

fig = go.Figure()
fig.add_trace(go.Scatter(x=historical_level.index, y=historical_level.values, name="Historical reservoir level"))
fig.update_layout(title="Historical reservoir level values", xaxis_title="Time", yaxis_title="masl")

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, it is no longer possible to drain the reservoir without violating the ramping constraints. Note that the historical reservoir levels have been appended to the optimized reservoir level results.

level = rsv.head.get()

#Combine historical and optimized reservoir levels to one Series
total_level = pd.concat((historical_level,level))
total_level = total_level.asfreq('1min').interpolate()

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

fig.update_layout(title="Reservoir level comparison", xaxis_title="Time", yaxis_title="masl")

fig.show()

To further visualize the different ramping constraints, we will validate that they are not broken by post-calculating the average reservoir levels in the fixed ramping periods and as a rolling average.

Period ramping validation

The figure below includes the reservoir level in addition to the maximum reservoir level reached within each fixed ramping period. The lower reservoir limit due to the ramping constraint is also calculated and shown as a black dashed line. The vertical dotted lines show the beginning of each fixed ramping period.

#Relevant reservoir levels for period ramping
total_level_period = total_level[total_level.index >= starttime - pd.Timedelta(hours = 16)]

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

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

t_max = []
max_period_level = []

for i in range(4):
    #The start and end time of each fixed ramping period
    period_start = starttime + pd.Timedelta(hours=-16 + 10*i)
    period_end = period_start + pd.Timedelta(hours=10)
    
    #Get the reservoir level in the ramping period, calculate the maximal level
    period_level = total_level_period[total_level_period.index >= period_start]
    period_level = period_level[period_level.index <= period_end]
    max_level = period_level.max()
    
    t_max += [period_start,period_end]    
    max_period_level += [max_level, max_level]     
    
    #Plot vertical dashed lines to indicate the change between two periods
    fig.add_trace(go.Scatter(x=[period_start,period_start], y=[93.5,95.5],mode="lines",showlegend=False, line = dict(color = 'rgb(0, 0, 0 )', dash = 'dot')))
    
t_max = np.array(t_max)   
max_period_level = np.array(max_period_level)

#The lower reservoir boundary is the max level in the previous period minus 75 cm
#Skip plotting lower level for times after the optimization
t_min = t_max[:-2] + pd.Timedelta(hours=10)
min_rsv_bound = max_period_level[:-2] - 0.75
    
#Plot the average level in each period
fig.add_trace(go.Scatter(x=t_max, y=max_period_level, name="Maximal period level"))
#Plot the lower limit for the average level in each period
fig.add_trace(go.Scatter(x=t_min, y=min_rsv_bound, name="Lower reservoir limit",line = dict(color = 'rgb(0, 0, 0)',dash = 'dash')))

fig.update_layout(title="Visualizing the downward period ramping constraints", xaxis_title="Time", yaxis_title="masl")

fig.show()

It is clear from the figure why at least 16 hours of historical values are needed to add periodic ramping constraints of 10 hour length in this case. The first visualized ramping period consists of purely historical values and can’t be optimized by SHOP. But the maximal value in this period limits the minimal value in the next ramping period, which includes the 4 first hours of the optimization. The remainder of the optimization is covered by two 10 hour ramping periods.

The lower limit due to the limit period ramping constraints is shown as a green dashed line, and calculated by taking the maximal reservoir value in the previous period and subtracting the ramping limit of 75 cm. The minimal reservoir level in the two last periods reach this lower limit, but it is never violated. Note that the start and end times of each ramping period is included when calculating the min and max reservoir level values in SHOP.

Rolling ramping validation

The average rolling reservoir level with a window size of 6 hours is calculated and shown in the figure below. The upper reservoir level limit due to this rolling average is also shown as a black dashed line:

#Relevant reservoir levels for 6 hour rolling ramping
total_level_rolling = total_level[total_level.index >= starttime - pd.Timedelta(minutes=360)]

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

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

#Upsample the reservoir level time series to 1 second frequency and calculate the rolling average
total_level_rolling = total_level_rolling.asfreq('1s').interpolate()
rolling_average = total_level_rolling.rolling(pd.Timedelta(minutes=360)).mean()
rolling_average = rolling_average[rolling_average.index >= starttime]

#Plot rolling 6 hour average and upper reservoir limit due to the rolling average
fig.add_trace(go.Scatter(x=rolling_average.index, y=rolling_average.values, name="Average rolling level, 6h window"))        
fig.add_trace(go.Scatter(x=rolling_average.index, y=rolling_average.values + 0.04, name="Upper reservoir limit", line = dict(color = 'rgb(0, 0, 0)', dash = 'dash')))        

fig.update_layout(title="Visualizing the upward rolling ramping constraints with a 6 hour window",xaxis_title="Time", yaxis_title="masl")

fig.show()

The optimized reservoir level always stays within the upper bound due to the rolling ramping constraint, and the constraint is binding for the 3 local maximums of the reservoir level. Note that the rolling average is post-calculated in this example by taking a simple average of the reservoir level time series with a rolling window. This is not identical to the true average since the reservoir level is assumed to change linearly over time steps. However, when the time series is upsampled to a frequency of 1 second, the difference between the true average found by integration and the simple average is vanishing. In SHOP, the true rolling average is used when building the ramping constraints in the optimization problem.

The average rolling ramping constraint with a 12 hour window is shown below:

#Relevant reservoir levels for 12 hour rolling ramping
total_level_rolling = total_level[total_level.index >= starttime - pd.Timedelta(minutes=720)]

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

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

#Upsample the reservoir level time series to 1 second frequency and calculate the rolling average
total_level_rolling = total_level_rolling.asfreq('1s').interpolate()
rolling_average = total_level_rolling.rolling(pd.Timedelta(minutes=720)).mean()
rolling_average = rolling_average[rolling_average.index >= starttime]

#Plot rolling 12 hour average and upper reservoir limit due to the rolling average
fig.add_trace(go.Scatter(x=rolling_average.index, y=rolling_average.values, name="Average rolling level, 12h window"))
fig.add_trace(go.Scatter(x=rolling_average.index, y=rolling_average.values + 0.05, name="Upper reservoir limit", line = dict(color = 'rgb(0, 0, 0)', dash = 'dash')))

fig.update_layout(title="Visualizing the upward rolling ramping constraints with a 12 hour window", xaxis_title="Time", yaxis_title="masl")

fig.show()

The 12 hour average changes more slowly than the 6 hour average, and is only binding at one point in time since the reservoir trajectory is mostly downward. Note that the shown upper reservoir limit is broken at the very beginning of the optimization. This is not actually the case in the optimization since the initial reservoir level is fixed by the user input, and so no rolling ramping constraints are added for t = 0. The optimization has hourly time steps, and the first optimized reservoir level at t = 1 is within the rolling ramping constraint. The post-calculated rolling average has a resolution of 1 second, and shows that the reservoir trajectory would violate the ramping constraint for the first 15 minutes or so. A finer time resolution in SHOP would remedy much of the initial problems, but the randomly generated historical values would violate the constraint around t = 0 in any case.

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.level_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)