Modelling stable creek intakes

Creek intakes are modelled as tiny reservoirs connected to a tunnel network in SHOP, typically connected directly on the main tunnel of a large plant. Since the reservoir volumes of creek intakes are so small compared to the flow in the tunnel network, it can be challenging to model a stable system when creek intakes are present. This example will describe general rules and best practice for modelling creek intakes in SHOP. Note that this example does not use the old creek_intake object in combination with a junction, as this way of modelling is deprecated.

The system that will be used to illustrate the creek intake modelling in SHOP is based on the code in the function build_model() found in the creek_model.py script:

from pyshop import ShopSession
import pandas as pd
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

from creek_model import build_model

The resulting topology is shown below:

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

A single large reservoir called Reservoir sits at the far end of a long tunnel stretch down to the plant. Three small creek intake reservoirs, Creek1 to 3, are connected to the main tunnel of the plant. All reservoirs have river objects connected to them that acts as spillways with up_head_flow_curve defining the flow out of the reservoirs.

Modelling the reservoir of a creek intake

The reservoirs representing creek intakes should be modelled using the following principles:

Storage capacity

The reservoir of a creek intake represents the volume of the tunnel segment in addition to any small dam at the top of the creek intake. In this model, the max_vol it is set to 0.001 Mm3 for all three creek intakes.

Upper and lower regulation levels

The hrl of the creek intake reservoirs should be the top of the dam or tunnel entrance of the creek intake. Above this level, water will spill out of the creek. The lrl should be set to the lrl of the large reservoir in the network. If the lrl of the creek intake reservoirs are set lower than the lrl of the main reservoir, problems can arise if the large reservoir is empty.

Volume to level conversion

One of the most important factors of modelling a creek intake reservoir is the vol_head curve. Since the reservoir capacity is so small, the level of a creek intake reservoir can change very rapidly over a time step in the optimization model. It is therefore important to model the vol_head curve with the same gradient in all segments. The simple way to do this is to set 0 Mm3 for lrl, max_vol for hrl, and then choose some additional height above hrl and calculate the volume that will ensure a constant gradient. Choosing the level of \(hrl + \Delta\) for the final point on the vol_head curve should result in a volume of \(V = (1+\frac{\Delta}{hrl-lrl})\cdot V_{max}\). This means that the reservoir level can potentially go far above the actual intake of the creek, but this will not happen if the spillway of the creek intake is modelled correctly.

Flow_descr

Before SHOP 16.4.0, the flow_descr attribute on the reservoir object would be used in the tunnel network simulations when a reservoir in the network had a level above hrl. SHOP now only uses the vol_head curve in these calculations, but for older SHOP versions it is important to define a flow_descr on the creek intake reservoirs even though a river is used to model the spillage.

The values for max_vol, hrl, lrl, and the vol_head curve of each reservoir in the model is printed below:

for rsv in shop.model.reservoir:
    print(rsv.get_name())
    print("max_vol:",rsv.max_vol.get())
    print("lrl:",rsv.lrl.get())
    print("hrl:",rsv.hrl.get())
    print("vol_head:")
    print(rsv.vol_head.get())
    print("")

Reservoir
max_vol: 50.0
lrl: 80.0
hrl: 100.0
vol_head:
0.0      80.0
5.0      85.0
12.0     90.0
20.0     95.0
50.0    100.0
60.0    102.0
Name: 0.0, dtype: float64

Creek1
max_vol: 0.001
lrl: 80.0
hrl: 110.0
vol_head:
0.000     80.0
0.001    110.0
0.002    140.0
Name: 0.0, dtype: float64

Creek2
max_vol: 0.001
lrl: 80.0
hrl: 105.0
vol_head:
0.000     80.0
0.001    105.0
0.002    130.0
Name: 0.0, dtype: float64

Creek3
max_vol: 0.001
lrl: 80.0
hrl: 101.0
vol_head:
0.000     80.0
0.001    101.0
0.002    121.0
Name: 0.0, dtype: float64

Modelling the tunnels connecting the creek intakes

