Time delay and environmental constraints

This example shows how rivers can be used to model environmental constraints in combination with advanced time delay functionality, natural flow out of reservoirs, and controlled discharge from power plants.

from advanced_river import build_model, run_commands
import pandas as pd
import plotly.graph_objs as go

The model definition is set up in the build_model() function in advanced_river.py. The time resolution of the model is 15 minutes.

shop = build_model()
time = shop.get_time_resolution()
starttime = time["starttime"]

The watercourse consists of 4 rivers, 2 reservoirs, and a single power plant. The rivers are connected in a cascade as shown in the figure below.

shop.model.build_connection_tree()

Input data

We will now specify additional input data on the 4 river objects.

r_inflow = shop.model.river.river_w_inflow
r_min_flow = shop.model.river.river_w_min_flow
r_flow_curve = shop.model.river.river_w_flow_curve
r_delay_curve = shop.model.river.river_w_delay_curve

Inflow is added to the r_inflow object, and a constant time delay of 3.5 hours is specified with the time_delay_const attribute. past_upstream_flow is also added to capture flow on its way down the river when the optimization starts. This past upstream flow will flow out of the bottom of the river for the first few hours before the new inflow arrives downstream after 3.5 hours.

r_inflow.inflow.set(10)
r_inflow.time_delay_const.set(3.5)
#Past time is measured in hours before optimization start, which means negative numbers
past_time = [starttime-pd.Timedelta(hours=10),starttime-pd.Timedelta(hours=2),starttime-pd.Timedelta(hours=1),starttime]
past_inflow = [15,14,13,12]
r_inflow.past_upstream_flow.set(pd.Series(past_inflow, index=past_time))

The minimum flow requirement of 100 m3/s is added to the r_min_flow object using the min_flow attribute. The requirement is only active between 72 and 96 hours into the optimization

min_flow = [0,100,0]
times = [starttime,starttime+pd.Timedelta(hours=72),starttime+pd.Timedelta(hours=96)]

r_min_flow.min_flow.set(pd.Series(min_flow,index=times))        

The flow out of the unregulated reservoir into the river object r_flow_curve below is determined by the XYN curve up_head_flow_curve. The x-values specify the upstream reservoir head, and the y-values is the resulting flow into the river. Several XY curves with reference values describing different gate positions may be provided in the up_head_flow_curve, but there is no controllable gate between the reservoir and the river in this example. Only a single XY curve is given (the reference value of 0 is unused in the optimization), which represents water flowing over a weir of constant height equal to the upstream_elevation of the river.

head = [100,102,104,106,108,110,111]
flow = [0,15,20,30,60,250,1000]

r_flow_curve.up_head_flow_curve.set([pd.Series(flow, index=head, name=0)])

fig = go.Figure()
fig.add_trace(go.Scatter(x=head, y=flow, name="Flow curve"))
fig.update_layout(xaxis_title="Reservoir level [masl]", yaxis_title="Flow [m3/s]")
fig.show()

Finally, we add a more complex time delay representation to the final river object r_delay_curve below the plant. The time_delay_curve is constructed to be exponentially decreasing over time. The first point is 3 hours after the water is discharged, while the last point is 13 hours after. Due to the exponential decay of the factors, most of the water will start arriving downstream after roughly 3 hours. Note that the final flow factor must be zero and that we scale the flow factors so that the sum is equal to 1. The scaling is also done automatically by SHOP, but a warning will be printed to the log file if the input data is not properly normalized.

time_delay = [3+i for i in range(11)]
flow_factor = [2**(-i) for i in range(1,11)]+[0]
flow_factor = [f/sum(flow_factor) for f in flow_factor]

r_delay_curve.time_delay_curve.set([pd.Series(flow_factor,index=time_delay,name=0)])

fig = go.Figure()
fig.add_trace(go.Scatter(x=time_delay, y=flow_factor, name="Flow fractions"))
fig.update_layout(xaxis_title="Time [hours]", yaxis_title="Flow fraction")
fig.show()

Now we can run the optimization and look at the results

