Energy Storage using the Powell Unified Framework (Part 2)

In Part 1 we used 1-hot encoding for the decision variable. We also optimized a buy-low-sell-high policy from the Policy Function Approximation (PFA) class. In addition, we provided sample-path functionality to the original code. In Part 2 we will maintain the 1-hot encoding but use a different policy - one in the Value Function Approximation (VFA) class. This policy will be based on Bellman’s exact dynamic programming approach making use of Backward Dynamic Programming (BDP).

0 STRUCTURE & FRAMEWORK

The overall structure of this project and report follows the traditional CRISP-DM format. However, instead of the CRISP-DM’S “4 Modeling” section, we inserted the “6 step modeling process” of Dr. Warren Powell in section 4 of this document.

The example explored (but also modified) in this report comes from Dr. Warren Powell (formerly at Princeton). It was chosen as a template to create a POC for a client. It was important to understand this example thoroughly before embarking on creating the POC. Dr Powell’s unified framework shows great promise for unifying the formalisms of at least a dozen different fields. Using his framework enables easier access to thinking patterns in these other fields that might be beneficial and informative to the sequential decision problem at hand. Traditionally, this kind of problem would be approached from the reinforcement learning perspective. However, using Dr. Powell’s wider and more comprehensive perspective almost certainly provides additional value.

Here is information on Dr. Powell’s perspective on Sequential Decision Analytics.

The original code for this example can be found here.

In order to make a strong mapping between the code in this notebook and the mathematics in the unified framework, many variable names in the original code have been adapted. Here is a summary of the mapping between the mathematics and the variable names in the code:

• Superscripts
  • variable names have a double underscore to indicate a superscript
  • \(X^{\pi}\): has code X__pi, is read X pi
• Subscripts
  • variable names have a single underscore to indicate a subscript
  • \(S_t\): has code S_t, is read ‘S at t’
  • \(M^{Spend}_t\) has code M__Spend_t which is read: “MSpend at t”
• Arguments
  • collection variable names may have argument information added
  • \(X^{\pi}(S_t)\): has code X__piIS_tI, is read ‘X pi in S at t’
  • the surrounding I’s are used to imitate the parentheses around the argument
• Next time/iteration
  • variable names that indicate one step in the future are quite common
  • \(R_{t+1}\): has code R_tt1, is read ‘R at t+1’
  • \(R^{n+1}\): has code R__nt1, is read ‘R at n+1’
• Rewards
  • State-independent terminal reward and cumulative reward
    • \(F\): has code F for terminal reward
    • \(\sum_{n}F\): has code cumF for cumulative reward
  • State-dependent terminal reward and cumulative reward
    • \(C\): has code C for terminal reward
    • \(\sum_{t}C\): has code cumC for cumulative reward
• Vectors where components use different names
  • \(S_t(R_t, p_t)\): has code S_t.R_t and S_t.p_t, is read ‘S at t in R at t, and, S at t in p at t’
  • the code implementation is by means of a named tuple
    • self.State = namedtuple('State', SNames) for the ‘class’ of the vector
    • self.S_t for the ‘instance’ of the vector
• Vectors where components reuse names
  • \(x_t(x_{t,GB}, x_{t,BL})\): has code x_t.x_t_GB and x_t.x_t_BL, is read ‘x at t in x at t for GB, and, x at t in x at t for BL’
  • the code implementation is by means of a named tuple
    • self.Decision = namedtuple('Decision', xNames) for the ‘class’ of the vector
    • self.x_t for the ‘instance’ of the vector
• Use of mixed-case variable names
  • to reduce confusion, sometimes the use of mixed-case variable names are preferred (even though it is not a best practice in the Python community), reserving the use of underscores and double underscores for math-related variables

1 BUSINESS UNDERSTANDING

The following business description comes from the free book by Dr. Powell, Sequential Decision Analytics and Modeling:

New Jersey is looking to develop 3,500 megawatts (MW) of offshore wind generating power. A challenge is that wind (and especially offshore wind) can be highly variable, due in part to the property that wind power (over intermediate ranges) increases with the cube of wind speed. This variability is depicted in figure 9.1, developed from a study of offshore wind. Energy from wind power has become popular in regions where wind is high, such as the midwestern United States, coastal regions off of Europe, the northeast of Brazil, and northern regions of China (to name just a few). Sometimes communities (and companies) have invested in renewables (wind or solar) to help reduce their carbon footprint and minimize their dependence on the grid. It is quite rare, however, that these projects allow a community to eliminate the grid from their portfolio. Common practice is to let the renewable source (wind or solar) sell directly to the grid, while a company may purchase from the grid. This can be useful as a hedge since the company will make a lot of money during price spikes (prices may jump from $20 per megawatt-hour (MWh) to $300 per MWh or more) that offsets the cost of purchasing power during those periods.

A challenge with renewables is handling the variability. While one solution is to simply pump any energy from a renewable source into the grid and use the tremendous power of the grid to handle this variability, there has been considerable interest in using storage (in particular, battery storage) to smooth out the peaks and valleys. In addition to handling the variability in the renewable source, there has also been interest in using batteries to take advantage of price spikes, buying power when it is cheap (prices can even go negative) and selling it back when they are high. Exploiting the variability in power prices on the grid is known as battery arbitrage.

This problem will provide insights into virtually any inventory/storage problem, including - Holding cash in mutual funds A bank has to determine how much of its investment capital to hold in cash to meet requests for redemptions, versus investing the money in loans, stocks and bonds. - Retailers (including both online as well as bricksandmortar stores) have to manage inventories of the hundreds of thousands of products. - Auto dealers have to decide how many cars to hold to meet customer demand. - Consulting firms have to decide how many employees to keep on staff to meet the variable demand of different consulting projects.

Energy storage is a particularly rich form of inventory problem.

2 DATA UNDERSTANDING

Next, we will start to look at the data that is provided for this problem.

# import pdb
from collections import namedtuple, defaultdict
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from copy import copy
import time
from scipy.ndimage.interpolation import shift
import pickle
from bisect import bisect
import math
from pprint import pprint
import matplotlib as mpl
pd.options.display.float_format = '{:,.4f}'.format
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', None)
pd.set_option('display.max_colwidth', None)
! python --version
# 3.8.10; found it changed on 3/10/23

Python 3.9.16

DeprecationWarning: Please use `shift` from the `scipy.ndimage` namespace, the `scipy.ndimage.interpolation` namespace is deprecated.
  from scipy.ndimage.interpolation import shift