The creek intake reservoirs are connected to the main tunnel with short tunnel objects with low loss coefficients. Note that setting the loss coefficients to zero may lead to instabilities when solving the tunnel network equations, so the loss factors should always be larger than zero. In this case, the tunnels connecting the creeks have a loss coefficient of 1e-6, while the main tunnel segments have much higher losses. 1e-6 is about 0.1% of the sum of all loss factors in the tunnel network:

sum_loss_factor = sum(t.loss_factor.get() for t in shop.model.tunnel)

for t in shop.model.tunnel:
    loss = t.loss_factor.get()
    print(f"Loss in {t.get_name()}: {loss}, {100*loss/sum_loss_factor:.2f} % of sum loss factor")
    print("")

Loss in Rsv_N1: 0.00015, 14.96 % of sum loss factor

Loss in N1_N2: 2.5e-05, 2.49 % of sum loss factor

Loss in N2_N3: 2.5e-05, 2.49 % of sum loss factor

Loss in N3_plant: 0.0008, 79.76 % of sum loss factor

Loss in Creek1_N1: 1e-06, 0.10 % of sum loss factor

Loss in Creek2_N2: 1e-06, 0.10 % of sum loss factor

Loss in Creek3_N3: 1e-06, 0.10 % of sum loss factor

Modelling the main tunnel

The main tunnel from the large reservoir to the plant is split into 4 tunnel segments to allow creek intakes to connect at different points in the system. 15% of the main tunnel loss happens in the first tunnel segment between the large reservoir and the first creek intake, then 2.5% loss is added between each creek before 80% of the loss happens after the last creek intake. The main tunnel is modelled as flat at the lrl of the large reservoir before it drops off down to the plant in the final tunnel segment. This is seen in the start_height and end_height of each tunnel:

for t in shop.model.tunnel:
    start = t.start_height.get()
    print(f"{t.get_name()}: starts at {t.start_height.get()} masl, ends at {t.end_height.get()} masl")
    print("")

Rsv_N1: starts at 80.0 masl, ends at 80.0 masl

N1_N2: starts at 80.0 masl, ends at 80.0 masl

N2_N3: starts at 80.0 masl, ends at 80.0 masl

N3_plant: starts at 80.0 masl, ends at 0.0 masl

Creek1_N1: starts at 80.0 masl, ends at 80.0 masl

Creek2_N2: starts at 80.0 masl, ends at 80.0 masl

Creek3_N3: starts at 80.0 masl, ends at 80.0 masl

Modelling the spillways

The overflow example gives a general overview of how to modell spillage in SHOP. For creek intakes, it is recommended to add a river object with a simple up_head_flow_curve to each creek reservoir. The upstream_elevation of the river should be the hrl of the reservoirs, and using a single segment in the up_head_flow_curve is usually sufficient. Setting a flow_cost on the spill rivers is also common practice, but not strictly necessary as long as the creek reservoirs has a water value consistent with the larger reservoir(s) in the tunnel system. In this model, no flow_cost is added to the rivers connected to the creek intake reservoirs, and the energy_value_input of all creeks are set to the same as the large reservoir.

General recommendations and considerations

• Start simple and expand the modelling of creeks in the system gradually. This allows the user to verify that the main components in the system function as expected before adding more complexity.

• Consider aggregating creek intakes that are close together with low losses between them. Large tunnel systems increase the complexity and calculation time of SHOP.

• Make sure the initial reservoir levels of the whole tunnel system are consistent to avoid numerical trouble.

Optimizing a dry system state

It is good practice to test the stability of the modelling in both wet and dry system conditions. First the a dry system state with all reservoirs at lrl and low inflow is optimized:

shop = build_model()

#Set initial state of the reservoirs to empty
for r in shop.model.reservoir:
    r.start_head.set(r.lrl.get())

#Create a random inflow profile that is used for all three creeks
np.random.seed(10)
t = [shop.get_time_resolution()["starttime"] + pd.Timedelta(minutes=60*i) for i in range(24)]
inflow_profile = pd.Series(np.random.rand(24), index=t)

