Marginal costs and optimal uploading#
The model setup for this example is available in the following formats:
pyshop
Introduction#
The best profit (BP) functionality is the modern marginal cost and optimal plant uploading functionality in SHOP. The theory is described in more detail in the BP tutorial.
#Necessary imports used in all examples
import pandas as pd
from pyshop import ShopSession
import plotly.graph_objs as go
pd.options.plotting.backend = "plotly"
#Functions used in this example for building a tunnel model, adding a gate to a tunnel and running the optimization
from bp import build_model, run_model
Run a standard optimization#
First, a SHOP model with one plant with two equal generators and a second plant with one large and three small generators is built and optimized in the regular way:
#Create a standard ShopSession
shop=ShopSession()
#Build a simple model
build_model(shop)
#Optimize model by calling "run_model" in tunnel_model.py
run_model(shop)
#Display topology to the screen
display(shop.model.build_connection_tree())
#Display results for generator production
pd.DataFrame([obj.production.get().rename(obj.get_name()) for obj in shop.model.generator]).transpose().plot.bar(title="Generator production")
Calculate best profit curves#
After the optimization is done, the create bp_curves command is used to create marginal cost, average cost, and best profit curves for both plants. The calculation is done for the 72 first hours, specified in the call to the command. The option “production” is given to signal that the production results from the optimization should be used as the reference production when calculating the BP curve.
shop.create_bp_curves(["production"],["0","71"])
The BP curves for each time step has now been calculated and saved to the best_profit_mc attribute. The attribute is a XYT that holds the best profit curve for each time step where the x-values are the production levels and the y-values are the prices. The ref_prod TXY shows the reference production for the plant in each time step.
Several unit_combination objects have been automatically created to hold the marginal cost (MC) curves and average cost (AC) curves of each generator combination that has been considered in the best profit calculation. Note that this includes the empty combination of no units running. The unit_combination objects are automatically connected to the plant object, and can be iterated over in pyshop as shown in the code below.
The BP curves for each plant in all hours are shown in the figures below, in addition to the MC and AC curves of each unit combination. Use the slider to change which hour to display. The BP curves follow the MC curves for different unit_combinations and jumps between them given certain criteria of always having a positive regulation profit, see the BP tutorial.
from scipy.interpolate import interp1d
import numpy as np
n_hours = 72
#Create one plot for each plant
for plant in shop.model.plant:
#Get the BP curve and reference production of the plant
bp_curves = plant.best_profit_mc.get()
ref_prod = plant.ref_prod.get()
comb_mc=[]
comb_ac=[]
#Get the marginal and acerage cost curves of all considered unit combinations
for comb in plant.unit_combinations:
comb_mc.append(comb.marginal_cost.get())
comb_ac.append(comb.average_cost.get())
fig = go.Figure()
fig.update_layout(title="BP curve for " + plant.get_name() + " at " + str(bp_curves[0].name),
xaxis_title="MW",yaxis_title="€/MWh")
#Plot the BP, MC, and AC curves for each hour, in addition to the reference production point
for h in range(n_hours):
bp = bp_curves[h]
fig.add_trace(go.Scatter(x=bp.index, y=bp.values, name='BP curve',mode='lines+markers', line=dict(width=8)))
#Calculate the price of the reference production point on the BP curve by interpolating
ref_price = interp1d(bp.index, bp.values, bounds_error=False, fill_value=(bp.values[0], bp.values[-1]))(ref_prod.iloc[h])
fig.add_trace(go.Scatter(x=[ref_prod.iloc[h]], y=[ref_price], name='Ref prod', mode='markers',marker=dict(color='rgba(0,0,0, 0.0)', size=20, line=dict(color='Red', width=4))))
for comb_no, comb in enumerate(plant.unit_combinations):
mc = comb_mc[comb_no][h]
ac = comb_ac[comb_no][h]
fig.add_trace(go.Scatter(x=mc.index, y=mc.values, name=comb.get_name().replace(plant.get_name(),'MC'), mode='lines+markers', line=dict(width=4)))
fig.add_trace(go.Scatter(x=ac.index, y=ac.values, name=comb.get_name().replace(plant.get_name(),'AC'), mode='lines+markers', line=dict(width=4)))
n_traces_per_timestep=2+2*len(plant.unit_combinations)
#Only show curves for hour zero on the first figure displayed
for n in range(n_traces_per_timestep, len(fig.data)):
fig.data[n].visible = False
# Create and add slider to choose which hour to show
steps = []
for h in range(n_hours):
step = dict(
method="update",
label = str(bp_curves[h].name),
args=[{"visible": [False] * n_hours*n_traces_per_timestep},
{"title": "BP curve for " + plant.get_name() + " at " + str(bp_curves[h].name)}], # layout attribute
)
#Show only curves for hour h
for n in range(n_traces_per_timestep):
step["args"][0]["visible"][h*n_traces_per_timestep+n] = True
steps.append(step)
sliders = [dict(
active=0,
currentvalue={"prefix": "Time: "},
pad={"t": 50},
steps=steps
)]
fig.update_layout(sliders=sliders)
fig.show()
Optimal uploading#
The attributes best_profit_p and best_profit_q on the generator object shows the optimal uploading of each generator in the BP calculations. Only the plant production levels on the final BP curve is included.
The figures below show the optimal uploading of both plants for the first hour of the BP calculation. Plant1 has identical generators, and is uploaded by first starting G1, then starting G2 and uploading the generators in tandem. Note that the production level of G1 is decreased when G2 is first started. On Plant2, we can see that it is best to regulate up the large generator G1 and turning on the identical G2, G3, and G4 as the plant production increases. The smaller generators are held on more or less the same production levels while the production on G1 is increased.
for p in shop.model.plant:
fig = go.Figure()
fig.update_layout(title=p.get_name()+" optimal uploading",xaxis_title="Plant production [MW]",yaxis_title="Generator production [MW]")
for g in p.generators:
gen_p = g.best_profit_p.get()
fig.add_trace(go.Scatter(x=gen_p[0].index, y=gen_p[0].values, name=g.get_name(),mode='lines+markers'))
fig.show()
Finally, the reference plant production and price from the BP curve in each hour is plotted in 3D:
for plant in shop.model.plant:
bp_curves = plant.best_profit_mc.get()
ref_prod = plant.ref_prod.get()
ref_price=pd.Series(index=ref_prod.index, dtype='float64')
for h in range(n_hours):
bp = bp_curves[h]
ref_price.values[h] = interp1d(bp.index, bp.values, bounds_error=False, fill_value=(bp.values[0], bp.values[-1]))(ref_prod.iloc[h])
df = pd.DataFrame()
fig = go.Figure()
fig.add_trace(go.Surface(x=df.columns.tolist(), y=df.index.tolist(), z=df.values.tolist(), colorscale="Viridis", hidesurface=True, contours_x_show=True, contours_y_show=True))
fig.update_traces(showscale=False)
fig.add_trace(go.Scatter3d(x=ref_prod.index, y=ref_prod.values, z=ref_price.values, mode='lines+markers', line=dict(width=4), marker=dict(size=2)))
fig.update_layout(
autosize=True,
title="Ref production and BP price for every hour for "+plant.get_name()
)
fig.update_scenes(
aspectratio=dict(x=1.5, y=1, z=1),
xaxis_title='Time',
yaxis_title='Production [MW]',
zaxis_title='Price [€/MWh]'
)
fig.show()
Files#
bp.py#
import pandas as pd
def build_model(shop):
starttime = pd.Timestamp('2018-01-23 00:00:00')
endtime = pd.Timestamp('2018-01-26')
shop.set_time_resolution(starttime=starttime, endtime=endtime, timeunit="hour", timeresolution=pd.Series(index=[starttime],data=[1]))
rsv1 = shop.model.reservoir.add_object('Reservoir1')
rsv1.max_vol.set(39)
rsv1.lrl.set(860)
rsv1.hrl.set(905)
rsv1.vol_head.set(pd.Series([860, 906, 907], index=[0, 39, 41.66], name=0))
rsv2 = shop.model.reservoir.add_object('Reservoir2')
rsv2.max_vol.set(97.5)
rsv2.lrl.set(650)
rsv2.hrl.set(679)
rsv2.vol_head.set(pd.Series([650, 679, 680], index=[0, 97.5, 104.15], name=0))
plant1 = shop.model.plant.add_object('Plant1')
plant1.outlet_line.set(672)
plant1.main_loss.set([0])
plant1.penstock_loss.set([0.001])
plant1.mip_flag.set(1)
for gen_no in range(2):
gen=shop.model.generator.add_object(f"{plant1.get_name()}_G{str(gen_no+1)}")
gen.connect_to(plant1)
gen.penstock.set(1)
gen.p_min.set(60)
gen.p_max.set(120)
gen.p_nom.set(120)
gen.startcost.set(300)
gen.gen_eff_curve.set(pd.Series([100, 100], index=[60, 120]))
gen.turb_eff_curves.set([pd.Series([86.2479, 87.5762, 88.7758, 89.8467, 90.7888, 91.6022, 92.2869, 92.8429, 93.2701, 93.5686, 93.7383, 93.7793, 93.6916, 93.4751, 93.1299, 92.6560, 92.0534],
index=[28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60],
name=170)])
plant2 = shop.model.plant.add_object('Plant2')
plant2.outlet_line.set(586)
plant2.main_loss.set([0])
plant2.penstock_loss.set([0.0001,0.0002])
plant2.mip_flag.set(1)
for gen_no in range(4):
gen=shop.model.generator.add_object(f"{plant2.get_name()}_G{str(gen_no+1)}")
gen.connect_to(plant2)
if gen_no == 0:
gen.penstock.set(1)
gen.p_min.set(100)
gen.p_max.set(180)
gen.p_nom.set(180)
gen.startcost.set(300)
gen.gen_eff_curve.set(pd.Series([100, 100], index=[100, 180]))
gen.turb_eff_curves.set([pd.Series([92.7201, 93.2583, 93.7305, 94.1368, 94.4785, 94.7525, 94.9606, 95.1028, 95.1790, 95.1892, 95.1335, 95.0118, 94.8232, 94.5191],
index=[126.54, 137.03, 147.51, 158.00, 168.53, 179.01, 189.50, 199.98, 210.47, 220.95, 231.44, 241.92, 252.45, 264.74],
name=60)])
else:
gen.penstock.set(2)
gen.p_min.set(30)
gen.p_max.set(55)
gen.p_nom.set(55)
gen.startcost.set(300)
gen.gen_eff_curve.set(pd.Series([100, 100], index=[30, 55]))
gen.turb_eff_curves.set([pd.Series([83.8700, 85.1937, 86.3825, 87.4362, 88.3587, 89.1419, 89.7901, 90.3033, 90.6815, 90.9248, 91.0331, 91.0063, 90.8436, 90.4817],
index=[40.82, 44.20, 47.58, 50.97, 54.36, 57.75, 61.13, 64.51, 67.89, 71.27, 74.66, 78.04, 81.44, 85.40],
name=60)])
rsv1.connect_to(plant1)
plant1.connect_to(rsv2)
rsv2.connect_to(plant2)
rsv1.start_head.set(900)
rsv2.start_head.set(670)
rsv1.energy_value_input.set(30)
rsv2.energy_value_input.set(10)
rsv2.inflow.set(pd.Series([60], [starttime]))
da = shop.model.market.add_object('Day_ahead')
da.sale_price.set(pd.DataFrame([32.992,31.122,29.312,28.072,30.012,33.362,42.682,74.822,77.732,62.332,55.892,46.962,42.582,40.942,39.212,39.142,41.672,46.922,37.102,32.992,31.272,29.752,28.782,28.082,27.242,26.622,25.732,25.392,25.992,27.402,28.942,32.182,33.082,32.342,30.912,30.162,30.062,29.562,29.462,29.512,29.672,30.072,29.552,28.862,28.412,28.072,27.162,25.502,26.192,25.222,24.052,23.892,23.682,26.092,28.202,30.902,31.572,31.462,31.172,30.912,30.572,30.602,30.632,31.062,32.082,36.262,34.472,32.182,31.492,30.732,29.712,28.982],
index=[starttime + pd.Timedelta(hours=i) for i in range(0,72)]))
da.max_sale.set(pd.Series([9999], [starttime]))
da.buy_price.set(da.sale_price.get()+0.002)
da.max_buy.set(pd.Series([9999], [starttime]))
def run_model(shop):
shop.start_sim([], ['5'])
shop.set_code(['incremental'], [])
shop.start_sim([], ['3'])