def process_raw_price_data(file, params):
    DISC_TYPE = "FROM_CUM"
    #DISC_TYPE = "OTHER"

    print("Processing raw price data. Constructing price change list and cdf using {}".format(DISC_TYPE))
    tS = time.time()

    # load energy price data from the Excel spreadsheet
    raw_data = pd.read_excel(file, sheet_name="Raw Data")

    # look at data spanning a week
    data_selection = raw_data.iloc[0:params['T'], 0:5]

    # rename columns to remove spaces (otherwise we can't access them)
    cols = data_selection.columns
    cols = cols.map(lambda x: x.replace(' ', '_') if isinstance(x, str) else x)
    data_selection.columns = cols

    # sort prices in ascending order
    sort_by_price = data_selection.sort_values('PJM_RT_LMP')
    #print(sort_by_price.head())

    hist_price = np.array(data_selection['PJM_RT_LMP'].tolist())
    #print(hist_price[0])

    max_price = pd.DataFrame.max(sort_by_price['PJM_RT_LMP'])
    min_price = pd.DataFrame.min(sort_by_price['PJM_RT_LMP'])
    print("Min price {:.2f} and Max price {:.2f}".format(min_price,max_price))
    
    # sort prices in ascending order
    sort_by_price = data_selection.sort_values('PJM_RT_LMP')

    # calculate change in price and sort values of change in price in ascending order
    data_selection['Price_Shift'] = data_selection.PJM_RT_LMP.shift(1)
    data_selection['Price_Change'] = data_selection['PJM_RT_LMP'] - data_selection['Price_Shift']
    sort_price_change = data_selection.sort_values('Price_Change')
    
    # discretize change in price and obtain f(p) for each price change
    max_price_change = pd.DataFrame.max(sort_price_change['Price_Change'])
    min_price_change = pd.DataFrame.min(sort_price_change['Price_Change'])
    print("Min price change {:.2f} and Max price change {:.2f}".format(min_price_change,max_price_change))
    
    # there are 191 values for price change
    price_changes_sorted = sort_price_change['Price_Change'].tolist()
    # remove the last NaN value
    price_changes_sorted.pop()

    if DISC_TYPE == "FROM_CUM":
    # discretize price change  by interpolating from cumulative distribution
        xp = price_changes_sorted
        fp = np.arange(len(price_changes_sorted) - 1) / (len(price_changes_sorted) - 1)
        cum_fn = np.append(fp, 1)

        # obtain 30 discrete prices
        discrete_price_change_cdf = np.linspace(0, 1, params['nPriceChangeInc'])
        discrete_price_change_list = []
        for i in discrete_price_change_cdf:
            interpolated_point = np.interp(i, cum_fn, xp)
            discrete_price_change_list.append(interpolated_point)
    else:
        price_change_range = max_price_change - min_price_change
        price_change_increment = price_change_range / params['nPriceChangeInc']
        discrete_price_change = np.arange(min_price_change, max_price_change, price_change_increment)
        discrete_price_change_list = list(np.append(discrete_price_change, max_price_change))

        f_p = np.arange(len(price_changes_sorted) - 1) / (len(price_changes_sorted) - 1)
        cum_fn = np.append(f_p, 1)
        discrete_price_change_cdf = []
        for c in discrete_price_change_list:
            interpolated_point = np.interp(c, price_changes_sorted, cum_fn)
            discrete_price_change_cdf.append(interpolated_point)

    price_changes_sorted = np.array(price_changes_sorted)
    discrete_price_change_list = np.array(discrete_price_change_list)
    discrete_price_change_cdf = np.array(discrete_price_change_cdf)
    discrete_price_change_pdf = discrete_price_change_cdf - shift(discrete_price_change_cdf,1,cval=0)

    mean_price_change = np.dot(discrete_price_change_list,discrete_price_change_pdf)

    #print("discrete_price_change_list ", discrete_price_change_list)
    #print("discrete_price_change_cdf", discrete_price_change_cdf)
    #print("discrete_price_change_pdf", discrete_price_change_pdf)

    print("Finishing processing raw price data in {:.2f} secs. Expected price change is {:.2f}. Hist_price len is {}".format(time.time()-tS,mean_price_change,len(hist_price)))
    #input("enter any key to continue...") 

    exog_params = {
        "hist_price": hist_price,
        "price_changes_sorted": price_changes_sorted,
        "discrete_price_change_list": discrete_price_change_list,
        "discrete_price_change_cdf": discrete_price_change_cdf}
    return exog_params

seed = 189654913
file = 'Parameters.xlsx'

# NOTE:
# R__max: unit is MWh which is an energy unit
# R_0:    unit is MWh which is an energy unit
parDf = pd.read_excel(f'{base_dir}/{file}', sheet_name='ParamsModel', index_col=0); print(f'{parDf}')
# parDict = parDf.set_index('Index').T.to_dict('list')#.
parDict = parDf.T.to_dict('list') #.
params = {key:v for key, value in parDict.items() for v in value}
params['seed'] = seed
params['T'] = min(params['T'], 192); print(f'{params=}')

                    0
Index                
Algorithm  BackwardDP
T                 195
eta            0.9000
R__max              1
R_0                 1
params={'Algorithm': 'BackwardDP', 'T': 192, 'eta': 0.9, 'R__max': 1, 'R_0': 1, 'seed': 189654913}

parDf = pd.read_excel(f'{base_dir}/{file}', sheet_name='GridSearch', index_col=0); print(parDf)
# parDict = parDf.set_index('Index').T.to_dict('list')#.
parDict = parDf.T.to_dict('list')
paramsPolicy = {key:v for key, value in parDict.items() for v in value}; print(f'{paramsPolicy=}')
params.update(paramsPolicy); pprint(f'{params=}')

                  0
Index              
theta_sell_min   10
theta_sell_max  100
theta_buy_min    10
theta_buy_max   100
theta_inc         1
paramsPolicy={'theta_sell_min': 10, 'theta_sell_max': 100, 'theta_buy_min': 10, 'theta_buy_max': 100, 'theta_inc': 1}
("params={'Algorithm': 'BackwardDP', 'T': 192, 'eta': 0.9, 'R__max': 1, 'R_0': "
 "1, 'seed': 189654913, 'theta_sell_min': 10, 'theta_sell_max': 100, "
 "'theta_buy_min': 10, 'theta_buy_max': 100, 'theta_inc': 1}")

parDf = pd.read_excel(f'{base_dir}/{file}', sheet_name='BackwardDP', index_col=0); print(parDf)
# parDict = parDf.set_index('Index').T.to_dict('list')#.
parDict = parDf.T.to_dict('list')
paramsPolicy = {key:v for key, value in parDict.items() for v in value}; print(f'{paramsPolicy=}')
params.update(paramsPolicy); pprint(f'{params=}')

                      0
Index                  
nPriceChangeInc      50
priceDiscSet     1, 0.5
run3D             False
paramsPolicy={'nPriceChangeInc': 50, 'priceDiscSet': '1, 0.5', 'run3D': False}
("params={'Algorithm': 'BackwardDP', 'T': 192, 'eta': 0.9, 'R__max': 1, 'R_0': "
 "1, 'seed': 189654913, 'theta_sell_min': 10, 'theta_sell_max': 100, "
 "'theta_buy_min': 10, 'theta_buy_max': 100, 'theta_inc': 1, "
 "'nPriceChangeInc': 50, 'priceDiscSet': '1, 0.5', 'run3D': False}")

if isinstance(params['priceDiscSet'], str):
    price_disc_list = params['priceDiscSet'].split(",")
    price_disc_list = [float(e) for e in price_disc_list]
else:
    price_disc_list = [float(params['priceDiscSet'])] 
params['price_disc_list'] = price_disc_list; print(f'{price_disc_list=}')

price_disc_list=[1.0, 0.5]

pprint(f"{params=}")
#input("enter any key to continue...")

("params={'Algorithm': 'BackwardDP', 'T': 192, 'eta': 0.9, 'R__max': 1, 'R_0': "
 "1, 'seed': 189654913, 'theta_sell_min': 10, 'theta_sell_max': 100, "
 "'theta_buy_min': 10, 'theta_buy_max': 100, 'theta_inc': 1, "
 "'nPriceChangeInc': 50, 'priceDiscSet': '1, 0.5', 'run3D': False, "
 "'price_disc_list': [1.0, 0.5]}")

#exog_params is a dictionary with  three lists: 
#   hist_price, price_changes_list, discrete_price_change_cdf
# exog_params = process_raw_price_data(f'{base_dir}/{file}', params)
exogParams = process_raw_price_data(f'{base_dir}/{file}', params) #.
# print(f"{exogParams['hist_price']=}")
# print(f"{exogParams['price_changes_sorted']=}")
# print(f"{exogParams['discrete_price_change_list']=}")
# print(f"{exogParams['discrete_price_change_cdf']=}")

Processing raw price data. Constructing price change list and cdf using FROM_CUM
Min price 0.00 and Max price 107.34
Min price change -54.40 and Max price change 65.81
Finishing processing raw price data in 1.76 secs. Expected price change is 1.28. Hist_price len is 192

print("Ingesting raw price data for understanding.")
# load energy price data from the Excel spreadsheet
raw_data = pd.read_excel(f'{base_dir}/{file}', sheet_name="Raw Data")

# look at data spanning a week
df_raw = raw_data.iloc[0:params['T'], 0:5]

# rename columns to remove spaces (otherwise we can't access them)
cols = df_raw.columns
cols = cols.map(lambda x: x.replace(' ', '_') if isinstance(x, str) else x)
df_raw.columns = cols
df_raw.head()

Ingesting raw price data for understanding.

	Time	Day	Hour	PJM_DA_LMP	PJM_RT_LMP
0	1.0000	2005-01-01 00:00:00	1.0000	25.7200	18.9100
1	2.0000	2005-01-01 01:00:00	2.0000	23.0800	14.2700
2	3.0000	2005-01-01 02:00:00	3.0000	21.0000	9.6700
3	4.0000	2005-01-01 03:00:00	4.0000	20.0000	17.6800
4	5.0000	2005-01-01 04:00:00	5.0000	20.0400	11.8800

