---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.2
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

(creek-intake-example)=
# Modelling stable creek intakes

Creek intakes are modelled as tiny [reservoirs](reservoir) connected to a tunnel network in SHOP, typically connected directly on the main tunnel of a large [plant](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](creek_intake) object in combination with a [junction](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:

```{code-cell} ipython3
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:

```{code-cell} ipython3
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](river) objects connected to them that acts as spillways with [up_head_flow_curve](river: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](reservoir:max_vol) it is set to 0.001 Mm3 for all three creek intakes.

### Upper and lower regulation levels
The [hrl](reservoir: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](reservoir: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](reservoir:vol_head). 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](reservoir: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:

```{code-cell} ipython3
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("")
```

## 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](tunnel:loss_factor). 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:

```{code-cell} ipython3
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("")
```

## 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](tunnel:start_height) and [end_height](tunnel:end_height) of each tunnel:

```{code-cell} ipython3
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("")
```

## Modelling the spillways

The [overflow example](reservoir-overflow) 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](river:up_head_flow_curve) to each creek reservoir. The [upstream_elevation](river: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](river:flow_cost) on the spill rivers is also common practice, but not strictly necessary as long as the creek reservoirs has a [water value](water-value-descriptions) 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](reservoir: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:

```{code-cell} ipython3
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)
```

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.

```{code-cell} ipython3
#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](tunnel:flow) and the post-calculated [physical_flow](tunnel:physical_flow) in the tunnels (note the scale on the y-axis).

```{code-cell} ipython3
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](river:flow) and post-calculated [physical_flow](river:physical_flow) is zero for all rivers.

```{code-cell} ipython3
#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.

```{code-cell} ipython3
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)
```

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.

```{code-cell} ipython3
#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.

```{code-cell} ipython3
#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](river:mip_flag) (or using the command [set universal_river_mip /on](set_universal_river_mip)) and adding a [flow_cost](river: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.

```{code-cell} ipython3
#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

```{code-cell} ipython3
:Collapsed: 'false'
:tags: [remove-input]

with open('creek_model.py', 'r') as f:
    print(f.read())
```