Bidding under uncertainty

The model setup is available in the following formats:

• pyshop

  • multi_price_bid_matrix_model.py

In this example we show how to import a spreadsheet containing multiple prices (up to 52 time series) into a dataframe. This price input will be used to create the same amount of scenarios which we in turn feed to SHOP. The output of this multiple price input run from SHOP will result into a joint bid matrix which consider all the stochastic price inputs to SHOP with all its price data intact, ready to make use of in a marked bidding situation with uncertainty regarding the price forecasts.

Imports and settings

The first thing we do is to import the needed packages. You can import whichever packages you like, however we use the following ones for this example:

• Pandas for structuring our data into dataframes

• Pyshop in order to create a SHOP session

• Plotly as backend for dynamic graph plotting

import pandas as pd
from pyshop import ShopSession
import plotly.express as px
pd.options.plotting.backend = "plotly"

Additionally, we import basic SHOP functions and data from a predefined demo dataset, see this section of the documentation:

from multi_price_bid_matrix_model import build_model, run_model

Instancing SHOP and building the model

In order to have SHOP receive our inputs, run the model we create and give us results, we need to instance a running SHOP session. You may create multiple SHOP sessions simultaneously if needed.

#Create a standard ShopSession
shop = ShopSession()

We then build our model using the existing function imported from bp.py

#Build a simple model with one plant with two equal generators and a second plant with one large and three small generators 
#by calling function "build_model" in bp.py
build_model(shop)

The imported model can now be visualized:

#Display topology to the screen
display(shop.model.build_connection_tree())

Data preperation

Since we we need to work a bit with time horizons given that we will consider part of a period with stochastic information, we retrive the start time from the imported data

# Retrieve the start time for adding stochastic time series later
timeres=shop.get_time_resolution()
shop.get_time_resolution()

{'starttime': Timestamp('2018-01-23 00:00:00'),
 'endtime': Timestamp('2018-01-26 00:00:00'),
 'timeunit': 'hour',
 'timeresolution': 2018-01-23 00:00:00    1.0
 2018-01-23 01:00:00    1.0
 2018-01-23 02:00:00    1.0
 2018-01-23 03:00:00    1.0
 2018-01-23 04:00:00    1.0
                       ... 
 2018-01-25 20:00:00    1.0
 2018-01-25 21:00:00    1.0
 2018-01-25 22:00:00    1.0
 2018-01-25 23:00:00    1.0
 2018-01-26 00:00:00    1.0
 Name: data, Length: 73, dtype: float64}

Importing prices from spreadsheet

After the model data has been imported and read, we move on to the prices we want to consider. We first define the number of prices we want to import

# Import number of prices
n_prices = 10

which will be equal to the scenarios we will create later on. We then define for how long of a period we want to consider the stochastic price input. This can not be a larger period than what you have data for in your price forecasts.

# Set the hours to import stochastic price
stochastic_prices_start = timeres['starttime']
# In this example, we choose 24 hours
stochastic_prices_end = stochastic_prices_start + pd.Timedelta(hours=23)

Next we import the stochastic prices accordingly into a dataframe

# Use pandas read_excel function to import n_prices to a dataframe, starting from the first column to the left in the spreadsheet
stochastic_prices_from_file = pd.read_excel ('52_hourly_prices_random.xls', header=None,usecols=list(range(0,n_prices)))

# Optionally, select which scenarios you want relative to the columns in your spreadsheet. They must then be equal to n_prices;
#stoch_price_from_file = pd.read_excel (r'20200123_scens_edit.xlsx', usecols=[1,2,5,7,9])

stochastic_prices_from_file

	0	1	2	3	4	5	6	7	8	9