shop.model.reservoir["Creek1"].inflow.set(1*inflow_profile)
shop.model.reservoir["Creek2"].inflow.set(3*inflow_profile)
shop.model.reservoir["Creek3"].inflow.set(2*inflow_profile)

#No inflow on the large reservoir
shop.model.reservoir["Reservoir"].inflow.set(0)

#Run a standard optimization

shop.set_universal_mip("off",[])
shop.start_sim([],[3])
shop.set_universal_mip("on",[])
shop.start_sim([],[2])
shop.set_code(['incremental'], [])
shop.start_sim([],[5])

dry_objective = shop.model.objective.average_objective.total.get()
print("Dry objective value:",dry_objective)

Dry objective value: -2418.8001375127965

As expected, the power plant produces very little since the system is dry. Enough water is saved to produce some power during the price peak towards the end of the optimization horizon.

All reservoir levels slowly increase in the beginning when the power plant is offline, and the inflow into the creeks flow into the larger reservoir. When the power plant draws water, the reservoir levels drop quickly in the closest creek intakes while it takes a little longer for the main upstream reservoir to be affected.

#Plot production and market price
fig = make_subplots(specs=[[{"secondary_y": True}]])
fig.update_layout(barmode="stack", title="Production and market price")
fig.update_yaxes(title_text="Energy price [€/MWh]", secondary_y=True)
fig.update_yaxes(title_text="Production [MW]", secondary_y=False)

for plant in shop.model.plant:

    prod = plant.production.get()
    fig.add_trace(go.Scatter(x=prod.index, y=prod.values, name=plant.get_name()))

price = shop.model.market.spot.sale_price.get()
fig.add_trace(go.Scatter(x=price.index, y=price.values, name="Price"),secondary_y=True)
fig.show()

#Plot reservoir levels
fig = go.Figure()
for rsv in shop.model.reservoir:
    level = rsv.head.get()
    fig.add_trace(go.Scatter(x=level.index, y=level.values, name=rsv.get_name()))

fig.update_layout(title="Reservoir levels")
fig.update_yaxes(title_text="Level [masl]")
fig.show()

#Plot reservoir inflow
fig = go.Figure()
for rsv in shop.model.reservoir:
    inflow = rsv.inflow.get()
    fig.add_trace(go.Scatter(x=inflow.index, y=inflow.values, name=rsv.get_name()))

fig.update_layout(title="Reservoir inflow")
fig.update_yaxes(title_text="Inflow [m3/s]")
fig.show()

The tunnel flow is shown in the figure below. As the reservoir volumes of the creek intakes are so small, most of the inflow goes directly into the larger reservoir. For these calm conditions, there are virtually no deviations between the optimized tunnel flow and the post-calculated physical_flow in the tunnels (note the scale on the y-axis).

fig = go.Figure()
for t in shop.model.tunnel:
    flow = t.flow.get()
    fig.add_trace(go.Scatter(x=flow.index, y=flow.values, name=t.get_name())) 

fig.update_layout(title="Tunnel flow")
fig.update_yaxes(title_text="Flow [m3/s]")
fig.show()    

fig = go.Figure()
for t in shop.model.tunnel:
    flow = t.flow.get()
    pflow = t.physical_flow.get()
    fig.add_trace(go.Scatter(x=flow.index, y=flow.values-pflow.values, name=t.get_name())) 

fig.update_layout(title="Tunnel flow deviation from post-calculated flow")
fig.update_yaxes(title_text="Flow [m3/s]")
fig.show()    

There is no flow in any of the spillways, as expected. As for the tunnel flow, the deviation between the optimized flow and post-calculated physical_flow is zero for all rivers.

#Plot river flow
fig = go.Figure()
for r in shop.model.river:
    flow = r.flow.get()
    fig.add_trace(go.Scatter(x=flow.index, y=flow.values, name=r.get_name()))

fig.update_layout(title="River flow")
fig.update_yaxes(title_text="Flow [m3/s]")
fig.show()

fig = go.Figure()
for r in shop.model.river:
    flow = r.flow.get()
    pflow = r.physical_flow.get()
    fig.add_trace(go.Scatter(x=flow.index, y=flow.values-pflow.values, name=r.get_name()))

