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
System topology and initial simulation#
First, a simple SHOP model with two reservoirs and two plants 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 for more information about reserves and the reserve_group object.
#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
#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:
#Run an optimizatio
shop.start_sim([],5)
shop.set_code("incremental",[])
shop.start_sim([],5)
The optimized production plan for each generator has now been set by the optimization, and the start shopsim command 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.
# 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 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:
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 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).
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 sim_head attribute is plotted together with the original reservoir trajectories from the optimization in the following figure:
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 sim_water_deficit shows how much water has to be created to avoid going below lrl in the simulation, and is shown in the figure below:
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 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:
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.
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.
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()