---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.14.6
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

+++ {"Collapsed": "false"}

(best-profit-basic-example)=
# Marginal costs and optimal uploading

+++ {"Collapsed": "false"}

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

- pyshop
    - [](bp-py)
    
# 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](best-profit).

```{code-cell} ipython3
:Collapsed: 'false'

#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](generator) and a second plant with one large and three small generators is built and optimized in the regular way:

```{code-cell} ipython3
:Collapsed: 'false'

#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](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.

```{code-cell} ipython3
:Collapsed: 'false'

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](plant:best_profit_mc) attribute. The attribute is a [XYT](datatype: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](plant: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](unit_combination:marginal_cost) and [average cost (AC) curves](unit_combination:average_cost) 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](best-profit).

```{code-cell} ipython3
:Collapsed: 'false'

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](generator:best_profit_p) and [best_profit_q](generator: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.

```{code-cell} ipython3
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:

```{code-cell} ipython3
:Collapsed: 'false'

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()
```

+++ {"Collapsed": "false"}

# Files

+++ {"Collapsed": "false"}

(bp-py)=
## bp.py

```{code-cell} ipython3
:Collapsed: 'false'
:tags: [remove-input]

with open('bp.py', 'r') as f:
    print(f.read())
```

```{code-cell} ipython3

```