Cut descriptions

The model setup for the two cut description examples are available in the following formats.

• pyshop

  • cuts.py

  • rhs_file.txt

  • wv_file.txt

• YAML

  • model.yaml

  • standard_cuts.yaml

  • price_dep_cuts.yaml

• ASCII

  • model.ascii

  • name_coupling.txt

  • cut_order.txt

  • standard_cuts.txt

  • price_dep_cuts.txt

Standard cut description

Adding a standard cut description to a simple two-reservoir system is shown in this example. The cuts used in the example were generated by the prototype cut creation functionality in SHOP, and is therefore not saved to the standard text file format used by ProdRisk and EOPS.

#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 and solving a simple model with cuts
from cuts import build_model, run_model, read_cuts, get_minimal_cut_surface

#Create a standard ShopSession
shop=ShopSession()
#Build a simple model with two reservoirs and two plants.
build_model(shop)
#Display topology to the screen
display(shop.model.build_connection_tree())

After the topology has been created, we read in the cuts from the text files. A new cut_group object must be created to hold the right-hand side constants of the cut constraints in the rhs attribute. The water_value_input attribute on the reservoir objects are used to store the cut coefficients with corresponding reference volumes. All reservoirs that are part of the cut constraints must be connected to the cut_group.

#Read the cut info
rhs_in,cut_coeffs = read_cuts()
n_cuts = len(rhs_in)

#Create a cut group object and add the RHS
my_cuts = shop.model.cut_group.add_object("my_cuts")
my_cuts.rhs.set([pd.Series(rhs_in, index=[0]*n_cuts, name=0)])

#Add the cut coefficients to the reservoirs and connect them to the cut_group
for rsv in shop.model.reservoir:
    name = rsv.get_name()
    
    v_ref = cut_coeffs[name][0]
    wv = cut_coeffs[name][1]
    
    rsv.water_value_input.set([pd.Series(wv, index=v_ref, name=0)])
    rsv.connect_to(my_cuts)    

The inequality constraints defined by each cut \(k\) is

\(c \leq \sum_i \textup{WV}_{ik}\cdot(v_i - V_i^{\textup{ref}}) + \textup{rhs}_k \qquad \text{for all } k,\)

when formulated as profit maximization. The single cut variable \(c\) is bounded from above by all the cuts, and the optimization problem tries to maximize its value. The cut coefficients \(\textup{WV}_{ik}\) for reservoir \(i\) is multiplied with the end volume variable \(v_i\) subtracted by the cut reference volume. When all \(v_i = V_i^{\textup{ref}}\), the cut is bounded by the right-hand side value \(\textup{rhs}_k\).

The cut constraints can be visualized as a 3D plot since there are only two reservoirs in this example. The x- and y-axes correspond to the end volume of each reservoir, while the z-axis shows the cut value in terms of future expected income. Each cut is a plane in 3D which is spanned around a reference point with the reference volumes as x and y coordinates and the rhs as the z-value. These points are plotted in 3D scatter plot below. The minimal surface spanned by all these cuts is also plotted, which is the effective bound that the optimization problem will see. The surface is found by taking the lowest (most restrictive) cut value over the set of cuts for different values of the reservoir volumes.

Note that the reference points in red are marked as “never binding” and hover above the calculated minimal surface. This should not happen, as a cut should at least be binding exactly at the reference point. Non-convexities in the optimization program used to create these sample cuts are to blame in this case, as the SHOP cut creation functionality is still experimental. Some of the binding cuts are also not binding at the reference point but somewhere else in the domain, which is also inaccurate.

import numpy as np
rsv1 = shop.model.reservoir.Reservoir1
rsv2 = shop.model.reservoir.Reservoir2

vmax1 = rsv1.max_vol.get()
vmax2 = rsv2.max_vol.get()

v_ref1 = np.array(cut_coeffs["Reservoir1"][0])
v_ref2 = np.array(cut_coeffs["Reservoir2"][0])