fig.update_layout(title="River flow deviation from post-calculated flow")
fig.update_yaxes(title_text="Flow [m3/s]")
fig.show()

Optimizing a wet system state

In the wet system state, all reservoirs start at hrl and the creek intakes receive 10 times the amount of inflow compared to the dry case.

shop = build_model()

#Set initial state of the reservoirs to full
for r in shop.model.reservoir:
    r.start_head.set(r.hrl.get())
        
#Create a random inflow profile that is used for all three creeks
np.random.seed(10)
t = [shop.get_time_resolution()["starttime"] + pd.Timedelta(minutes=60*i) for i in range(24)]
inflow_profile = pd.Series(np.random.rand(24), index=t)

shop.model.reservoir["Creek1"].inflow.set(10*inflow_profile)
shop.model.reservoir["Creek2"].inflow.set(30*inflow_profile)
shop.model.reservoir["Creek3"].inflow.set(20*inflow_profile)

#No inflow on the large reservoir
shop.model.reservoir["Reservoir"].inflow.set(0)

#Run a standard optimization

shop.set_universal_mip("off",[])
shop.start_sim([],[3])
shop.set_universal_mip("on",[])
shop.start_sim([],[2])
shop.set_code(['incremental'], [])
shop.start_sim([],[5])

wet_objective = shop.model.objective.average_objective.total.get()
print("Wet objective value:",wet_objective)

Wet objective value: -521677.9081555805

The power plant has to operate in almost all hours to avoid spillage losses in the system, but is still able to ramp down when the prices are low.

The initial reservoir levels of all creek intakes are high above the large reservoir, but this initial difference is quickly eliminated due to the operation of the power plant and the relatively low losses in the system. The reservoir levels in the three creeks vary rapidly from hour to hour as inflow and power plant discharge changes.

#Plot production and market price
fig = make_subplots(specs=[[{"secondary_y": True}]])
fig.update_layout(barmode="stack", title="Production and market price")
fig.update_yaxes(title_text="Energy price [€/MWh]", secondary_y=True)
fig.update_yaxes(title_text="Production [MW]", secondary_y=False)

for plant in shop.model.plant:

    prod = plant.production.get()
    fig.add_trace(go.Scatter(x=prod.index, y=prod.values, name=plant.get_name()))

price = shop.model.market.spot.sale_price.get()
fig.add_trace(go.Scatter(x=price.index, y=price.values, name="Price"),secondary_y=True)
fig.show()

#Plot reservoir levels
fig = go.Figure()
for rsv in shop.model.reservoir:
    level = rsv.head.get()
    fig.add_trace(go.Scatter(x=level.index, y=level.values, name=rsv.get_name()))

fig.update_layout(title="Reservoir levels")
fig.update_yaxes(title_text="Level [masl]")
fig.show()

#Plot reservoir inflow
fig = go.Figure()
for rsv in shop.model.reservoir:
    inflow = rsv.inflow.get()
    fig.add_trace(go.Scatter(x=inflow.index, y=inflow.values, name=rsv.get_name()))

fig.update_layout(title="Reservoir inflow")
fig.update_yaxes(title_text="Inflow [m3/s]")
fig.show()

The flow in the tunnel system is much more dramatic in the wet system state, and the direction of flow into and out of the large reservoir changes multiple times during the optimization. This is due to the changes in the operation of the power plant and the large but random changes in inflow from hour to hour. There are some slight deviations between the flow and physical_flow in this run, which means that there is a small error in the linearization of the non-linear tunnel network equations. However, the deviations are very small compared to the flow in the tunnels.

#Plot tunnel flow
fig = go.Figure()
for t in shop.model.tunnel:
    flow = t.flow.get()
    fig.add_trace(go.Scatter(x=flow.index, y=flow.values, name=t.get_name()))

fig.update_layout(title="Tunnel flow")
fig.update_yaxes(title_text="Flow [m3/s]")
fig.show()

