---
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
---

(sim-prod-plan-example)=
# Simulating deviations from a production plan
This documentation demonstrates how to use the SHOP simulator to simulate the validity of an optimized production plan and how changes to the production plan will affect the physical hydropower system. This could be useful when ensuring that a bid in an intra-day or real-time market is physically feasible without re-running the whole optimization. 

The model setup for this example is available in the following format:

- pyshop
    - [](sim_two_rsv.py)

## System topology and initial simulation
First, a simple SHOP model with two [reservoirs](reservoir) and two [plants](plant) are built. In addition, a reserve obligation for spinning FRR up and down reserves must be fulfilled by the generators in the system. See the [reserve tutorial](reserves) for more information about reserves and the [](reserve_group) object.

```{code-cell} ipython3
#Necessary imports used in all examples
import pandas as pd
from pyshop import ShopSession
import plotly.graph_objects as go
pd.options.plotting.backend = "plotly"

from sim_two_rsv import build_model
```

```{code-cell} ipython3
#Create a standard ShopSession
shop=ShopSession()
#Build a simple model with two reservoirs, two plants, and 6 generators.
build_model(shop)

#Add a reserve_group object and 15 MW of FRR up and down requirements
frr = shop.model.reserve_group.add_object("frr_group")
frr.frr_up_obligation.set(15)
frr.frr_down_obligation.set(15)

#Connect all generators to the reserve group
for gen in shop.model.generator:
    gen.connect_to(frr)

#Display topology to the screen
display(shop.model.build_connection_tree())
```

To get a baseline production plan for each generator, a basic optimization of the system is performed:

```{code-cell} ipython3
#Run an optimizatio
shop.start_sim([],5)    
shop.set_code("incremental",[])
shop.start_sim([],5)    
```

The optimized [production](generator:production) plan for each generator has now been set by the optimization, and the [start shopsim command](start_shopsim) can be used with the option "gen_mw_result". This will perform a simulation of the system based on the initial reservoir state, reservoir inflow, and the generator production plans. 

The output time resolution of the simulated timeseries attributes is set to 60 seconds with the optional value given to the shopsim command.

```{code-cell} ipython3
# Run a simulation based on the optimized production levels for each generator
shop.start_shopsim(["gen_mw_result"], 60)
```

After running the simulation, it is possible to retrieve various simulation result attributes on the regular SHOP objects. The simulated generator production values are stored as the attribute [](generator:sim_production), which are identical to the optimized production plans since the "gen_mw_result" option was used. This is confirmed in the figure below, which shows the sum plant production in the simulator and optimization:

```{code-cell} ipython3
for plant in shop.model.plant:
    fig = go.Figure()
    name = plant.get_name()
    sim_dis = plant.sim_production.get()
    opt_dis = plant.production.get()
    opt_dis[sim_dis.index[-1]] = opt_dis.values[-1]
    opt_dis = opt_dis.resample("1min").ffill()
    
    fig.add_trace(go.Scatter(x=opt_dis.index, y=opt_dis.values, name="optimized"))
    fig.add_trace(go.Scatter(x=sim_dis.index, y=sim_dis.values, name="simulated", line=dict(dash='dash')))

    fig.update_layout(xaxis_title="Time", yaxis_title="Power [MW]",title=f"Production: {name}")
    fig.show()
```

However, the simulated and optimized *discharge* is not identical. The [](plant:sim_discharge) attribute is systematically higher than the optimized discharge for Plant2, end exhibits a sawtooth pattern over each clock hour. This deviation comes from the fact that the reservoir levels change continuously inside each hourly step in the optimization. The PQ curves used in the optimization assume that the generator head is constant over the time step, while the simulator updates the reservoir levels much more rapidly (every 20 seconds in this case, see the [set simtimestep command](set_simtimestep)).

```{code-cell} ipython3
for plant in shop.model.plant:
    fig = go.Figure()
    name = plant.get_name()
    sim_dis = plant.sim_discharge.get()
    opt_dis = plant.discharge.get()
    opt_dis[sim_dis.index[-1]] = opt_dis.values[-1]
    opt_dis = opt_dis.resample("1min").ffill()
    
    fig.add_trace(go.Scatter(x=opt_dis.index, y=opt_dis.values, name="optimized"))
    fig.add_trace(go.Scatter(x=sim_dis.index, y=sim_dis.values, name="simulated", line=dict(dash='dash')))

    fig.update_layout(xaxis_title="Time", yaxis_title="Discharge [m3/s]",title=f"Discharge: {name}")
    fig.show()
```