v1 = np.arange(0,vmax1,0.25)
v2 = np.arange(0,vmax2,0.25)

#Add the reference volumes to the volume ranges to have a good representation for the surface around the reference points
v1 = np.concatenate((v1,v_ref1))
v1.sort()
v2 = np.concatenate((v2,v_ref2))
v2.sort() 

#Calculate the minimal cut surface and the binding cuts
binding_cuts, minimal_surface = get_minimal_cut_surface(rhs_in,cut_coeffs,v1,v2)

#Plot the cut surface
fig = go.Figure()
fig.add_trace(go.Surface(z=minimal_surface, x=v1*100/vmax1, y=v2*100/vmax2,hoverinfo="skip"))

colors = ["red"]*n_cuts
text = [f"Cut {i+1}, never binding" for i in range(n_cuts)]
for k in binding_cuts:
    colors[k] = "black"
    text[k] = f"Cut {k+1}, binding"
    
#Plot the reference points
fig.add_trace(go.Scatter3d(z=rhs_in, x=v_ref1*100/vmax1, y=v_ref2*100/vmax2,mode='markers',marker_color=colors,hovertext=text))

fig.update_scenes(
    aspectratio=dict(x=1, y=1, z=1),
    xaxis_title='Rsv1 vol [%]',
    yaxis_title='Rsv2 vol [%]',
    zaxis_title='Cut value'
)
fig.update_layout(
    width=1000,
    height=900,
    autosize=True)
    
fig.show()

Now we can run the optimization and look at the results. The final value in the TimeSeries attribute end_value on the cut_group object holds the optimized cut value (future expected (negative) income). The binding_cut_up attribute holds the cut constraint that was the most constraining cut in this optimization.

#Optimize model by calling "run_model"
run_model(shop)

#Get the final end value of the cut group, which is the future expected (negative) profit
cut_val = -my_cuts.end_value.get().iloc[-1]

#Get the binding cut constraint and print the cut information
cut_index = my_cuts.binding_cut_up.get().iloc[-1]
cut_index = int(cut_index)

rhs_out = my_cuts.rhs.get()[0]
active_rhs = rhs_out.values[cut_index]

print(f"The future expeced income is {cut_val:.2f} € and the binding cut is number {cut_index+1} out of {n_cuts} with the following values:")
print(f"RHS is {active_rhs:.2f} €")
for rsv in shop.model.reservoir:
    name = rsv.get_name()
    
    cut_coeffs = rsv.water_value_input.get()[0]
    v_ref = cut_coeffs.index[cut_index]
    wv = cut_coeffs.values[cut_index]

    end_vol = rsv.storage.get().iloc[-1]

    print(f"Cut coefficient {wv:.2f} €/Mm3 and reference volume {v_ref:.2f} Mm3 for {name} with end volume {end_vol:.2f} Mm3")

The future expeced income is 1743877.53 € and the binding cut is number 8 out of 9 with the following values:
RHS is 1591330.63 €
Cut coefficient 20069.49 €/Mm3 and reference volume 33.33 Mm3 for Reservoir1 with end volume 35.50 Mm3
Cut coefficient 5270.82 €/Mm3 and reference volume 6.81 Mm3 for Reservoir2 with end volume 27.52 Mm3

Note that the attribute binding_cut_down is not used when the cuts are not price dependent, and have the default value of:

i = my_cuts.binding_cut_down.get().iloc[-1]
print(f"'binding_cut_down' attribute for standard cuts: {i}")

'binding_cut_down' attribute for standard cuts: -1.0

The evolution of the reservoir storage and the global water value result time series’ are shown in the figures below. The final value of the global water value time series is equal to the binding cut coefficient for each reservoir.

pd.DataFrame([rsv.storage.get().rename(rsv.get_name()) for rsv in shop.model.reservoir]).transpose().plot(title="Reservoir storage")
pd.DataFrame([-rsv.water_value_global_result.get().rename(rsv.get_name()) for rsv in shop.model.reservoir]).transpose().plot(title="Reservoir global water value")