0	38.688925	30.841437	39.303183	43.122221	31.399708	38.211738	42.479485	29.838857	29.932279	46.548664
1	38.916484	28.478524	33.843385	36.856333	40.447140	34.473711	33.751629	33.699346	36.987757	38.796374
2	37.133973	29.923425	31.280585	40.117271	33.285538	34.728074	34.096047	32.085033	36.426590	45.003989
3	28.567145	25.738590	33.428984	39.937462	39.570224	31.558648	38.290842	37.527192	28.990410	40.885543
4	32.909629	30.642448	34.668521	47.708656	36.963273	33.749174	38.703468	31.802267	29.553703	46.706921
5	34.438556	34.623470	33.969606	41.136482	35.480291	41.834732	33.788576	28.953164	39.690668	46.564680
6	35.470862	30.479294	35.031910	41.902961	37.544836	33.300794	41.630947	42.224127	40.853218	38.904588
7	43.884363	51.591232	43.362872	48.165823	43.549434	48.468647	53.586396	43.276922	46.691875	51.359444
8	46.485712	49.473627	44.295022	54.174682	47.947510	53.740830	50.935949	54.031818	53.416069	58.500872
9	51.831571	48.735697	56.201060	50.879699	54.291574	54.693474	58.224800	53.485508	64.635484	54.939225
10	42.770787	41.493970	37.828462	42.263994	49.966896	54.106208	44.128320	49.391066	60.575104	45.339350
11	37.870600	37.425262	44.116926	41.932840	42.416888	46.175475	45.121995	49.683231	40.777990	45.914877
12	38.720126	42.771099	41.507051	42.146804	38.803972	48.316620	40.811780	47.991236	53.715754	41.338486
13	36.119370	43.060226	38.300180	35.631930	39.999911	46.296312	45.363284	32.378365	42.465181	34.729646
14	38.194049	46.108626	37.375087	36.575645	44.014637	48.480447	50.672425	40.079081	41.311692	49.711962
15	38.260227	51.175027	39.216375	38.053676	46.202766	47.608828	44.795150	39.291503	51.563778	37.664475
16	38.058048	47.077703	40.799196	45.250958	43.690508	44.294230	50.905649	46.362381	46.659761	44.224564
17	56.710367	47.710875	55.554464	49.066695	54.546479	52.044637	56.701394	53.911895	60.111022	53.380406
18	59.004344	56.797886	60.172981	58.237505	57.513100	55.195022	54.109290	73.961415	56.121528	55.876500
19	55.776245	58.245847	58.270617	54.601773	57.317745	56.954792	53.898015	70.852715	59.414145	53.165951
20	43.244210	50.687627	49.486841	49.714059	37.706558	46.243266	47.601519	58.956010	43.540229	49.330592
21	48.317103	41.461231	45.007857	43.138390	40.507536	44.108610	41.700816	51.102177	40.245707	39.758460
22	38.237187	47.593982	42.574141	36.780849	49.548371	38.097889	33.731909	51.826269	46.920595	42.529905
23	42.371404	35.496953	46.265499	38.247478	43.228214	36.816092	42.987433	46.249756	39.407231	41.273543

In order to later combine the stochastic prices with a deterministic price into the future, we make sure the dataframe has the correct timestamp indexing

# Add correct indexing to the dataframe
stochastic_prices_from_file.index=pd.date_range(stochastic_prices_start,stochastic_prices_end,freq='h')

Then we can plot and review the newly imported prices

# Plot all stochastic prices imported from Excel
stochastic_prices_from_file.plot(title="Stochastic prices imported from spreadsheet", labels=dict(index="Time", value="EUR/MWh", variable=""))

Creating scenarios

Now it is time to create scenarios that we can populate with the imported prices. We make sure to create just as many scenarios as we have prices imported.

# Generate as many scenarios as prices
n_scenarios = n_prices
for i in range (1, n_scenarios+1):
    scenario_name='S'+str(i)
    # The first scenario always exists in SHOP and should not be added again
    if i>1:
        scenario = shop.model.scenario.add_object(scenario_name)
    else:
        scenario=shop.model.scenario[scenario_name]
    scenario.scenario_id.set(i)
    
    # Set each scenario equally probable
    scenario.probability.set(1.0/n_scenarios)
    
    # Branch immediately, i.e. at 'starttime'
    scenario.common_scenario.set(pd.Series([i], index=[timeres['starttime']]))
    
    # Optionally set branching to start after given number of hours (all scenarios are set equal to scenario 2 before this) 
    #scen.common_scenario.set(pd.Series([2, i], index=[timeres['starttime'], branching_start_time]))

Creating the new price array from stochastic and deterministic prices

Since we have chosen to only consider the first 24 hours as stochastic when it comes to the price, but have longer total time horizon, we need to combine the stochastic and deterministic prices into a joint price dataframe. We have already defined the start and end time for the stochastic price, but need to make sure that we also define when the deterministic price should be valid and thus overlap each other.

# Use first (and only) market as index for setting stochastic data and getting results
name_list = shop.model.market.get_object_names()

# The deterministic model only has one market, so we give stochastic price to that (da is just a name we use to indicate that we regard it to be a day ahead energy market)
day_ahead_market = shop.model.market[name_list[0]]

# Get the deterministic price from ASCII import for correct indexing 
deterministic_price = day_ahead_market.sale_price.get()

# Create a new array for deterministic price with n_prices columns in order to combine it with the stochastic price
deterministic_price_multi_dimension_array = pd.DataFrame(index=deterministic_price.index, columns=list(range(0,n_prices)))

# Define when to use deterministic price
deterministic_price_start = timeres['starttime']+pd.Timedelta(hours=24)
deterministic_price_end = timeres['endtime']-pd.Timedelta(hours=1)

We print out the start and end times to make sure they are correct.

# Print out start and end set point for prices
print("Stochastic price(s) start time",stochastic_prices_start)
print("Stochastic price(s) end time",stochastic_prices_end)
print("Deterministic prices start time",deterministic_price_start)
print("Deterministic prices end time",deterministic_price_end)

Stochastic price(s) start time 2018-01-23 00:00:00
Stochastic price(s) end time 2018-01-23 23:00:00
Deterministic prices start time 2018-01-24 00:00:00
Deterministic prices end time 2018-01-25 23:00:00

Lastly, we need to combine the stochastic and deterministic price inputs. We do this by creating a new dataframe with the right dimensions and indexes, and then populate it according to the start and end times defined above.