from certifi.core import where
def plot_data(df):
  n_charts = 2
  ylabelsize = 16
  fig, axs = plt.subplots(n_charts, sharex=True)
  fig.set_figwidth(10); fig.set_figheight(7)
  fig.suptitle('PRICE OF ELECTRICITY [$/MWh]', fontsize=20)
  i = 0
  axs[i].set_ylim(auto=True); axs[i].spines['top'].set_visible(False); axs[i].spines['right'].set_visible(True); axs[i].spines['bottom'].set_visible(True)
  axs[i].plot(df['PJM_DA_LMP'], 'r')
  axs[i].set_ylabel('$\mathrm{PJM\ DA\ LMP}$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);
  i = 1
  axs[i].set_ylim(auto=True); axs[i].spines['top'].set_visible(False); axs[i].spines['right'].set_visible(True); axs[i].spines['bottom'].set_visible(True)
  axs[i].plot(df['PJM_RT_LMP'], 'g')
  axs[i].set_ylabel('$\mathrm{PJM\ RT\ LMP}$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);
  axs[i].set_xlabel('$t\ \mathrm{[hours]}$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);
plot_data(df_raw)  

3 DATA PREPARATION

The most important aspect of data preparation will be to create an empirical distribution rather than to try and fit the data to a parametric distribution. This involves the computation of a cumulative distribution from the data. We first sort the prices from smallest to largest. Let this ordered sequence be represented by \(\tilde{p}_t\) where \(\tilde{p}_{t-1} \le \tilde{p}_t\). Also, let \(T\) be the number of time periods. The percentage of time periods with a price less than \(\tilde{p}_t\) is then \(t/T\). The empirical cumulative distribution is given by \[ F_P(\tilde{p}_t) = \frac{t}{T} \]

# hide
df_hist_price = pd.DataFrame({'hist_price': exogParams['hist_price']});
print(df_hist_price.shape)
# print(df_hist_price)
# df_hist_price.plot('hist_price');
# df_hist_price.plot();

(192, 1)

df_price_changes_sorted = pd.DataFrame({'price_changes_sorted': exogParams['price_changes_sorted']});
print(df_price_changes_sorted.shape)
# print(df_price_changes_sorted)
df_price_changes_sorted.plot(figsize=(10,5));

(191, 1)

df_discrete_price_change_list = pd.DataFrame(
    {'discrete_price_change_list': exogParams['discrete_price_change_list']})
print(df_discrete_price_change_list.shape)
# print(df_discrete_price_change_list)
df_discrete_price_change_list.plot(figsize=(10,5));

(50, 1)

df_discrete_price_change_cdf = pd.DataFrame(
    {'discrete_price_change_cdf': exogParams['discrete_price_change_cdf']})
print(df_discrete_price_change_cdf.shape)
# print(df_discrete_price_change_cdf)
# df_discrete_price_change_cdf.plot(figsize=(10,5));

(50, 1)

4 MODELING

4.1 Narrative

There has been an interest making use of batteries to exploit electricity price spikes. Power is bought when the price is low and power is sold when the price is high. This is known as battery arbitrage.

This project hightlights some interesting aspects of this kind of problem. Among these are:

• Electricity prices on the grid can be very volatile (varies between about $20/MWh to about $1,000/MWh or even $10,000/MWh for brief periods)
• Prices change every 5 minutes
• We assume that we can observe the price and then decide if we want to buy or sell
• We can only buy or sell for an entire 5-minute interval at a maximum rate of 0.01 MW
• Wind energy can be forecasted but not very well. Rolling forecasts can be used.
• Energy demand is variable but quite predictable. It depends on:
  • temperature (primarily)
  • humidity (less)
• Energy may be bought from or sold to the grid at current grid prices.
• There is a 5% to 10% loss when converting power from AC (grid) to DC (to store in battery).

4.2 Core Elements

This section attempts to answer three important questions: - What metrics are we going to track? - What decisions do we intend to make? - What are the sources of uncertainty?

For this problem, the only metric we are interested in is the amount of profit we make after each decision window. A single type of decision needs to be made at the start of each window - how much energy to buy or sell. The only source of uncertainty is the price of energy.

4.3 Mathematical Model | SUM Design

A Python class is used to implement the model for the SUM (System Under Management):

class EnergyStorageModel():
        def __init__(
        self, SVarNames, xVarNames, S_0, params, exogParams, possibleDecisions,
        W_fn=None, S__M_fn=None, C_fn=None):
        ...
        ...

4.3.1 State variables

The state variables represent what we need to know.

• \(R_t\)
  • amount of energy stored in the battery at time \(t\)
  • measured in kWh or MWh
• \(p_t\)
  • price of energy on the grid

The state is:

\(S_t = (R_t,p_t)\)

The state variables are represented by the following variables in the EnergyStorageModel class:

        self.SVarNames = SVarNames
        self.State = namedtuple('State', SVarNames) # 'class'
        self.S_t = self.build_state(self.S_0) # 'instance'

where

SVarNames = ['R_t', 'p_t']

4.3.2 Decision variables

The decision variables represent what we control.

• \(x_t\)
  • amount of energy bought from (\(x_t>0\)) or sold to (\(x_t<0\)) the grid
• Losses
  • When energy is transferred to or from the battery a loss is incurred
  • Assume the loss into the battery is the same as the loss out of the battery
  • \(\eta\) is the fraction of the energy retained
  • \(1 - \eta\) is the loss
• Constraints
  • buy from grid \((x_t>0)\):
    • \(x_t \le \frac{1}{\eta}(R^{max} - R_t)\)
  • sell to grid \((x_t<0)\):
    • \(-x_t \le R_t\), which comes to
    • \(x_t \ge -R_t\)
• Decisions are made with a policy (TBD below):
  • \(X^{\pi}(S_t)\)

The decision variables are represented by the following variables in the EnergyStorageModel class:

self.Decision = namedtuple('Decision', xVarNames) # 'class'

where

xVarNames = ['x_t_buy', 'x_t_hold', 'x_t_sell']

**NOTE: We use a 1-hot representation of the decision variable. For this reason, some of our method names will have a suffix ’_1hot’ to point this out. In followup notebooks we will move away from the 1-hot representation.**

4.3.3 Exogenous information variables

The exogenous information variables represent what we did not know (when we made a decision). These are the variables that we cannot control directly. The information in these variables become available after we make the decision \(x_t\).

When we assume that the price in each time period is revealed, without any model to predict the price based on past prices, we have:

\[ W_{t+1} = p_{t+1} \]

Alternatively, when we assume that we observe the change in price \(\hat{p}=p_t-p_{t-1}\), we have:

\[ W_{t+1} = \hat{p}_{t+1} \]

The exogenous information have been preprocessed and is captured in:

• exogParams['hist_price']
• exogParams['price_changes_sorted']
• exogParams['discrete_price_change_list']
• exogParams['discrete_price_change_cdf']

The latest exogenous information can be accesses by calling the following method from class EnergyStorageModel():

def W_fn(self, t):
  W_tt1 = self.exogParams['hist_price'][t]
  return W_tt1

4.3.4 Transition function

The transition function describe how the state variables evolve over time. Because we currently have two state variables in the state, \(S_t=(R_t,p_t)\), we have the equations:

\[ \begin{aligned} R_{t+1} &= R_t + \eta x_t \\ p_{t+1} &= p_t + \hat{p}_{t+1} \end{aligned} \]

Collectively, they represent the general transition function:

\[ S_{t+1} = S^M(S_t,X^{\pi}(S_t)) \]

The transition function is implemented by the following method in class EnergyStorageModel():

def S__M_fn_1hot(self, t, x_t): #.
    p_tt1 = self.W_fn(t)
    R_tt1 = self.S_t.R_t + (self.initArgs['eta']*x_t.x_t_buy) - x_t.x_t_sell
    if len(self.SVarNames) == 2:
        S_tt1 = self.build_state({'R_t': R_tt1, 'p_t': p_tt1})
    elif len(self.SVarNames) == 3:
        S_tt1 = self.build_state({'R_t': R_tt1, 'p_t': p_tt1, 'p_t_1': self.S_t.p_t})
    return S_tt1

4.3.5 Objective function

The objective function captures the performance metrics of the solution to the problem. It is given by:

\[ \max_{\pi} \mathbb{E} \sum_{t=0}^T C(S_t,x_t) \]

where \[ C(S_t,x_t) = -p_tx_t \]

The contribution (reward) function is implemented by the following method in class EnergyStorageModel():

def C_fn_1hot(self, x_t):
    C_t = self.S_t.p_t*( self.initArgs['eta']*x_t.x_t_sell - x_t.x_t_buy )#. C(S_t, x_t, W_t+1)
    return C_t

4.3.6 Implementation of SUM Model

Here is the complete implementation of the EnergyStorageModel class:

class EnergyStorageModel():
    def __init__(
        self, SVarNames, xVarNames, S_0, params, exogParams, possibleDecisions,
        W_fn=None, S__M_fn=None, C_fn=None):
        self.initArgs = params
        self.prng = np.random.RandomState(params['seed'])
        self.exogParams = exogParams
        self.S_0 = S_0
        self.SVarNames = SVarNames
        self.xVarNames = xVarNames
        self.possibleDecisions = possibleDecisions
        self.State = namedtuple('State', SVarNames) # 'class'
        self.S_t = self.build_state(self.S_0) # 'instance'
        self.Decision = namedtuple('Decision', xVarNames) # 'class'
        self.cumC = 0.0 #. cumulative reward; use F or cumF for final (i.e. no cumulative) reward
        # self.states = [self.S_t] # keep a list of states visited

    def reset(self):
        self.cumC = 0.0
        self.S_t = self.build_state(self.S_0)
        # self.states = [self.S_t]

    def build_state(self, info):
        return self.State(*[info[sn] for sn in self.SVarNames])

    def build_decision_1hot(self, x_t_info, R_t):
        x_t_info_constrained = {'buy': 0, 'hold': 0, 'sell': 0} #.
        #amount of power that can be bought/sold is limited by constraints
        for xn in self.xVarNames:
            if xn == 'x_t_buy' and x_t_info[xn] > 0:
                x_t_info_constrained[xn] = (self.initArgs['R__max'] - R_t)/self.initArgs['eta']
            elif xn == 'x_t_sell' and x_t_info[xn] > R_t:
                x_t_info_constrained[xn] = R_t
            else: # 'x_t_hold'
                x_t_info_constrained[xn] = x_t_info[xn]
        return self.Decision(*[x_t_info_constrained[xn] for xn in self.xVarNames])

    def W_fn(self, t):
        W_tt1 = self.exogParams['hist_price'][t]
        return W_tt1

    def S__M_fn_1hot(self, t, x_t):
        p_tt1 = self.W_fn(t)
        R_tt1 = self.S_t.R_t + (self.initArgs['eta']*x_t.x_t_buy) - x_t.x_t_sell
        if len(self.SVarNames) == 2:
            S_tt1 = self.build_state({'R_t': R_tt1, 'p_t': p_tt1})
        elif len(self.SVarNames) == 3:
            S_tt1 = self.build_state({'R_t': R_tt1, 'p_t': p_tt1, 'p_t_1': self.S_t.p_t})
        return S_tt1

    def C_fn_1hot(self, x_t):
        C_t = self.S_t.p_t*( self.initArgs['eta']*x_t.x_t_sell - x_t.x_t_buy )# C(S_t, x_t, W_t+1)
        return C_t

    def step_1hot(self, t, x_t):
        self.cumC += self.C_fn_1hot(x_t)
        self.S_t = self.S__M_fn_1hot(t, x_t)
        return (self.S_t, self.cumC, x_t) #. for plotting

4.4 Uncertainty Model

As mentioned above in section 3, we make use of an empirical distribution rather than to try and fit the data to a parametric distribution. This involves the computation of a cumulative distribution from the data. We first sort the prices from smallest to largest. Let this ordered sequence be represented by \(\tilde{p}_t\) where \(\tilde{p}_{t-1} \le \tilde{p}_t\). Also, let \(T\) be the number of time periods. The percentage of time periods with a price less than \(\tilde{p}_t\) is then \(t/T\). The empirical cumulative distribution is given by \[ F_P(\tilde{p}_t) = \frac{t}{T} \]

4.5 Policy Design

There are two main meta-classes of policy design. Each of these has two subclasses: - Policy Search - Policy Function Approximations (PFAs) - Cost Function Approximations (CFAs) - Lookahead - Value Function Approximations (VFAs) - Direct Lookaheads (DLAs)

In this project we will use: - A Belmann type policy, making use of Backward Dynamic Programming (BDP). This policy falls in the VFA class.

The Belmann policy is implemented by the following method in class EnergyStoragePolicy():

def X__bellman_1hot(self, t, S_t, bellmanModel, T): #. SDAM-9.4.2
    p_t = S_t.p_t
    R_t = S_t.R_t
    vMax = -np.inf
    dMax = None
    for d in self.model.possibleDecisions: #. for each x_t
        x_t = self.model.build_decision_1hot(d, R_t)
        C_t = p_t*(x_t.x_t_sell - x_t.x_t_buy)
        Sum_w = 0
        w_idx = 0
        for w in bellmanModel.discretePriceChanges: #. for each w
            f__W = bellmanModel.f_p[w_idx] if w_idx == 0 else bellmanModel.f_p[w_idx] - bellmanModel.f_p[w_idx - 1]
            S_tt1 = bellmanModel.S__M(S_t, x_t, w)
            V_tt1 = bellmanModel.V_tIS_tI[t+1][S_tt1] if t < bellmanModel.time \
                else bellmanModel.terminalContribution
            Sum_w += f__W*V_tt1
            w_idx += 1
        # print(f"{w_idx=}")
        v = C_t + Sum_w
        if (v > vMax):
            vMax = v
            dMax = d
    return dMax

4.5.1 Implementation of Policy Design

The EnergyStoragePolicy() class implements the policy design.

import random
class EnergyStoragePolicy():
    def __init__(self, model, piNames):
        self.model = model
        self.piNames = piNames
        self.Policy = namedtuple('Policy', piNames)

    def X__bellman_1hot(self, t, S_t, bellmanModel, T): #. SDAM-9.4.2
        p_t = S_t.p_t
        R_t = S_t.R_t
        vMax = -np.inf
        dMax = None
        for d in self.model.possibleDecisions: #. for each x_t
            x_t = self.model.build_decision_1hot(d, R_t)
            C_t = p_t*(x_t.x_t_sell - x_t.x_t_buy)
            Sum_w = 0
            w_idx = 0
            for w in bellmanModel.discretePriceChanges: #. for each w
                f__W = bellmanModel.f_p[w_idx] if w_idx == 0 else bellmanModel.f_p[w_idx] - bellmanModel.f_p[w_idx - 1]
                S_tt1 = bellmanModel.S__M_1hot(S_t, x_t, w)
                V_tt1 = bellmanModel.V_tIS_tI[t+1][S_tt1] if t < bellmanModel.time \
                    else bellmanModel.terminalContribution
                Sum_w += f__W*V_tt1
                w_idx += 1
            # print(f"{w_idx=}")
            v = C_t + Sum_w
            if (v > vMax):
                vMax = v
                dMax = d
        return dMax
        
    def run_policy_1hot(self, piInfo, piName, params):
        model_copy = copy(self.model)
        T = params['T']
        # nTrades = {'buy': 0, 'sell': 0, 'hold': 0}
        # nTrades = {'x_t_buy': 0, 'x_t_sell': 0, 'x_t_hold': 0}
        buyList = []
        sellList = []
        for t in range(T): #for each transition/step
            x_t_info = getattr(self, piName)(t, model_copy.S_t, piInfo, T)
            if t == model_copy.initArgs['T'] - 1: #Last time period - sell
                x_t_info = {'x_t_buy': 0, 'x_t_hold': 0, 'x_t_sell': 1}   
            x_t = model_copy.build_decision_1hot(x_t_info, model_copy.S_t.R_t)#.

            # nTrades['x_t_buy'] += x_t.x_t_buy; #print(f'### {x_t.buy=}')
            # nTrades['x_t_sell'] += x_t.x_t_sell; #print(f'### {x_t.sell=}')
            # nTrades['x_t_hold'] += model_copy.S_t.R_t
            if x_t.x_t_buy > 0:
                buyList.append((t, model_copy.S_t.R_t, model_copy.S_t.p_t, x_t.x_t_buy, params['R__max']))
            elif x_t.x_t_sell > 0:
                sellList.append((t, model_copy.S_t.R_t, model_copy.S_t.p_t, x_t.x_t_sell, params['R__max']))
            _, _, _ = model_copy.step_1hot(t, x_t) # step the model forward
        cumC = model_copy.cumC
        
        # print(f"Energy traded - Buy: {nTrades['x_t_buy']:.2f} - Hold: {nTrades['x_t_hold']:.2f} - Sell: {nTrades['x_t_sell']:.2f}")
        # print("Buy times and prices ")
        # for i in range(len(buyList)):
        #     print(f"t = {buyList[i][0]:.2f} and price = {buyList[i][1]:.2f}") #.
        # print("Sell times and prices ")
        # for i in range(len(sellList)):
        #     print(f"t = {sellList[i][0]:.2f} and price = {sellList[i][1]:.2f}") #.
        return cumC

    def run_policy_sample_paths_1hot(self, theta, piName, params):
        FhatIomega__lI = []
        for l in range(1, params['L'] + 1): #for each sample-path
            model_copy = copy(self.model)
            record_l = [piName, theta, l]
            T = params['T']
            # nTrades = {'x_t_buy': 0, 'x_t_sell': 0, 'x_t_hold': 0}
            buyList = []
            sellList = []
            for t in range(T): #for each transition/step
                x_t_info = getattr(self, piName)(t, model_copy.S_t, theta)
                if t == model_copy.initArgs['T'] - 1: #Last time period - sell
                    x_t_info = {'x_t_buy': 0, 'x_t_hold': 0, 'x_t_sell': 1}
                x_t = model_copy.build_decision_1hot(x_t_info, model_copy.S_t.R_t)
                # nTrades['x_t_buy'] += x_t.x_t_buy; #print(f'### {x_t.buy=}')
                # nTrades['x_t_sell'] += x_t.x_t_sell; #print(f'### {x_t.sell=}')
                # nTrades['x_t_hold'] += model_copy.S_t.R_t
                if x_t.x_t_buy > 0:
                    buyList.append((t, model_copy.S_t.R_t, model_copy.S_t.p_t, x_t.x_t_buy, params['R__max']))
                elif x_t.x_t_sell > 0:
                    sellList.append((t, model_copy.S_t.R_t, model_copy.S_t.p_t, x_t.x_t_sell, params['R__max']))
                _, _, _ = model_copy.step_1hot(t, x_t)
            FhatIomega__lI.append(model_copy.cumC) #. just above (SDAM-eq2.9) #. Fhat for this sample-path is in model_copy.cumC
            # print(f"Inventory traded - Buy: {nTrades['x_t_buy']:.2f} - Hold: {nTrades['x_t_hold']:.2f} - Sell: {nTrades['x_t_sell']:,}")
            # print("Buy times ")
            # for i in range(len(buy_list)):
            #     print(f"t={buy_list[i][0]}, R_t={buy_list[i][1]:.2f}, D_t={buy_list[i][2]:.2f}, x_t.x_t_buy={buy_list[i][3]:.2f}, R__max={buy_list[i][4]:,}")
            # print("Sell times ")
            # for i in range(len(sell_list)):
            #     print(f"t={sell_list[i][0]}, R_t={sell_list[i][1]:.2f}, D_t={sell_list[i][2]:.2f}, x_t.x_t_sell={sell_list[i][3]:.2f}, R__max={buy_list[i][4]:,}")
            # return cumC
        # end L
        return FhatIomega__lI

    def perform_bellman_search_1hot(self, params, discreteEnergy, minPrice, maxPrice, incValues): #.
        tS = time.time()
        cumCI_inc_I = {}
        bdpI_inc_I = {}
        best_inc = -np.inf
        best_cumC = -np.inf
        best_bdp = None
        for inc in incValues: #for each discretization
            discretePrices = np.arange(minPrice, maxPrice + inc, inc)
            print(f"\nStarting BackwardDP 2D for price discretization with inc={inc}")
            test_2D_bdp = BDP(
                discretePrices=discretePrices, 
                discreteEnergy=discreteEnergy, 
                priceChanges=exogParams['price_changes_sorted'], 
                discretePriceChanges=exogParams['discrete_price_change_list'], 
                f_p=exogParams['discrete_price_change_cdf'], 
                stop_time=params['T'], 
                model=copy(M))
            # 2D states - time the process with a 2D state variable
            t0 = time.time()
            V_tIS_tI = test_2D_bdp.bellman_1hot()
            t1 = time.time()
            time_elapsed = t1 - t0
            print("Time_elapsed_2D_model={:.2f} secs.".format(time_elapsed))

            print("Starting policy evaluation for the actual sample path")
            tS = time.time()
            cumC = self.run_policy_1hot(test_2D_bdp, "X__bellman_1hot", params)
            cumCI_inc_I[inc] = cumC
            bdpI_inc_I[inc] = test_2D_bdp
            if(cumC > best_cumC):
                best_cumC = cumC
                best_inc = inc
                best_bdp = test_2D_bdp
            print(f"Finishing inc {inc} with cumC {cumC:.2f}. Best inc so far {best_inc}. Best contribution {cumCI_inc_I[best_inc]:.2f}")
        print(f"Finishing BellmanSearch in {time.time() - tS:.2f} secs")
        return cumCI_inc_I, best_bdp, bdpI_inc_I

    def plot_Fhat_map(self, FhatI_theta_I, thetasX, thetasY, labelX, labelY, title):
        Fhat_values = [FhatI_theta_I[(thetaX,thetaY)]  for thetaY in thetasY for thetaX in thetasX]
        Fhats = np.array(Fhat_values)
        increment_count = len(thetasX)
        Fhats = np.reshape(Fhats, (-1, increment_count))

        fig, ax = plt.subplots()
        im = ax.imshow(Fhats, cmap='hot', origin='lower', aspect='auto')
        # create colorbar
        cbar = ax.figure.colorbar(im, ax=ax)
        # cbar.ax.set_ylabel(cbarlabel, rotation=-90, va="bottom")
        # we want to show all ticks...
        ax.set_xticks(np.arange(0,len(thetasX), 5))
        ax.set_yticks(np.arange(0,len(thetasY), 5))
        # ... and label them with the respective list entries
        ax.set_xticklabels(thetasX[::5])
        ax.set_yticklabels(thetasY[::5])
        # rotate the tick labels and set their alignment.
        #plt.setp(ax.get_xticklabels(), rotation=45, ha="right",rotation_mode="anchor")

        ax.set_title(title)

        ax.set_xlabel(labelX)
        ax.set_ylabel(labelY)

        #fig.tight_layout()
        plt.show()
        return True

    def plot_Fhat_chart(self, FhatI_theta_I, thetasX, labelX, labelY, title, color_style):
        mpl.rcParams['lines.linewidth'] = 1.2
        xylabelsize = 18
        plt.figure(figsize=(25, 8))
        plt.title(title, fontsize=20)
        Fhats = FhatI_theta_I.values()
        plt.plot(thetasX, Fhats, color_style)
        plt.xlabel(labelX, rotation=0, ha='right', va='center', fontweight='bold', size=xylabelsize)
        plt.ylabel(labelY, rotation=0, ha='right', va='center', fontweight='bold', size=xylabelsize)
        plt.show()        

class BDP(): #. SDAM-9.4.2: Backward Dynamic Programming; SDAM-fig9.6; RLSO-fig14.3
    def __init__(self, discretePrices, discreteEnergy, priceChanges, discretePriceChanges,
                 f_p, stop_time, model):
        self.discreteEnergy = discreteEnergy
        self.discretePrices = discretePrices
        self.priceChanges = priceChanges
        self.discretePriceChanges = discretePriceChanges
        self.f_p = f_p
        self.time = stop_time - 1
        self.model = model
        self.terminalContribution = 0
        self.V_tIS_tI = None #this will store the vfas - it will be computed by the method bellman_2D or bellman_3D

    def S__M_1hot(self, S_t, x_t, exogInfo):
        R_tt1 = S_t.R_t + (self.model.initArgs['eta']*x_t.x_t_buy) - x_t.x_t_sell
        R_tt1_adjusted = math.ceil(R_tt1)

        if len(S_t) == 2:
            p_tt1 = S_t.p_t + exogInfo
        elif len(S_t) == 3:
            p_tt1 = 0.5*S_t.p_t_1 + 0.5*S_t.p_t + exogInfo
        
        if p_tt1 <= min(self.discretePrices):
            p_tt1_adjusted = min(self.discretePrices)
        elif p_tt1 >= max(self.discretePrices):
            p_tt1_adjusted = max(self.discretePrices)
        else:
            index = bisect(self.discretePrices, p_tt1)
            p_tt1_adjusted = self.discretePrices[index]

        if len(S_t) == 2:
            S_tt1 = self.model.build_state({'R_t': R_tt1_adjusted, 'p_t': p_tt1_adjusted})
        elif len(S_t) == 3:
            p_t = S_t.p_t #.
            if p_t <= min(self.discretePrices):
                p_t_adjusted = min(self.discretePrices)
            elif p_t >= max(self.discretePrices):
                p_t_adjusted = max(self.discretePrices)
            else:
                index = bisect(self.discretePrices, p_t)
                p_t_adjusted = self.discretePrices[index]
            S_tt1 = self.model.build_state({
                'R_t': R_tt1_adjusted,
                'p_t': p_tt1_adjusted,
                'p_t_1': p_t_adjusted
            })
        return S_tt1

    def S__M(self, S_t, x_t, exogInfo):
        R_tt1 = self.model.S_t.R_t + self.model.initArgs['eta']*x_t #. x_t here is number
        R_tt1_adjusted = math.ceil(R_tt1)

        if len(S_t) == 2:
            p_tt1 = S_t.p_t + exogInfo
        elif len(S_t) == 3:
            p_tt1 = 0.5*S_t.p_t_1 + 0.5*S_t.p_t + exogInfo
        
        if p_tt1 <= min(self.discretePrices):
            p_tt1_adjusted = min(self.discretePrices)
        elif p_tt1 >= max(self.discretePrices):
            p_tt1_adjusted = max(self.discretePrices)
        else:
            index = bisect(self.discretePrices, p_tt1)
            p_tt1_adjusted = self.discretePrices[index]

        if len(S_t) == 2:
            S_tt1 = self.model.build_state({'R_t': R_tt1_adjusted, 'p_t': p_tt1_adjusted})
        elif len(S_t) == 3:
            p_t = S_t.p_t
            if p_t <= min(self.discretePrices):
                p_t_adjusted = min(self.discretePrices)
            elif p_t >= max(self.discretePrices):
                p_t_adjusted = max(self.discretePrices)
            else:
                index = bisect(self.discretePrices, p_t)
                p_t_adjusted = self.discretePrices[index]
            S_tt1 = self.model.build_state({
                'R_t': R_tt1_adjusted,
                'p_t': p_tt1_adjusted,
                'p_t_1': p_t_adjusted
            })
        return S_tt1

    def bellman_1hot(self):
        # make list of all possible 2D states using discretized values
        self.possibleStates = []
        if len(self.model.SVarNames) == 2:
            for price in self.discretePrices:
                for energy in self.discreteEnergy:
                    S_t = self.model.build_state({'R_t': energy, 'p_t': price})
                    self.possibleStates.append(S_t)
        else:
            for p in self.discretePrices:
                for prev_p in self.discretePrices:
                    for energy in self.discreteEnergy:
                        S_t = self.model.build_state({'R_t': energy, 'p_t': p, 'p_t_1': prev_p})
                        self.possibleStates.append(S_t)

        print("State dimension: {}. State space size: {}. Exogenous info size: {}".format(len(self.model.SVarNames),len(self.possibleStates),len(self.discretePriceChanges)))
        print(f'{self.discretePrices=}')
        print(f'{self.discreteEnergy=}')
        # print(f'{self.possibleStates=}')

        t = self.time
        V_tIS_tI = defaultdict(dict)
        while t != -1: #. for each t
            maxDict = {}
            for S_t in self.possibleStates: #. for each S_t
                p_t = S_t.p_t
                R_t = S_t.R_t
                vList = []
                for d in self.model.possibleDecisions: #. for each x_t
                    x_t = self.model.build_decision_1hot(d, R_t)
                    C_t = p_t*(self.model.initArgs['eta']*x_t.x_t_sell - x_t.x_t_buy)
                    Sum_w = 0
                    w_idx = 0
                    for w in self.discretePriceChanges: #. for each w
                        f__W = self.f_p[w_idx] if w_idx == 0 else self.f_p[w_idx] - self.f_p[w_idx - 1]
                        S_tt1 = self.S__M_1hot(S_t, x_t, w); #print(f'### {next_state=}')
                        V_tt1 = V_tIS_tI[t + 1][S_tt1] if t < self.time \
                            else self.terminalContribution
                        Sum_w += f__W*V_tt1
                        w_idx += 1
                    v = C_t + Sum_w
                    vList.append(v)
                maxValue = max(vList)
                # decList = ["Buy", "Sell", "Hold"]
                # print(f"Time: {t} State: p_t={p_t:.2f}, R_t={R_t:.2f} - Buy: {v_list[0]:.2f} Sell: {v_list[1]:.2f} Hold: {v_list[2]:.2f} - Max_value {max_value:.2f} - maxDec {decList[v_list.index(max(v_list))]}")
                maxDict.update({S_t: maxValue})
            V_tIS_tI[t] = maxDict
            t -= 1
        pass
        
        self.V_tIS_tI = V_tIS_tI
        return V_tIS_tI

4.6 Evaluating Policies

4.6.1 Training/Tuning

# UPDATE PARAMETERS
params.update({'Algorithm': 'BackwardDP'})
params.update({'price_disc_list': [1.0, 0.5]})
params.update({'R__max': 1})
params

{'Algorithm': 'BackwardDP',
 'T': 192,
 'eta': 0.9,
 'R__max': 1,
 'R_0': 1,
 'seed': 189654913,
 'theta_sell_min': 10,
 'theta_sell_max': 100,
 'theta_buy_min': 10,
 'theta_buy_max': 100,
 'theta_inc': 1,
 'nPriceChangeInc': 50,
 'priceDiscSet': '1, 0.5',
 'run3D': False,
 'price_disc_list': [1.0, 0.5]}

# create a model and a policy
piNames = ['X__bellman_1hot']
SVarNames = ['R_t', 'p_t']
initialState = {
    'R_t': params['R_0'],
    'p_t': exogParams['hist_price'][0],
}
xVarNames = ['x_t_buy', 'x_t_hold', 'x_t_sell']
possibleDecisions = [
    {'x_t_buy': 1, 'x_t_hold': 0, 'x_t_sell': 0}, 
    {'x_t_buy': 0, 'x_t_hold': 0, 'x_t_sell': 1},
    {'x_t_buy': 0, 'x_t_hold': 1, 'x_t_sell': 0}]
M = EnergyStorageModel(
    SVarNames, 
    xVarNames, 
    initialState, 
    params, 
    exogParams,
    possibleDecisions
)
P = EnergyStoragePolicy(M, piNames)

%%time
#################################################################################
#BackwardDP
if params['Algorithm'] == 'BackwardDP':
    #Constructing the state space
    # make list of possible energy amount stored at a time
    discreteEnergy = np.array([0., 1.])

    # make list of prices with different increments
    minPrice = np.min(exogParams['hist_price'])
    maxPrice = np.max(exogParams['hist_price'])
            
    incValues = params['price_disc_list']
    # incValues = [1.0] #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! REMOVE
    cumCI_inc_I, best_bdp, bdpI_inc_I = \
        P.perform_bellman_search_1hot(
            params, discreteEnergy, minPrice, maxPrice, incValues)
# #########################################################################
# SUMMARY
# P.perform_bellman_search()
#     for inc in inc_values: #for each discretization
#         test_2D_bdp.bellman()|bellman_1hot()
#             while t != -1: #. for each t
#                 for S_t in self.possible_states: #. for each S_t
#                     for d in self.model.possible_decisions: #. for each x_t
#                         for w in self.discrete_price_changes: #. for each w
#             return V_tIS_tI
#     return cumCI_inc_I, best_bdp, bdpI_inc_I

Starting BackwardDP 2D for price discretization with inc=1.0
State dimension: 2. State space size: 218. Exogenous info size: 50
self.discretePrices=array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,
        11.,  12.,  13.,  14.,  15.,  16.,  17.,  18.,  19.,  20.,  21.,
        22.,  23.,  24.,  25.,  26.,  27.,  28.,  29.,  30.,  31.,  32.,
        33.,  34.,  35.,  36.,  37.,  38.,  39.,  40.,  41.,  42.,  43.,
        44.,  45.,  46.,  47.,  48.,  49.,  50.,  51.,  52.,  53.,  54.,
        55.,  56.,  57.,  58.,  59.,  60.,  61.,  62.,  63.,  64.,  65.,
        66.,  67.,  68.,  69.,  70.,  71.,  72.,  73.,  74.,  75.,  76.,
        77.,  78.,  79.,  80.,  81.,  82.,  83.,  84.,  85.,  86.,  87.,
        88.,  89.,  90.,  91.,  92.,  93.,  94.,  95.,  96.,  97.,  98.,
        99., 100., 101., 102., 103., 104., 105., 106., 107., 108.])
self.discreteEnergy=array([0., 1.])
Time_elapsed_2D_model=147.04 secs.
Starting policy evaluation for the actual sample path
Finishing inc 1.0 with cumC 249.66. Best inc so far 1.0. Best contribution 249.66

Starting BackwardDP 2D for price discretization with inc=0.5
State dimension: 2. State space size: 432. Exogenous info size: 50
self.discretePrices=array([  0. ,   0.5,   1. ,   1.5,   2. ,   2.5,   3. ,   3.5,   4. ,
         4.5,   5. ,   5.5,   6. ,   6.5,   7. ,   7.5,   8. ,   8.5,
         9. ,   9.5,  10. ,  10.5,  11. ,  11.5,  12. ,  12.5,  13. ,
        13.5,  14. ,  14.5,  15. ,  15.5,  16. ,  16.5,  17. ,  17.5,
        18. ,  18.5,  19. ,  19.5,  20. ,  20.5,  21. ,  21.5,  22. ,
        22.5,  23. ,  23.5,  24. ,  24.5,  25. ,  25.5,  26. ,  26.5,
        27. ,  27.5,  28. ,  28.5,  29. ,  29.5,  30. ,  30.5,  31. ,
        31.5,  32. ,  32.5,  33. ,  33.5,  34. ,  34.5,  35. ,  35.5,
        36. ,  36.5,  37. ,  37.5,  38. ,  38.5,  39. ,  39.5,  40. ,
        40.5,  41. ,  41.5,  42. ,  42.5,  43. ,  43.5,  44. ,  44.5,
        45. ,  45.5,  46. ,  46.5,  47. ,  47.5,  48. ,  48.5,  49. ,
        49.5,  50. ,  50.5,  51. ,  51.5,  52. ,  52.5,  53. ,  53.5,
        54. ,  54.5,  55. ,  55.5,  56. ,  56.5,  57. ,  57.5,  58. ,
        58.5,  59. ,  59.5,  60. ,  60.5,  61. ,  61.5,  62. ,  62.5,
        63. ,  63.5,  64. ,  64.5,  65. ,  65.5,  66. ,  66.5,  67. ,
        67.5,  68. ,  68.5,  69. ,  69.5,  70. ,  70.5,  71. ,  71.5,
        72. ,  72.5,  73. ,  73.5,  74. ,  74.5,  75. ,  75.5,  76. ,
        76.5,  77. ,  77.5,  78. ,  78.5,  79. ,  79.5,  80. ,  80.5,
        81. ,  81.5,  82. ,  82.5,  83. ,  83.5,  84. ,  84.5,  85. ,
        85.5,  86. ,  86.5,  87. ,  87.5,  88. ,  88.5,  89. ,  89.5,
        90. ,  90.5,  91. ,  91.5,  92. ,  92.5,  93. ,  93.5,  94. ,
        94.5,  95. ,  95.5,  96. ,  96.5,  97. ,  97.5,  98. ,  98.5,
        99. ,  99.5, 100. , 100.5, 101. , 101.5, 102. , 102.5, 103. ,
       103.5, 104. , 104.5, 105. , 105.5, 106. , 106.5, 107. , 107.5])
self.discreteEnergy=array([0., 1.])
Time_elapsed_2D_model=480.19 secs.
Starting policy evaluation for the actual sample path
Finishing inc 0.5 with cumC 278.75. Best inc so far 0.5. Best contribution 278.75
Finishing BellmanSearch in 0.95 secs
CPU times: user 10min 21s, sys: 1.17 s, total: 10min 23s
Wall time: 10min 29s

cumCI_inc_I

{1.0: 249.66233333333327, 0.5: 278.75477777777775}

best_bdp

<__main__.BDP at 0x7fb9a72a3df0>

4.6.2 Evaluation

params['R_0']

1

# create a model and a policy
piNames = ['X__bellman_1hot']
SNames = ['R_t', 'p_t']
initialState = {
    # 'R_t': params['R_0'],
    'R_t': 1,
    'p_t': exogParams['hist_price'][0],
}
xVarNames = ['x_t_buy', 'x_t_hold', 'x_t_sell']
possibleDecisions = M.possibleDecisions
M_evalu = EnergyStorageModel(
    SVarNames, 
    xVarNames, 
    initialState, 
    params, 
    exogParams,
    possibleDecisions
)
P_evalu = EnergyStoragePolicy(M_evalu, piNames)
params

{'Algorithm': 'BackwardDP',
 'T': 192,
 'eta': 0.9,
 'R__max': 1,
 'R_0': 1,
 'seed': 189654913,
 'theta_sell_min': 10,
 'theta_sell_max': 100,
 'theta_buy_min': 10,
 'theta_buy_max': 100,
 'theta_inc': 1,
 'nPriceChangeInc': 50,
 'priceDiscSet': '1, 0.5',
 'run3D': False,
 'price_disc_list': [1.0, 0.5]}

def run_policy_evalu_1hot(piInfo_evalu, piName_evalu, stop_time_evalu):
    model_copy = copy(M_evalu)
    nTrades = {'x_t_buy': 0, 'x_t_sell': 0, 'x_t_hold': 0}
    buyList = []
    sellList = []
    record = []
    for t in range(stop_time_evalu):
        x_t_info = getattr(P_evalu, piName_evalu)(t, model_copy.S_t, piInfo_evalu, stop_time_evalu)
        if t == stop_time_evalu - 1: #Last time period - sell
            x_t_info = {'x_t_buy': 0, 'x_t_hold': 0, 'x_t_sell': 1}   
        x_t = model_copy.build_decision_1hot(x_t_info, model_copy.S_t.R_t)
        # nTrades['x_t_buy'] += x_t.x_t_buy; #print(f'### {x_t.buy=}')
        # nTrades['x_t_sell'] += x_t.x_t_sell; #print(f'### {x_t.sell=}')
        # nTrades['x_t_hold'] += model_copy.S_t.R_t#.
        if x_t.x_t_buy > 0:
            buyList.append((t, model_copy.S_t.p_t))
        elif x_t.x_t_sell > 0:
            sellList.append((t, model_copy.S_t.p_t))
        res = model_copy.step_1hot(t, x_t) # step the model forward one iteration
        record.append([res[0].R_t, res[0].p_t, res[1], res[2].x_t_buy, res[2].x_t_hold, res[2].x_t_sell])
    cumC = model_copy.cumC
    # print("Energy traded - Sell: {:.2f} - Buy: {:.2f} - Hold % : {:.2f}".format(nTrades['sell'], nTrades['buy'], nTrades['hold']/model_copy.init_args['T']))
    # print(f"Energy traded - Sell: {nTrades['x_t_sell']:.2f} - Buy: {nTrades['x_t_buy']:.2f} - Hold % : {nTrades['x_t_hold']/model_copy.init_args['T']:.2f}") #.
    # print("Sell times and prices ")
    # for i in range(len(sell_list)):
    #     print(f"t = {sell_list[i][0]} and price = {sell_list[i][1]:.2f}") #.
    # print("Buy times and prices ")
    # for i in range(len(buy_list)):
    #     print(f"t = {buy_list[i][0]} and price = {buy_list[i][1]:.2f}") #.
    return cumC, record

# EVALUATION
piName_evalu = 'X__bellman_1hot'
stop_time_evalu = 180

4.6.2.1 Optimal policy

piName_evalu

'X__bellman_1hot'

# inc = 1.0
inc = 0.5 #increment for optimal
discreteEnergy = np.array([0., 1.])
minPrice = np.min(exogParams['hist_price'])
maxPrice = np.max(exogParams['hist_price'])
discrete_prices = np.arange(minPrice, maxPrice + inc, inc)
bdp_2D_evalu = best_bdp
cumC, record = run_policy_evalu_1hot(bdp_2D_evalu, piName_evalu, stop_time_evalu=params['T'])
labels = ["R_t", "p_t", "cumC", "x_t_buy", "x_t_hold", "x_t_sell"]
print(f'{bdp_2D_evalu=}')
df = pd.DataFrame.from_records(data=record, columns=labels); df[:20]

bdp_2D_evalu=<__main__.BDP object at 0x7fb9a72a3df0>

	R_t	p_t	cumC	x_t_buy	x_t_hold	x_t_sell
0	1.0000	18.9100	0.0000	0.0000	0	0.0000
1	1.0000	14.2700	0.0000	0.0000	0	0.0000
2	1.0000	9.6700	0.0000	0.0000	0	0.0000
3	1.0000	17.6800	0.0000	0.0000	0	0.0000
4	1.0000	11.8800	0.0000	0.0000	0	0.0000
5	1.0000	2.0800	0.0000	0.0000	0	0.0000
6	1.0000	0.0000	0.0000	0.0000	0	0.0000
7	1.0000	18.0000	0.0000	0.0000	0	0.0000
8	1.0000	14.4800	0.0000	0.0000	0	0.0000
9	1.0000	17.8800	0.0000	0.0000	0	0.0000
10	1.0000	24.6400	0.0000	0.0000	0	0.0000
11	1.0000	19.5700	0.0000	0.0000	0	0.0000
12	1.0000	20.2600	0.0000	0.0000	0	0.0000
13	1.0000	22.2100	0.0000	0.0000	0	0.0000
14	1.0000	18.5200	0.0000	0.0000	0	0.0000
15	1.0000	18.7800	0.0000	0.0000	0	0.0000
16	1.0000	23.0900	0.0000	0.0000	0	0.0000
17	1.0000	30.6900	0.0000	0.0000	0	0.0000
18	1.0000	22.8800	0.0000	0.0000	0	0.0000
19	1.0000	24.6000	0.0000	0.0000	0	0.0000

df.tail()

	R_t	p_t	cumC	x_t_buy	x_t_hold	x_t_sell
187	0.0000	53.2500	278.7548	0.0000	0	0.0000
188	0.0000	48.2300	278.7548	0.0000	0	0.0000
189	0.0000	33.4200	278.7548	0.0000	0	0.0000
190	0.0000	32.3800	278.7548	0.0000	0	0.0000
191	0.0000	33.1900	278.7548	0.0000	0	0.0000

4.6.2.2 Non-optimal policy

piName_evalu = 'X__bellman_1hot'
bdp_2D_evalu = bdpI_inc_I[1.0] #non-optimal
cumC, record = run_policy_evalu_1hot(bdp_2D_evalu, piName_evalu, stop_time_evalu=params['T'])
labels = ["R_t", "p_t", "cumC", "x_t_buy", "x_t_hold", "x_t_sell"]
print(f'{bdp_2D_evalu=}')
df_non = pd.DataFrame.from_records(data=record, columns=labels); df_non[:20]

bdp_2D_evalu=<__main__.BDP object at 0x7fb9a732c5e0>

	R_t	p_t	cumC	x_t_buy	x_t_hold	x_t_sell
0	1.0000	18.9100	0.0000	0.0000	0	0.0000
1	1.0000	14.2700	0.0000	0.0000	0	0.0000
2	1.0000	9.6700	0.0000	0.0000	0	0.0000
3	1.0000	17.6800	0.0000	0.0000	0	0.0000
4	1.0000	11.8800	0.0000	0.0000	0	0.0000
5	1.0000	2.0800	0.0000	0.0000	0	0.0000
6	1.0000	0.0000	0.0000	0.0000	0	0.0000
7	1.0000	18.0000	0.0000	0.0000	0	0.0000
8	1.0000	14.4800	0.0000	0.0000	0	0.0000
9	1.0000	17.8800	0.0000	0.0000	0	0.0000
10	1.0000	24.6400	0.0000	0.0000	0	0.0000
11	1.0000	19.5700	0.0000	0.0000	0	0.0000
12	1.0000	20.2600	0.0000	0.0000	0	0.0000
13	1.0000	22.2100	0.0000	0.0000	0	0.0000
14	1.0000	18.5200	0.0000	0.0000	0	0.0000
15	1.0000	18.7800	0.0000	0.0000	0	0.0000
16	1.0000	23.0900	0.0000	0.0000	0	0.0000
17	1.0000	30.6900	0.0000	0.0000	0	0.0000
18	1.0000	22.8800	0.0000	0.0000	0	0.0000
19	1.0000	24.6000	0.0000	0.0000	0	0.0000

4.6.2.3 Visualization

from certifi.core import where
def plot_output_1hot(df1, df2):
  legendlabels =     [r'$x_{t,buy}$', r'$x_{t,hold}$', r'$x_{t,sell}$']
  n_charts = 4
  ylabelsize = 16
  mpl.rcParams['lines.linewidth'] = 1.2
  fig, axs = plt.subplots(n_charts, sharex=True)
  # fig.set_figwidth(50); fig.set_figheight(20)
  fig.set_figwidth(13); fig.set_figheight(9)
  fig.suptitle('PERFORMANCE OF OPTIMIZED Bellman POLICY', fontsize=20)
  i = 0 #buy/hold/sell
  axs[i].set_ylim(auto=True); axs[i].spines['top'].set_visible(False); axs[i].spines['right'].set_visible(True); axs[i].spines['bottom'].set_visible(False)
  axs[i].step(df1['x_t_buy'], 'r')
  axs[i].step(df1['x_t_hold'], 'g')
  axs[i].step(-df1['x_t_sell'], 'b')
  axs[i].axhline(y=0, color='k', linestyle=':')  
  axs[i].set_ylabel('$x_{t,buy}$\n$x_{t,hold}$\n$x_{t,sell}$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);
  i = 1 #R_t
  axs[i].set_ylim(auto=True); axs[i].spines['top'].set_visible(False); axs[i].spines['right'].set_visible(True); axs[i].spines['bottom'].set_visible(False)
  axs[i].step(df1['R_t'], 'c')
  axs[i].axhline(y=0, color='k', linestyle=':')
  axs[i].set_ylabel('$R_t$'+'\n'+'$\mathrm{[MWh]}$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);
  i = 2 #p_t
  axs[i].set_ylim(auto=True); axs[i].spines['top'].set_visible(False); axs[i].spines['right'].set_visible(True); axs[i].spines['bottom'].set_visible(False)
  axs[i].step(df1['p_t'], 'm')
  # axs[i].axhline(y=theta_evalu[0], color='m', linestyle=':')
  # axs[i].text(0, theta_evalu[0], r'$\theta^{buy}$', size=16)
  # axs[i].axhline(y=theta_evalu[1], color='m', linestyle='--')
  # axs[i].text(0, theta_evalu[1], r'$\theta^{sell}$', size=16)
  axs[i].set_ylabel('$p_t$'+'\n'+'$\mathrm{[\$/MWh]}$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);
  i = 3 #cumC
  axs[i].set_ylim(auto=True); axs[i].spines['top'].set_visible(False); axs[i].spines['right'].set_visible(True); axs[i].spines['bottom'].set_visible(False)
  axs[i].step(df1['cumC'], 'k')
  axs[i].step(df2['cumC'], 'k:')
  axs[i].text(0, 100, r'$Optimal\ (solid)$', size=16)
  axs[i].text(0, 50, r'$Non-Optimal\ (dotted)$', size=16)
  axs[i].set_ylabel('$\mathrm{cumC}$'+'\n'+'$\mathrm{(Profit)}$'+'\n'+''+'$\mathrm{[\$]}$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);
  axs[i].set_xlabel('$t\ \mathrm{[hours]}$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);
  fig.legend(labels=legendlabels, loc='center', fontsize=18)
plot_output_1hot(df, df_non)