run_commands(shop)

Results

First, we check that the flow in the bottom river meets the minimum flow requirement that we specified:

flow = r_min_flow.flow.get()
min_flow = r_min_flow.min_flow.get()

fig = go.Figure()
fig.add_trace(go.Scatter(x=flow.index, y=flow.values, name="Flow"))
fig.add_trace(go.Scatter(x=min_flow.index, y=min_flow.values, name="Minimum flow requirement"))
fig.update_layout(xaxis_title="Time", yaxis_title="Flow [m3/s]")
fig.show()

The flow in the bottom river is the combined flow of all upstream rivers. It turns out that most of the water flowing in the bottom river comes from the natural flow out of the unregulated reservoir, which is determined by the reservoir level and the input up_head_flow_curve on the r_flow_curve river object. Notice that distinctive kinks in the river flow appears when the reservoir level reaches 104 masl and 106 masl since the up_head_flow_curve changes for these values.

from plotly.subplots import make_subplots

flow = r_flow_curve.flow.get()
head = shop.model.reservoir.rsv_unreg.head.get()

fig = make_subplots(specs=[[{"secondary_y": True}]])
fig.add_trace(go.Scatter(x=flow.index, y=flow.values, name="Flow"))

fig.add_trace(go.Scatter(x=head.index, y=head.values, name="Reservoir level"),secondary_y=True)
fig.update_layout(xaxis_title="Time", yaxis_title="Flow [m3/s]")
fig.update_yaxes(title_text="Reservoir level [masl]", secondary_y=True)
fig.show()

The water coming from the inflow river is also unregulated, and the upstream river flow is equal to the inflow for all time steps. Unlike the flow in the river connected to the unregulated reservoir, there is a time delay in the inflow river. This means that the downstream_flow is different from the upstream flow for the first hours due to the arrival of the past_upstream_flow.

flow = r_inflow.flow.get()
flow_down = r_inflow.downstream_flow.get()
inflow = r_inflow.inflow.get()

fig = go.Figure()
fig.add_trace(go.Scatter(x=flow.index, y=flow.values, name="Upstream flow"))
fig.add_trace(go.Scatter(x=flow_down.index, y=flow_down.values, name="Downstream flow"))
fig.add_trace(go.Scatter(x=inflow.index, y=inflow.values, name="Natural inflow"))
fig.update_layout(xaxis_title="Time", yaxis_title="Flow [m3/s]")
fig.show()

The combined unregulated flow from the reservoir and the inflow river is not enough to satisfy the minimum flow constraint. The only controllable flow in this example comes from the discharge of the power plant into the river with a time_delay_curve, and it is needed to cover the last part of the environmental constraint. The water value in the regulated reservoir is high compared to the market price in this example, and the plant only produces power to help cover the minimum flow constraint in the bottom river.

The upstream flow in the river is the same as the plant discharge in this case since there is no inflow or other upstream connections to the river. The upstream flow time series has a rectangular box shape in this case, and the flow ramps up and down between 0 and about 50 m3/s in a single time step. The ramping in the downstream flow is smoothed out due to the exponential decay used in the time_delay_curve. We note that plant has to time its discharge so that the water will arrive in time to cover the minimum flow constraint downstream.

flow = r_delay_curve.flow.get()
flow_down = r_delay_curve.downstream_flow.get()

fig = go.Figure()
fig.add_trace(go.Scatter(x=flow.index, y=flow.values, name="Upstream flow"))
fig.add_trace(go.Scatter(x=flow_down.index, y=flow_down.values, name="Downstream flow"))

fig.update_layout(xaxis_title="Time", yaxis_title="Flow [m3/s]")
fig.show()

advanced_river.py

from pyshop import ShopSession
import pandas as pd
from pandas import Timestamp
from numpy import nan