# Creating a new price dataframe which will combine both stochastic and deterministic price
combined_price = pd.DataFrame(index=deterministic_price.index, columns=list(range(0,n_prices)))

# Loop over the new price array and assign stochastic and deterministic prices

for j in range(0,n_prices):
    deterministic_price_multi_dimension_array[j]=deterministic_price
for i in pd.date_range(stochastic_prices_start,stochastic_prices_end,freq='h'):
       combined_price.loc[i] = stochastic_prices_from_file.loc[i]
for i in pd.date_range(deterministic_price_start,(deterministic_price_end),freq='h'):
        combined_price.loc[i] = deterministic_price_multi_dimension_array.loc[i]

Then we can plot out the resulting combined price.

# Plot the combined price
combined_price.plot(title="Combined price input", labels=dict(index="Time", value="EUR/MWh", variable=""))

This price is then set as the day ahead market price to consider for all scenarios. We also define a slightly higher buy-back price

day_ahead_market.sale_price.set(combined_price)

# Make sure that the buy price is slightly higher, so no arbitrage occurs

day_ahead_market.buy_price.set(day_ahead_market.sale_price.get()+0.1)

Create bid groups and configure limits

Now we need to create a bid group, which is the connection between a set of hydropower plants and their bid matrix.

# Create a bid group
bg=shop.model.bid_group.add_object('bg')

# Add all plants in the system to the bid group
for plant in shop.model.plant:
    bg.connect_to(plant)

# Defining which periods the bid curve should consider
day_ahead_market.bid_flag.set(pd.Series([1,0],index=[stochastic_prices_start,stochastic_prices_end]))

Running multi price scenarios in SHOP

It is time to optimize. We call the predefined function run_model, which calls for five full and three incremental iterations. Since we have defined scenarios in the SHOP instance earlier, all scenarios will be computed and optimized with a single call.

#Optimize model by calling "run_model" bp.py
run_model(shop)

Creating and plotting the bid matrix

Once SHOP has finished optimizing all scenarios, we can retrive and start processing the bid matrix the way we want it.

# Get bid matrix as an array
bid_result=bg.bid_curves.get()

# Convert bid matrix to a dataframe structure of choice
bid_matrix=pd.DataFrame(index=bid_result[0].index)
for t in range (0, 23):
    bid_matrix[t]=bid_result[t].values

# Transpose the bid matrix if necessary depending on the viewing needs 
bid_matrix_transposed=bid_matrix.transpose()

When we are satisfied with a data structure that fulfill our needs, we can either print it or plot it to review it

# Plotting bid matrices with different properties
fig = px.imshow(bid_matrix.transpose(), color_continuous_scale='spectral')
fig.update_layout(title="Bid matrix heatmapped (and inverted compared to normal bid matrix)", xaxis_title ="Market price", yaxis_title="Hours")
fig.update_yaxes(autorange=True) 
fig.show()

# Plotting bid matrices with different properties
bid_matrix_transposed.plot.bar(title="Bid matrix per hour", labels=dict(index="Hour [h]", value="Market price", variable="")).update_layout(barmode='group')

Plotting other results

If we want to see how each scenario suggest to produce, we can of course have a look at that as well

# Get the optimal production value(s)
optimal_production=day_ahead_market.sale.get()

# Plot optimal production in different plot types
optimal_production.plot.bar(title="Bid volumes per stochastic price #", labels=dict(index="Time", value="Sale/Production [MW]", variable=""))

# Plot optimal production in different plot types
optimal_production.plot.bar(title="Bid volumes per stochastic price #, stacked", labels=dict(index="Time", value="Sale/Production [MW]", variable="")) # bargap=.1, 

Files

multi_price_bid_matrix_model.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([85.8733, 87.0319, 88.0879, 89.0544, 89.9446, 90.7717, 91.5488, 92.2643, 92.8213, 93.1090, 93.2170, 93.0390, 92.6570, 92.1746],
                                    index=[28.12, 30.45, 32.78, 35.11, 37.45, 39.78, 42.11, 44.44, 46.77, 49.10, 51.43, 53.76, 56.10, 58.83],
                                    name=170),
                          pd.Series([86.7321, 87.9022, 88.9688, 89.9450, 90.8441, 91.6794, 92.4643, 93.1870, 93.7495, 94.0401, 94.1492, 93.9694, 93.5836, 93.0964],
                                    index=[28.12, 30.45, 32.78, 35.11, 37.45, 39.78, 42.11, 44.44, 46.77, 49.10, 51.43, 53.76, 56.10, 58.83],
                                    name=200),
                          pd.Series([87.5908, 88.7725, 89.8497, 90.8355, 91.7435, 92.5871, 93.3798, 94.1096, 94.6777, 94.9712, 95.0813, 94.8998, 94.5101, 94.0181],
                                    index=[28.12, 30.45, 32.78, 35.11, 37.45, 39.78, 42.11, 44.44, 46.77, 49.10, 51.43, 53.76, 56.10, 58.83],
                                    name=230)])

    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'])