The deviation in discharge are also be present in the simulated reservoir levels. The [](reservoir:sim_head) attribute is plotted together with the original reservoir trajectories from the optimization in the following figure:

```{code-cell} ipython3
for rsv in shop.model.reservoir:
    fig = go.Figure()
    name = rsv.get_name()
    sim_level = rsv.sim_head.get()
    opt_level = rsv.head.get()
    fig.add_trace(go.Scatter(x=opt_level.index, y=opt_level.values, name="optimized"))
    fig.add_trace(go.Scatter(x=sim_level.index, y=sim_level.values, name="simulated", line=dict(dash='dash')))

    fig.update_layout(xaxis_title="Time", yaxis_title="Level [masl]",title=f"Reservoir levels: {name}")
    fig.show()
```

The bottom reservoir is completely emptied at the end of the horizon in the optimization, but the simulation shows that the reservoir would actually empty a few minutes earlier due to a slightly higher water consumption to keep the production plan from the optimization. The attribute [](reservoir:sim_water_deficit) shows how much water has to be created to avoid going below [](reservoir:lrl) in the simulation, and is shown in the figure below:

```{code-cell} ipython3
fig = go.Figure()
for rsv in shop.model.reservoir:
    name = rsv.get_name()
    sim_vol_deficit = rsv.sim_water_deficit.get()
    fig.add_trace(go.Scatter(x=sim_vol_deficit.index, y=sim_vol_deficit.values, name=name))

fig.update_layout(xaxis_title="Time", yaxis_title="Volume [Mm3]",title=f"Created water during simulation")
fig.show()
```

## Simulating the activation of reserve capacity
The original optimization model included a reserve capacity requirement of 15 MW FRR up and down in every hour. What would happen to the system if all of the FRR up reserve capacity on the upstream plant was activated in all hours? To simulate this scenario, a [](generator:production_schedule) has to be set manually for all generators in the system. The option "gen_mw_schedule" must be used in this new simulation as the "gen_mw_result" option would ignore the production schedule attributes:

```{code-cell} ipython3
for gen in shop.model.generator:
    
    name = gen.get_name()
    prod = gen.production.get()
    
    # Set the production plan to the optimized production and 100% of the allocated FRR up capacity on Plant1 generators
    if "Plant1" in name:
        frr_up = gen.frr_up_delivery.get()
        gen.production_schedule.set(prod + frr_up)
    else:        
        gen.production_schedule.set(prod)
    
# Simulate the system once more with the new production schedules as input
shop.start_shopsim(["gen_mw_schedule"], 60)
```

As the generators of Plant1 have to fully activate the optimized FRR up capacity allocation, the simulated discharge is now considerably higher compared to the original discharge plan. The simulated discharge on Plant2 is still not identical to the optimized results, but the added water from Plant1 actually makes the deviation smaller compared to the first simulation.

```{code-cell} ipython3
for plant in shop.model.plant:
    fig = go.Figure()
    name = plant.get_name()
    sim_dis = plant.sim_discharge.get()
    opt_dis = plant.discharge.get()
    opt_dis[sim_dis.index[-1]] = opt_dis.values[-1]
    opt_dis = opt_dis.resample("1min").ffill()
    
    fig.add_trace(go.Scatter(x=opt_dis.index, y=opt_dis.values, name="optimized"))
    fig.add_trace(go.Scatter(x=sim_dis.index, y=sim_dis.values, name="simulated with activation", line=dict(dash='dash')))

    fig.update_layout(xaxis_title="Time", yaxis_title="Discharge [m3/s]",title=f"Discharge: {name}")
    fig.show()
```

The simulated reservoir level of the upper reservoir has a fairly large deviation from the optimized reservoir levels, but is still above the minimum level. Since more water is discharged from the upstream reservoir in the simulation, Reservoir2 now ends up slightly above the minimum reservoir level at the end of the horizon. This means that the activation of reserve capacity on Plant1 alleviated the water deficit seen in the original production plan.

```{code-cell} ipython3
for rsv in shop.model.reservoir:
    fig = go.Figure()
    name = rsv.get_name()
    sim_level = rsv.sim_head.get()
    opt_level = rsv.head.get()
    fig.add_trace(go.Scatter(x=opt_level.index, y=opt_level.values, name="optimized"))
    fig.add_trace(go.Scatter(x=sim_level.index, y=sim_level.values, name="simulated with activation", line=dict(dash='dash')))

    fig.update_layout(xaxis_title="Time", yaxis_title="Level [masl]",title=f"Reservoir levels: {name}")
    fig.show()
```