Price-dependent cut description

To illustrate how price-dependent cuts can be added to SHOP, we will use the same cuts as in the example above and synthesize some price variations.

import numpy as np

#Create a standard ShopSession
shop=ShopSession()
#Build a simple model with two reservoirs and two plants.
build_model(shop)

#Read the cut info
rhs_in,cut_coeffs = read_cuts()
n_cuts = len(rhs_in)

#Create a cut group object
my_cuts = shop.model.cut_group.add_object("my_cuts")
for rsv in shop.model.reservoir:
    rsv.connect_to(my_cuts)

We use the average price in the SHOP run and create four arbitrary price levels around it. Only the two closest price levels to the average price will actually be used in SHOP, so the highest and lowest price levels are not necessary for this particular model run. The rhs and cut coefficients are then scaled linearly for each price level.

Since SHOP also uses the average price to interpolate the cuts in the two closest price levels, the solution we get from using the scaled price-dependent cuts will be identical to the solution in the previous example with standard cuts. This is assuming that the models are set up in the exact same way and that a MIP gap of zero is reached in each iteration in both models. To get some different results in this example, we add some “random” noise to the scaling in the code below. Try to comment out the code adding the noise and see for yourself if the results are the same as in the original case!

#Find the average price and make some new price levels around it
market = shop.model.market.Day_ahead
average_price = np.mean(market.sale_price.get())

prices = [average_price*0.5,average_price*0.8,average_price*1.12,average_price*1.3]

scaling = [p/average_price for p in prices]

#Comment out the two lines below to remove the noise applied to the cut scaling
noise = [0.95,0.96,1.04,1.05]
scaling = [s*n for s,n in zip(scaling,noise)]

#Set the scaled RHS, the price is set as a reference value
my_cuts.rhs.set([pd.Series(np.array(rhs_in)*scale, index=[0]*n_cuts, name=price) for scale,price in zip(scaling,prices)])

#Add the scaled cut coefficients, the price is set as a reference value
for rsv in shop.model.reservoir:
    name = rsv.get_name()
    
    v_ref = cut_coeffs[name][0]
    wv = cut_coeffs[name][1]
    
    rsv.water_value_input.set([pd.Series(np.array(wv)*scale, index=v_ref, name=price) for scale,price in zip(scaling,prices)])

The scatter plots below shows the scaled and original cuts in terms of the rhs and cut coefficients. The values are plotted in 2D for simplicity:

#Plot the scaled and original RHS values

rhs = my_cuts.rhs.get()

fig = go.Figure(layout={'title':f"RHS",'xaxis_title':"Cut number",'yaxis_title':"RHS [€]"})
for r in rhs:
    fig.add_trace(go.Scatter(name=f"Price {r.name:.2f}",x=list(range(1,n_cuts+1)),y=r.values,mode='markers'))    

fig.add_trace(go.Scatter(name=f"Original",x=list(range(1,n_cuts+1)),y=rhs_in,mode='markers'))    
fig.show()

#Plot the scaled and original (ref volume, cut coefficient) pairs
for rsv in shop.model.reservoir:    
    wv_in = rsv.water_value_input.get()
    name = rsv.get_name()
    
    fig = go.Figure(layout={'title':f"{name}",'xaxis_title':"Ref volume [Mm3]",'yaxis_title':r"Cut coefficient [€/Mm3]"})
    for wv in wv_in:
        fig.add_trace(go.Scatter(name=f"Price {wv.name:.2f}",x=wv.index,y=wv.values,text=[f"Cut number {i+1}" for i in range(n_cuts)],mode='markers'))    
    
    fig.add_trace(go.Scatter(name=f"Original",x=cut_coeffs[name][0],y=cut_coeffs[name][1],text=[f"Cut number {i+1}" for i in range(n_cuts)],mode='markers'))    
    
    fig.show()

