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 show 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 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 axis 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 are 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 bidning 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()