#Plot tunnel flow
fig = go.Figure()
for t in shop.model.tunnel:
    flow = t.flow.get()
    pflow = t.physical_flow.get()
    fig.add_trace(go.Scatter(x=flow.index, y=flow.values-pflow.values, name=t.get_name()))

fig.update_layout(title="Tunnel flow deviation from post-calculated flow")
fig.update_yaxes(title_text="Flow [m3/s]")
fig.show()

SHOP is able to avoid spillage in the wet system state, but spillage from creek intakes are not uncommon. Specifying a mip_flag (or using the command set universal_river_mip /on) and adding a flow_cost on the creek spillways can be considered when there is a lot of spillage in the system. However, it is important to consider that adding high spillage penalties in a system where spillage is unavoidable can simply make the calculation time longer. This can happen when the objective function is increased a lot due to (spillage) penalties.

#Plot river flow
fig = go.Figure()
for r in shop.model.river:
    flow = r.flow.get()
    fig.add_trace(go.Scatter(x=flow.index, y=flow.values, name=r.get_name()))

fig.update_layout(title="River flow")
fig.update_yaxes(title_text="Flow [m3/s]")
fig.show()

fig = go.Figure()
for r in shop.model.river:
    flow = r.flow.get()
    pflow = r.physical_flow.get()
    fig.add_trace(go.Scatter(x=flow.index, y=flow.values-pflow.values, name=r.get_name()))

fig.update_layout(title="River flow deviation from post-calculated flow")
fig.update_yaxes(title_text="Flow [m3/s]")
fig.show()

creek_model.py

from pyshop import ShopSession
import pandas as pd
import numpy as np