Now we run the model and take a look at the final cut value and the binding cuts:

#Optimize model by calling "run_model"
run_model(shop)

#Get the final end value of the cut group, which is the future expected (negative) profit
cut_val = -my_cuts.end_value.get().iloc[-1]
print(f"The future expeced income is {cut_val:.2f} €\n")

#Get the binding cut constraints in the two price levels and print the cut information
cut_index_up = my_cuts.binding_cut_up.get().iloc[-1]
cut_index_up = int(cut_index_up)

cut_index_down = my_cuts.binding_cut_down.get().iloc[-1]
cut_index_down = int(cut_index_down)

rhs_out = my_cuts.rhs.get()
active_rhs_up = rhs_out[2].values[cut_index_up]
active_rhs_down = rhs_out[1].values[cut_index_down]

print(f"The binding cut in the set of cuts over the SHOP price is number {cut_index_up+1} out of {n_cuts} with the following values:")
print(f"\tRHS is {active_rhs_up:.2f} €")
for rsv in shop.model.reservoir:
    name = rsv.get_name()
    
    cut_coeffs = rsv.water_value_input.get()[2]
    v_ref = cut_coeffs.index[cut_index_up]
    wv = cut_coeffs.values[cut_index_up]

    end_vol = rsv.storage.get().iloc[-1]

    print(f"\tCut coefficient {wv:.2f} €/Mm3 and reference volume {v_ref:.2f} Mm3 for {name} with end volume {end_vol:.2f} Mm3")
print("")

print(f"The binding cut in the set of cuts below the SHOP price is number {cut_index_down+1} out of {n_cuts} with the following values:")
print(f"\tRHS is {active_rhs_down:.2f} €")
for rsv in shop.model.reservoir:
    name = rsv.get_name()
    
    cut_coeffs = rsv.water_value_input.get()[1]
    v_ref = cut_coeffs.index[cut_index_down]
    wv = cut_coeffs.values[cut_index_down]

    end_vol = rsv.storage.get().iloc[-1]

    print(f"\tCut coefficient {wv:.2f} €/Mm3 and reference volume {v_ref:.2f} Mm3 for {name} with end volume {end_vol:.2f} Mm3")

The future expeced income is 1788304.21 €

The binding cut in the set of cuts over the SHOP price is number 8 out of 9 with the following values:
	RHS is 1853581.92 €
	Cut coefficient 23376.95 €/Mm3 and reference volume 33.33 Mm3 for Reservoir1 with end volume 36.20 Mm3
	Cut coefficient 6139.45 €/Mm3 and reference volume 6.81 Mm3 for Reservoir2 with end volume 27.92 Mm3

The binding cut in the set of cuts below the SHOP price is number 8 out of 9 with the following values:
	RHS is 1222141.93 €
	Cut coefficient 15413.37 €/Mm3 and reference volume 33.33 Mm3 for Reservoir1 with end volume 36.20 Mm3
	Cut coefficient 4047.99 €/Mm3 and reference volume 6.81 Mm3 for Reservoir2 with end volume 27.92 Mm3

The final value of the global water value result time series is equal to the weighted sum of the binding cut coefficients in the price levels above and below the average SHOP price. The interpolation weights used in SHOP are:

w_up = (average_price - prices[1])/(prices[2]-prices[1])
w_down = 1 - w_up
print(f"Weight for cuts with price {(prices[1]):.2f}: {w_down:.3f}")
print(f"Weight for cuts with price {(prices[2]):.2f}: {w_up:.3f}")

Weight for cuts with price 26.53: 0.375
Weight for cuts with price 37.14: 0.625

pd.DataFrame([rsv.storage.get().rename(rsv.get_name()) for rsv in shop.model.reservoir]).transpose().plot(title="Reservoir storage")
pd.DataFrame([-rsv.water_value_global_result.get().rename(rsv.get_name()) for rsv in shop.model.reservoir]).transpose().plot(title="Reservoir global water value")