def build_model() -> ShopSession:

    #Initialize a new ShopSession
    shop = ShopSession()

    #Set the time resolution of the optimization
    start_time = Timestamp('2018-02-27 00:00:00')
    end_time = Timestamp('2018-03-05 00:00:00')
    t = [Timestamp('2018-02-27 00:00:00')]
    y = [15.0]
    step_length = pd.Series(y,index=t)
    shop.set_time_resolution(start_time, end_time, 'minute', step_length)

    #Add all objects and set all attributes
    river_w_flow_curve = shop.model.river.add_object('river_w_flow_curve')
    river_w_flow_curve.upstream_elevation.set(100.0)

    river_w_delay_curve = shop.model.river.add_object('river_w_delay_curve')
    river_w_delay_curve.upstream_elevation.set(60.0)

    river_w_min_flow = shop.model.river.add_object('river_w_min_flow')
    river_w_min_flow.upstream_elevation.set(50.0)

    river_w_inflow = shop.model.river.add_object('river_w_inflow')
    river_w_inflow.upstream_elevation.set(60.0)

    rsv_unreg = shop.model.reservoir.add_object('rsv_unreg')
    rsv_unreg.max_vol.set(12.0)
    rsv_unreg.lrl.set(100.0)
    rsv_unreg.hrl.set(110.0)
    x = [0.0, 12.0, 14.0]
    y = [100.0, 110.0, 111.0]
    rsv_unreg.vol_head.set(pd.Series(y,index=x,name=0.0))
    rsv_unreg.start_head.set(102.0)
    rsv_unreg.inflow.set(50.0)
    rsv_unreg.energy_value_input.set(40.0)

    rsv_reg = shop.model.reservoir.add_object('rsv_reg')
    rsv_reg.max_vol.set(120.0)
    rsv_reg.lrl.set(100.0)
    rsv_reg.hrl.set(110.0)
    x = [0.0, 120.0, 124.0]
    y = [100.0, 110.0, 111.0]
    rsv_reg.vol_head.set(pd.Series(y,index=x,name=0.0))
    rsv_reg.start_head.set(102.0)
    rsv_reg.inflow.set(0.0)
    rsv_reg.energy_value_input.set(41.0)

    plant = shop.model.plant.add_object('plant')
    plant.prod_area.set(1.0)
    plant.outlet_line.set(60.0)
    plant.main_loss.set([0.0002])
    plant.penstock_loss.set([0.0001])

    market_day_ahead = shop.model.market.add_object('Day_ahead')
    market_day_ahead.market_type.set('ENERGY')
    market_day_ahead.max_buy.set(9999.0)
    market_day_ahead.max_sale.set(9999.0)
    market_day_ahead.buy_price.set(40.01)
    market_day_ahead.sale_price.set(39.99)

    g1 = shop.model.generator.add_object('plant_G1')
    g1.penstock.set(1)
    g1.p_min.set(5.0)
    g1.p_max.set(20.0)
    g1.p_nom.set(20.0)
    x = [0.0, 20.0]
    y = [95.0, 98.0]
    g1.gen_eff_curve.set(pd.Series(y,index=x,name=0.0))
    xy_curves = []
    x = [15.0, 48.3, 57.0]
    y = [80.0, 95.0, 90.0]
    xy_curves.append(pd.Series(y,index=x,name=40.0))
    x = [12.0, 37.5, 45.0]
    y = [82.0, 98.0, 92.0]
    xy_curves.append(pd.Series(y,index=x,name=50.0))
    g1.turb_eff_curves.set(xy_curves)
    g1.startcost.set(500.0)


    #Connect all objects
    river_w_flow_curve.connect_to(river_w_min_flow)

    river_w_delay_curve.connect_to(river_w_min_flow)

    river_w_inflow.connect_to(river_w_min_flow)

    rsv_unreg.connect_to(river_w_flow_curve)

    rsv_reg.connect_to(plant)

    plant.connect_to(river_w_delay_curve)
    plant.connect_to(g1)

    return shop

def run_commands(shop:ShopSession) -> None:
    
    shop.set_universal_mip("off",[])
    shop.start_sim([],4)
    shop.set_universal_mip("on",[])
    shop.start_sim([],2)
    shop.set_code("incremental",[])
    shop.start_sim([],5)