def build_model() -> ShopSession:

    shop = ShopSession()

    shop.set_time_resolution(pd.Timestamp("2024-01-01"), pd.Timestamp("2024-01-02"), "hour") #,"minute",timeresolution=pd.Series([15],index=[pd.Timestamp("2024-01-01")]))

    rsv = shop.model.reservoir.add_object("Reservoir")
    rsv.lrl.set(80)
    rsv.hrl.set(100)
    rsv.max_vol.set(50)
    rsv.vol_head.set(pd.Series([80,85,90,95,100,102], index=[0,5,12,20,50,60]))

    creek1 = shop.model.reservoir.add_object("Creek1")
    creek1.lrl.set(80)
    creek1.hrl.set(110)
    creek1.max_vol.set(0.001)
    creek1.vol_head.set(pd.Series([80,110,140], index=[0,0.001,0.002]))

    creek2 = shop.model.reservoir.add_object("Creek2")
    creek2.lrl.set(80)
    creek2.hrl.set(105)
    creek2.max_vol.set(0.001)
    creek2.vol_head.set(pd.Series([80,105,130], index=[0,0.001,0.002]))

    creek3 = shop.model.reservoir.add_object("Creek3")
    creek3.lrl.set(80)
    creek3.hrl.set(101)
    creek3.max_vol.set(0.001)
    creek3.vol_head.set(pd.Series([80,101,121], index=[0,0.001,0.002]))

    tunnel1 = shop.model.tunnel.add_object("Rsv_N1")
    tunnel1.start_height.set(80)
    tunnel1.end_height.set(80)
    tunnel1.loss_factor.set(0.00015)

    tunnel2 = shop.model.tunnel.add_object("N1_N2")
    tunnel2.start_height.set(80)
    tunnel2.end_height.set(80)
    tunnel2.loss_factor.set(0.000025)

    tunnel3 = shop.model.tunnel.add_object("N2_N3")
    tunnel3.start_height.set(80)
    tunnel3.end_height.set(80)
    tunnel3.loss_factor.set(0.000025)

    tunnel4 = shop.model.tunnel.add_object("N3_plant")
    tunnel4.start_height.set(80)
    tunnel4.end_height.set(0)
    tunnel4.loss_factor.set(0.0008)

    tunnelC1 = shop.model.tunnel.add_object("Creek1_N1")
    tunnelC1.start_height.set(80)
    tunnelC1.end_height.set(80)
    tunnelC1.loss_factor.set(1e-6)

    tunnelC2 = shop.model.tunnel.add_object("Creek2_N2")
    tunnelC2.start_height.set(80)
    tunnelC2.end_height.set(80)
    tunnelC2.loss_factor.set(1e-6)

    tunnelC3 = shop.model.tunnel.add_object("Creek3_N3")
    tunnelC3.start_height.set(80)
    tunnelC3.end_height.set(80)
    tunnelC3.loss_factor.set(1e-6)

    river1 = shop.model.river.add_object("Spill_rsv")
    river1.upstream_elevation.set(100)
    river1.up_head_flow_curve.set([pd.Series([0,50],index=[100,102], name=0)])
    river1.flow_cost.set(1000)

    river2 = shop.model.river.add_object("Spill_C1")
    river2.upstream_elevation.set(110)
    river2.up_head_flow_curve.set([pd.Series([0,50],index=[110,112], name=0)])

    river3 = shop.model.river.add_object("Spill_C2")
    river3.upstream_elevation.set(105)
    river3.up_head_flow_curve.set([pd.Series([0,50],index=[105,107], name=0)])

    river4 = shop.model.river.add_object("Spill_C3")
    river4.upstream_elevation.set(101)
    river4.up_head_flow_curve.set([pd.Series([0,50],index=[101,102], name=0)])

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

    q = [5,7.5,10,12.5,15,17.5,20,22.5,25,27.5,30,32.5,35,37.5,40,42.5,45,47.5,50]
    eff = [80.0, 82.15033783783784, 84.12612612612612, 85.92736486486487, 87.55405405405406, 89.00619369369369, 90.28378378378379, 91.38682432432432, 92.31531531531532, 93.06925675675676, 93.64864864864865, 94.053490990991, 94.28378378378379, 94.33952702702703, 94.22072072072072, 93.92736486486487, 93.45945945945945, 92.81700450450451, 92.0]

    gen1 = shop.model.generator.add_object("G1")
    gen1.penstock.set(1)
    gen1.p_min.set(5)
    gen1.p_max.set(35)
    gen1.p_nom.set(35)
    gen1.gen_eff_curve.set(pd.Series([95,95],index=[5,35]))
    gen1.turb_eff_curves.set([pd.Series(eff,index=q, name=70), pd.Series(np.array(eff)+2,index=q,name=100)])
    gen1.startcost.set(500)

    gen2 = shop.model.generator.add_object("G2")
    gen2.penstock.set(2)
    gen2.p_min.set(5)
    gen2.p_max.set(35)
    gen2.p_nom.set(35)
    gen2.gen_eff_curve.set(pd.Series([95,95],index=[5,35]))
    gen2.turb_eff_curves.set([pd.Series(eff,index=q, name=70), pd.Series(np.array(eff)+2,index=q,name=100)])
    gen2.startcost.set(500)

    rsv.connect_to(tunnel1)
    tunnel1.connect_to(tunnel2)
    tunnel2.connect_to(tunnel3)
    tunnel3.connect_to(tunnel4)
    tunnel4.connect_to(plant)

    creek1.connect_to(tunnelC1)
    tunnelC1.connect_to(tunnel2)

    creek2.connect_to(tunnelC2)
    tunnelC2.connect_to(tunnel3)

    creek3.connect_to(tunnelC3)
    tunnelC3.connect_to(tunnel4)

    rsv.connect_to(river1)
    creek1.connect_to(river2)
    creek2.connect_to(river3)
    creek3.connect_to(river4)


    plant.connect_to(gen1)
    plant.connect_to(gen2)

    market = shop.model.market.add_object("spot")

    t = [pd.Timestamp("2024-01-01") + pd.Timedelta(hours=i) for i in range(24)]
    p = [29.56,28.46,26.66,24.48,24.01,21.23,22.62,25.04,26.24,32.21,41.34,43.51,43.02,44.29,46.24,50.61,59.47,64.99,61.74,55.07,48.01,44.01,45.20,38.00]
    market.sale_price.set(pd.Series(p, index=t))
    market.max_sale.set(10000)

    for r in shop.model.reservoir:
        r.energy_value_input.set(40)

    return shop