Dealership Inventory Management using the Powell Unified Framework (Part 1)

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. 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.

In order to make a strong mapping between the code in this notebook and the mathematics in the Powell Unified Framework (PUF), we follow the following convention for naming Python identifier names:

• 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', SVarNames) 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', xVarNames) 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

Inventory management is a critical component of any business, whether it be a small retail store or a multinational corporation. At its core, inventory management is the process of tracking and controlling a company’s inventory, from raw materials to finished products. Proper inventory management is important for several reasons.

First and foremost, inventory management helps businesses avoid stock overages and underages (overstocks and stockouts). By tracking inventory levels and forecasting demand, businesses can ensure that they always have the right amount of product on hand to meet customer needs without overbuying and tying up capital in excess inventory. This helps businesses maintain a healthy cash flow and avoid costly stockouts that can result in lost sales and dissatisfied customers.

In addition, effective inventory management can help businesses streamline their operations and improve their overall efficiency. By reducing excess inventory and optimizing order quantities and lead times, businesses can minimize waste and improve their supply chain management. This can lead to cost savings, improved profitability, and increased customer satisfaction.

Finally, inventory management is critical for businesses that need to comply with regulatory requirements, such as those in the pharmaceutical or food industries. Proper inventory tracking and documentation can help businesses meet these requirements and avoid costly fines and penalties.

Overall, inventory management is an essential function for any business that wants to operate efficiently, meet customer demand, and maximize profitability. Effective inventory management requires careful planning, accurate data, and the right tools and processes to ensure that businesses always have the right amount of product on hand, at the right time, and at the right cost.

In this project the client had a need to be convinced of the benefits of formal optimized sequential decision making. This was provided in the form of a series of POCs.

2 DATA UNDERSTANDING

Next, we will look at how we will simulate the data for this problem.

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

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

We will simulate the inventory demand, \(D\), by: \[ D_{t+1} = \theta_0 D_t + \theta_1 D_{t-1} + \theta_2 D_{t-2} + \epsilon_{pois(5)} \] where \[ \begin{aligned} \theta_0 &= 0.14 \\ \theta_1 &= 0.31 \\ \theta_2 &= 0.25 \\ \end{aligned} \]

class DemandSimulator():
    def __init__(self, th_0, th_1, th_2, D_t, D_t_1, D_t_2):
        self.th_0 = th_0
        self.th_1 = th_1
        self.th_2 = th_2
        self.D_t = D_t
        self.D_t_1 = D_t_1
        self.D_t_2 = D_t_2

    def simulate(self):
        D_tt1 = self.th_0*self.D_t + self.th_1*self.D_t_1 + self.th_2*self.D_t_2 + np.random.poisson(5)
        d = int(max(0, D_tt1))
        self.D_t_2 = self.D_t_1
        self.D_t_1 = self.D_t
        self.D_t = D_tt1
        return d

dem_sim = DemandSimulator(.14,.31,.25, 5,5,5)
DemandData = []
for i in range(100):
  d = dem_sim.simulate()
  DemandData.append([d])

labels = ['demand']
df = pd.DataFrame.from_records(data=DemandData, columns=labels); df[:10]

	demand
0	6
1	6
2	9
3	10
4	10
5	8
6	14
7	12
8	12
9	13

def plot_output(df1, df2):
  n_charts = 2
  ylabelsize = 16
  mpl.rcParams['lines.linewidth'] = 1.2
  fig, axs = plt.subplots(n_charts, sharex=True)
  fig.set_figwidth(13); fig.set_figheight(9)
  fig.suptitle('Demand Simulation', fontsize=20)

  i = 0 #Demand
  axs[i].set_title('Demanded units')
  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['demand'], 'k')

  i = 1 #parameters
  axs[i].set_title('Parameter values')
  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].axhline(y=dem_sim.th_0, color='r', linestyle=':')
  axs[i].axhline(y=dem_sim.th_1, color='g', linestyle=':')
  axs[i].axhline(y=dem_sim.th_2, color='b', linestyle=':')
plot_output(df, None)

seed = 189654913
file = 'Parameters.xlsx'

# NOTE:
# R__max: maximum number of inventory units
# R_0:    initial number of inventory units
parDf = pd.read_excel(f'{base_dir}/{file}', sheet_name='ParamsModel', index_col=0); print(f'{parDf}')
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  GridSearch
T                 195
eta            0.9000
R__max              1
R_0                 1
params={'Algorithm': 'GridSearch', '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.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': 'GridSearch', '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}")

pprint(f"{params=}")

("params={'Algorithm': 'GridSearch', '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}")

3 DATA PREPARATION

We will use the data provided by the simulator directly. There is no need to perform additional data preparation.

4 MODELING

4.1 Narrative

This first project in the Inventory Series tackles the simplest of inventory problems.

We have the following setting: Mr. Optimal is car lot inventory manager for the largest dealership in a big city. For now, he will only be responsible to manage the inventory level of a single car model, that of the Hyunday Elantra. He has a maximum number of lot spaces assigned to him (which is 57). He has a choice to strive to always keep these spaces occupied by new cars. This way he is unlikely to run out of stock and lose a sale due to that. However, capital is tied up by the unsold inventory in his lot space.

At the other extreme, he may choose to work on a just-in-time principle: Each time a potential customer expresses interest in this model, the customer will have to wait until he obtains a new car from the supplier. Of course, he will likely lose the sale, but the upside is that no capital is tied up in his inventory.

It seems intuitive that the optimal level of inventory will be somewhere between these extremes. The challenge is to find that optimal level. For now, we will assume that the buy and sell prices will remain constant. The only random variable will be the demand for this specific model. Another assumption is that ordered inventory will arrive immediately. In addition, unsatisfied demands are lost, i.e. there will be no ability to backlog unsatisfied demands.

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 many new cars to order. The only source of uncertainty is the level of demand for this model.

4.3 Mathematical Model | SUM Design

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

class InventoryStorageModel():
    def __init__(
        self, SVarNames, xVarNames, S_0, params, exogParams, possibleDecisions,
        p__buy, p__sell, 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\)
  • the inventory on hand at time \(t\) before we make a new ordering decision, and before we have satisfied any demands arising in time interval \(t\)
  • measured in inventory units
• \(D_t\)
  • the demand
  • measured in inventory units

The state is:

\(S_t = (R_t,D_t)\)

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

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

where

SVarNames = ['R_t', 'D_t']

4.3.2 Decision variables

The decision variables represent what we control.

• \(x_t = (x_{t,ELA})\)
  • number of Elantras ordered (\(x_t\ge0\)) where \(x_t\) is a positive integer
  • we represent this scalar by a vector with a single component to allow for easy expansion later in the series
• Constraints
  • \(x_{t,ELA} \le (R^{max} - R_t)\)
  • \(R^{max}\) is the number of lot units (i.e. parking spaces) assigned to Mr. Optimal
• Decisions are made with a policy (TBD below):
  • \(X^{\pi}(S_t)\)

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

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

where

xVarNames = ['x_t_ELA']

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\).

We assume that any unsatisfied demand is lost. Additionally, we assume that the demand in each time period is revealed, so that we have:

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

The exogenous information is obtained by a call to

DemandSimulator.simulate_demand(...)

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

def W_fn(self, t):
  W_tt1 = dem_sim.simulate()
  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,D_t)\), we have the equations:

\[ \begin{aligned} R_{t+1} &= R_t - \mathrm{min}\{R_t,D_t\}+x_t \quad (Eq. 1) \\ D_{t+1} &= \hat{D}_{t+1} \quad\quad\quad\quad\quad\quad\quad\quad (Eq. 2) \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 InventoryStorageModel():

def S__M_fn(self, t, x_t): #.
    D_tt1 = self.W_fn(t)
    R_tt1 = max( 0, self.S_t.R_t - min(self.S_t.R_t, self.S_t.D_t) + x_t.x_t_ELA ) #max to keep >0
    
    if len(self.SVarNames) == 2:
        S_tt1 = self.build_state({'R_t': R_tt1, 'D_t': D_tt1})
    elif len(self.SVarNames) == 3:
        S_tt1 = self.build_state({'R_t': R_tt1, 'D_t': D_tt1, 'D_t_1': self.S_t.D_t})
    return S_tt1

Note that x_t.x_t_ELA is currently the only component of \(x_t\). The ELA stands for Elantra.

4.3.5 Objective function

The objective function captures the performance metrics of the solution to the problem.

We can write the state-based reward (also called contribution) based on what we will receive between \(t-1\) and \(t\) (i.e. looking backward relative to \((S_t,x_t)\)):

\[ \begin{aligned} C(S_t,x_t) = p^{sell}\min\{R_t,D_t\} - p^{buy}x_t \end{aligned} \] This is a deterministic expression.

Alternatively, we can write the state-based reward based on what we will receive between \(t\) and \(t+1\) (i.e. looking forward relative to \((S_t,x_t)\)):

\[ \begin{aligned} C(S_t,x_t,\hat{D}_{t+1}) &= p^{sell}\min\{R_{t+1},D_{t+1}\} - p^{buy}x_t \\ &= p^{sell}\mathrm{min}\{(R_t-\mathrm{min}\{R_t,D_t\}+x_t), \hat{D}_{t+1}\} - p^{buy}x_t \quad \mathrm{(RLSO-Eq8.5)} \end{aligned} \] because, from (Eq. 1) and (Eq. 2) above: \[ \begin{aligned} R_{t+1} &= R_t - \mathrm{min}\{R_t,D_t\}+x_t \quad (Eq. 1) \\ D_{t+1} &= \hat{D}_{t+1} \quad\quad\quad\quad\quad\quad\quad\quad (Eq. 2) \end{aligned} \] This is a stochastic expression due to the dependence on the random variable \(\hat{D}_{t+1}\). It is random because it comes from a stochastic process but it is also in the future.

This second form leads to the objective function:

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

where \[ C(S_t,x_t,W_{t+1}) = p^{sell}\mathrm{min}\{(R_t-\mathrm{min}\{R_t,D_t\}+x_t), \hat{D}_{t+1}\} - p^{buy}x_t \]

Note that \(x_t\) is a single-component vector. In the present case \(x_t = (x_{t,ELA})\)

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

def C_fn(self, x_t):
    Dhat_tt1 = dem_sim.simulate()
    C_t = \
    self.p__sell*min((self.S_t.R_t - min(self.S_t.R_t, self.S_t.D_t) + x_t.x_t_ELA), Dhat_tt1) \
      - self.p__buy*x_t.x_t_ELA
    return C_t

4.3.6 Implementation of SUM Model

Here is the complete implementation of the InventoryStorageModel class:

class InventoryStorageModel():
    def __init__(
        self, SVarNames, xVarNames, S_0, params, exogParams, possibleDecisions,
        p__buy, p__sell, 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.p__buy = p__buy
        self.p__sell = p__sell
        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

    def reset(self):
        self.cumC = 0.0
        self.S_t = self.build_state(self.S_0)

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

    def build_decision(self, info):
        return self.Decision(*[info[xn] for xn in self.xVarNames])

    def W_fn(self, t): #.
        W_tt1 = dem_sim.simulate()
        return W_tt1

    def S__M_fn(self, t, x_t): #.
        D_tt1 = self.W_fn(t) 
        R_tt1 = max( 0, self.S_t.R_t - min(self.S_t.R_t, self.S_t.D_t) + x_t.x_t_ELA ) #max to keep >0
        if len(self.SVarNames) == 2:
            S_tt1 = self.build_state({'R_t': R_tt1, 'D_t': D_tt1})
        elif len(self.SVarNames) == 3:
            S_tt1 = self.build_state({'R_t': R_tt1, 'D_t': D_tt1, 'D_t_1': self.S_t.D_t})
        return S_tt1

    # based on what we will receive between t and t+1 (i.e. looking *forward* relative to (S_t,x_t) #.
    # RLSO-Eq8.5
    def C_fn(self, x_t):
        Dhat_tt1 = dem_sim.simulate()
        C_t = \
        self.p__sell*min((self.S_t.R_t - min(self.S_t.R_t, self.S_t.D_t) + x_t.x_t_ELA), Dhat_tt1) \
          - self.p__buy*x_t.x_t_ELA
        return C_t

    def step(self, t, x_t):
        self.cumC += self.C_fn(x_t)
        self.S_t = self.S__M_fn(t, x_t)
        return (self.S_t, self.cumC, x_t) #. for plotting

4.4 Uncertainty Model

As described in section 2, we will simulate the inventory demand, \(D\), by: \[ D_{t+1} = \theta_0 D_t + \theta_1 D_{t-1} + \theta_2 D_{t-2} + \epsilon_{pois(5)} \]

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 only use one approach: - A simple buy below parameterized policy (from the PFA class)

The buy below policy is implemented by the following method in class InventoryStoragePolicy():

def X__BuyBelow(self, t, S_t, theta, T):
    info = {
        'x_t_ELA': 0, #number of Elantras ordered
    }
    if t >= T:
        print(f"ERROR: t={t} should not reach or exceed the max steps ({T})")
        return self.model.build_decision(info)
    theta__buy = theta[0] # theta__buy
    if S_t.R_t <= theta__buy: # BUY if R_t <= theta__buy
        info['x_t_ELA'] = self.model.initArgs['R__max'] - S_t.R_t
    return self.model.build_decision(info)

4.5.1 Implementation of Policy Design

The InventoryStoragePolicy() class implements the policy design.

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

    def X__BuyBelow(self, t, S_t, theta, T):
        info = {
            'x_t_ELA': 0, #number of Elantras ordered
        }
        if t >= T:
            print(f"ERROR: t={t} should not reach or exceed the max steps ({T})")
            return self.model.build_decision(info)
        theta__buy = theta[0] # theta__buy
        if S_t.R_t <= theta__buy: # BUY if R_t <= theta__buy
            info['x_t_ELA'] = self.model.initArgs['R__max'] - S_t.R_t
        return self.model.build_decision(info)

    def run_policy(self, piInfo, piName, params):
        model_copy = copy(self.model)
        T = params['T']
        for t in range(T):
            x_t = getattr(self, piName)(t, model_copy.S_t, piInfo, T) # piInfo is theta value
            _, _, _ = model_copy.step(t, x_t)
        cumC = model_copy.cumC        
        return cumC

    def perform_grid_search(self, params, thetaValues):
        tS = time.time()
        cumCI_theta_I = {}
        bestTheta = None
        for theta in thetaValues:
            print(f'\n=== {theta=} ===')
            cumC = self.run_policy(theta, "X__BuyBelow", params)
            cumCI_theta_I[theta] = cumC
            best_theta = max(cumCI_theta_I, key=cumCI_theta_I.get)
            print(f"Finishing theta {theta} with cumC {cumC:.2f}. Best theta so far {best_theta}. Best cumC {cumCI_theta_I[best_theta]:.2f}")
        print(f"Finishing GridSearch in {time.time() - tS:.2f} secs")
        print(f"Best theta: {best_theta}. Best cumC: {cumCI_theta_I[best_theta]:.2f}")
        return cumCI_theta_I

    def run_policy_sample_paths(self, theta, piName, params):
        FhatIomega__lI = []
        for l in range(1, params['L'] + 1):
            model_copy = copy(self.model)
            record_l = [piName, theta, l]
            T = params['T']
            for t in range(T):
                x_t = getattr(self, piName)(t, model_copy.S_t, theta, T)
                _, _, _ = model_copy.step(t, x_t)
            FhatIomega__lI.append(model_copy.cumC) # just above (SDAM-eq2.9); Fhat for this sample-path is in model_copy.cumC
        return FhatIomega__lI

    def perform_grid_search_sample_paths(self, params, thetaValues):
        tS = time.time()
        Fhat_mean = None
        Fhat_var = None
        Fhat__meanI_th_I = {}
        Fhat__stdvI_th_I = {}
        Fhat_mean = None
        for theta in thetaValues:
            print(f'\n=== {theta=} ===')
            FhatIomega__lI = self.run_policy_sample_paths(theta, "X__BuyBelow", params)
            Fhat_mean = np.array(FhatIomega__lI).mean() #. (SDAM-eq2.9); call Fbar in future
            Fhat_var = np.sum(np.square(np.array(FhatIomega__lI) - Fhat_mean))/(params['L'] - 1)
            Fhat__meanI_th_I[theta] = Fhat_mean
            Fhat__stdvI_th_I[theta]= np.sqrt(Fhat_var/params['L'])
            best_theta = max(Fhat__meanI_th_I, key=Fhat__meanI_th_I.get)
            print(f"Finishing theta {theta} with cumC {Fhat__meanI_th_I[best_theta]:.2f}. Best theta so far {best_theta}. Best cumC {Fhat__meanI_th_I[best_theta]:.2f}")
        print(f"Finishing GridSearch in {time.time() - tS:.2f} secs")
        print(f"Best theta: {best_theta}. Best cumC: {Fhat__meanI_th_I[best_theta]:.2f}")
        return Fhat__meanI_th_I, Fhat__stdvI_th_I

    def grid_search_theta_values(self, params):
        theta_buy_values = np.arange(params['theta_buy_min'], params['theta_buy_max'], params['theta_inc'])
        # theta_sell_values = np.arange(params['theta_sell_min'], params['theta_sell_max'], params['theta_inc'])

        # theta_values = [(x, y) for x in theta_buy_values for y in theta_sell_values]        
        theta_values = [(x, None) for x in theta_buy_values]
        
        # return theta_values, theta_buy_values, theta_sell_values
        return theta_values, theta_buy_values

    def plot_heat_map(self, contribution_dict, theta_buy_values, theta_sell_values):
        contribution_values = [contribution_dict[(theta_buy,theta_sell)]  for theta_sell in theta_sell_values for theta_buy in theta_buy_values]
        contributions = np.array(contribution_values)
        increment_count = len(theta_buy_values)
        contributions = np.reshape(contributions, (-1, increment_count))

        fig, ax = plt.subplots()
        im = ax.imshow(contributions, 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(theta_buy_values),5))
        ax.set_yticks(np.arange(0,len(theta_sell_values),5))
        # ... and label them with the respective list entries
        ax.set_xticklabels(theta_buy_values[::5])
        ax.set_yticklabels(theta_sell_values[::5])
        # rotate the tick labels and set their alignment.
        #plt.setp(ax.get_xticklabels(), rotation=45, ha="right",rotation_mode="anchor")
        ax.set_title("Heatmap of contribution values across different values of theta")

        ax.set_ylabel('Theta sell high values') 
        ax.set_xlabel('Theta buy low values')

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

    def plot_cumC_chart(self, contribution_dict, theta_buy_values): # 1 theta
        plt.plot(theta_buy_values, contribution_dict.values())

    def plot_cumC_chart_sample_paths(self, Fhat__meanI_th_I, Fhat__stdvI_th_I, theta_buy_values, params): # 1 theta
        mpl.rcParams['lines.linewidth'] = 1.2
        ylabelsize = 16
        fig, axs = plt.subplots(1, 2)
        fig.set_figwidth(25); fig.set_figheight(8)
        fig.suptitle(f"Inventory Model: Cumulative Reward\n (L = {params['L']}, T = {params['T']})", fontsize=20)

        FhatMeans = [v for k,v in Fhat__meanI_th_I.items()] #.
        # axs[0].plot(theta_buy_values, Fhat__meanI_th_I, 'ro-')
        axs[0].plot(theta_buy_values, FhatMeans, 'ro-')#.
        axs[0].set_ylabel(r'$\bar{F}$'+'\n$\hat{F}^{mean}$'+'\n$[\$]$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);
        axs[0].set_xlabel(r'$\theta$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);

        FhatStdvs = [v for k,v in Fhat__stdvI_th_I.items()] #.
        # axs[1].plot(theta_buy_values, Fhat__stdvI_th_I, 'bo-')
        axs[1].plot(theta_buy_values, FhatStdvs, 'bo-')#.
        axs[1].set_ylabel('$\hat{F}^{stdv}$'+'\n$[\$]$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);
        axs[1].set_xlabel(r'$\theta$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);

        plt.show()

4.6 Policy Evaluation

4.6.1 Training/Tuning

4.6.1.1 With a few sample-paths

# UPDATE PARAMETERS
params.update({'Algorithm': 'GridSearch'}); pprint(f'{params=}')
params.update({'R__max': 57})
params.update({'R_0': 0})
params.update({'eta': None})

params.update({'theta_buy_min': 10})
params.update({'theta_buy_max': 50})
params.update({'theta_sell_min': None})
params.update({'theta_sell_max': None})

params.update({'L': 2}) #number of sample-paths
params.update({'T': 100_000}) #number of transitions/steps in each sample-path
params

("params={'Algorithm': 'GridSearch', '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}")

{'Algorithm': 'GridSearch',
 'T': 100000,
 'eta': None,
 'R__max': 57,
 'R_0': 0,
 'seed': 189654913,
 'theta_sell_min': None,
 'theta_sell_max': None,
 'theta_buy_min': 10,
 'theta_buy_max': 50,
 'theta_inc': 1,
 'L': 2}

# create a model and a policy
piNames = ['X__BuyBelow']
SVarNames = ['R_t', 'D_t']
initialState = {
    'R_t': params['R_0'],
    'D_t': dem_sim.simulate(),
}
xVarNames = ['x_t_ELA']
possibleDecisions = None
p__buy = 100
p__sell = 150
exogParams = {} # we use simulation
M = InventoryStorageModel(
    SVarNames, 
    xVarNames, 
    initialState, 
    params, 
    exogParams,
    possibleDecisions,
    p__buy,
    p__sell
)
P = InventoryStoragePolicy(M, piNames)

%%time
##########################################################################
#GridSearch #. SDAM-9.4.1
if params['Algorithm'] == 'GridSearch':
    # obtain the theta values to carry out a full grid search
    grid_search_theta_values = P.grid_search_theta_values(params)
    print(grid_search_theta_values)
    
    # cumCI_theta_I = \
    #   P.perform_grid_search(params, grid_search_theta_values[0])
    # P.plot_cumC_chart( # 1 theta
    #     cumCI_theta_I, #cumC values in here
    #     grid_search_theta_values[1]) #theta values    

    Fhat__meanI_th_I, Fhat__stdvI_th_I = \
      P.perform_grid_search_sample_paths(params, grid_search_theta_values[0])
    # P.plot_cumC_chart_sample_paths( # 1 theta
    #     Fhat__meanI_th_I, 
    #     Fhat__stdvI_th_I, 
    #     grid_search_theta_values[1], params)
##################################################################################

([(10, None), (11, None), (12, None), (13, None), (14, None), (15, None), (16, None), (17, None), (18, None), (19, None), (20, None), (21, None), (22, None), (23, None), (24, None), (25, None), (26, None), (27, None), (28, None), (29, None), (30, None), (31, None), (32, None), (33, None), (34, None), (35, None), (36, None), (37, None), (38, None), (39, None), (40, None), (41, None), (42, None), (43, None), (44, None), (45, None), (46, None), (47, None), (48, None), (49, None)], array([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]))

=== theta=(10, None) ===
Finishing theta (10, None) with cumC 64286200.00. Best theta so far (10, None). Best cumC 64286200.00

=== theta=(11, None) ===
Finishing theta (11, None) with cumC 64871450.00. Best theta so far (11, None). Best cumC 64871450.00

=== theta=(12, None) ===
Finishing theta (12, None) with cumC 65702325.00. Best theta so far (12, None). Best cumC 65702325.00

=== theta=(13, None) ===
Finishing theta (13, None) with cumC 66559125.00. Best theta so far (13, None). Best cumC 66559125.00

=== theta=(14, None) ===
Finishing theta (14, None) with cumC 67560675.00. Best theta so far (14, None). Best cumC 67560675.00

=== theta=(15, None) ===
Finishing theta (15, None) with cumC 68566925.00. Best theta so far (15, None). Best cumC 68566925.00

=== theta=(16, None) ===
Finishing theta (16, None) with cumC 69824550.00. Best theta so far (16, None). Best cumC 69824550.00

=== theta=(17, None) ===
Finishing theta (17, None) with cumC 70624400.00. Best theta so far (17, None). Best cumC 70624400.00

=== theta=(18, None) ===
Finishing theta (18, None) with cumC 71601000.00. Best theta so far (18, None). Best cumC 71601000.00

=== theta=(19, None) ===
Finishing theta (19, None) with cumC 72235400.00. Best theta so far (19, None). Best cumC 72235400.00

=== theta=(20, None) ===
Finishing theta (20, None) with cumC 73062800.00. Best theta so far (20, None). Best cumC 73062800.00

=== theta=(21, None) ===
Finishing theta (21, None) with cumC 73854200.00. Best theta so far (21, None). Best cumC 73854200.00

=== theta=(22, None) ===
Finishing theta (22, None) with cumC 74185475.00. Best theta so far (22, None). Best cumC 74185475.00

=== theta=(23, None) ===
Finishing theta (23, None) with cumC 74981250.00. Best theta so far (23, None). Best cumC 74981250.00

=== theta=(24, None) ===
Finishing theta (24, None) with cumC 75915675.00. Best theta so far (24, None). Best cumC 75915675.00

=== theta=(25, None) ===
Finishing theta (25, None) with cumC 77126775.00. Best theta so far (25, None). Best cumC 77126775.00

=== theta=(26, None) ===
Finishing theta (26, None) with cumC 78204950.00. Best theta so far (26, None). Best cumC 78204950.00

=== theta=(27, None) ===
Finishing theta (27, None) with cumC 79429625.00. Best theta so far (27, None). Best cumC 79429625.00

=== theta=(28, None) ===
Finishing theta (28, None) with cumC 80083750.00. Best theta so far (28, None). Best cumC 80083750.00

=== theta=(29, None) ===
Finishing theta (29, None) with cumC 80118925.00. Best theta so far (29, None). Best cumC 80118925.00

=== theta=(30, None) ===
Finishing theta (30, None) with cumC 80314225.00. Best theta so far (30, None). Best cumC 80314225.00

=== theta=(31, None) ===
Finishing theta (31, None) with cumC 80433650.00. Best theta so far (31, None). Best cumC 80433650.00

=== theta=(32, None) ===
Finishing theta (32, None) with cumC 80433650.00. Best theta so far (31, None). Best cumC 80433650.00

=== theta=(33, None) ===
Finishing theta (33, None) with cumC 80433650.00. Best theta so far (31, None). Best cumC 80433650.00

=== theta=(34, None) ===
Finishing theta (34, None) with cumC 80559325.00. Best theta so far (34, None). Best cumC 80559325.00

=== theta=(35, None) ===
Finishing theta (35, None) with cumC 80559325.00. Best theta so far (34, None). Best cumC 80559325.00

=== theta=(36, None) ===
Finishing theta (36, None) with cumC 80674400.00. Best theta so far (36, None). Best cumC 80674400.00

=== theta=(37, None) ===
Finishing theta (37, None) with cumC 80674400.00. Best theta so far (36, None). Best cumC 80674400.00

=== theta=(38, None) ===
Finishing theta (38, None) with cumC 80823275.00. Best theta so far (38, None). Best cumC 80823275.00

=== theta=(39, None) ===
Finishing theta (39, None) with cumC 80823275.00. Best theta so far (38, None). Best cumC 80823275.00

=== theta=(40, None) ===
Finishing theta (40, None) with cumC 80823275.00. Best theta so far (38, None). Best cumC 80823275.00

=== theta=(41, None) ===
Finishing theta (41, None) with cumC 80883000.00. Best theta so far (41, None). Best cumC 80883000.00

=== theta=(42, None) ===
Finishing theta (42, None) with cumC 80883000.00. Best theta so far (41, None). Best cumC 80883000.00

=== theta=(43, None) ===
Finishing theta (43, None) with cumC 80899400.00. Best theta so far (43, None). Best cumC 80899400.00

=== theta=(44, None) ===
Finishing theta (44, None) with cumC 80899400.00. Best theta so far (43, None). Best cumC 80899400.00

=== theta=(45, None) ===
Finishing theta (45, None) with cumC 80899400.00. Best theta so far (43, None). Best cumC 80899400.00

=== theta=(46, None) ===
Finishing theta (46, None) with cumC 80899400.00. Best theta so far (43, None). Best cumC 80899400.00

=== theta=(47, None) ===
Finishing theta (47, None) with cumC 80899400.00. Best theta so far (43, None). Best cumC 80899400.00

=== theta=(48, None) ===
Finishing theta (48, None) with cumC 80899400.00. Best theta so far (43, None). Best cumC 80899400.00

=== theta=(49, None) ===
Finishing theta (49, None) with cumC 80899400.00. Best theta so far (43, None). Best cumC 80899400.00
Finishing GridSearch in 90.31 secs
Best theta: (43, None). Best cumC: 80899400.00
CPU times: user 1min 29s, sys: 0 ns, total: 1min 29s
Wall time: 1min 30s

P.plot_cumC_chart_sample_paths( # 1 theta
    Fhat__meanI_th_I, 
    Fhat__stdvI_th_I, 
    grid_search_theta_values[1], params)

4.6.1.2 With many sample-paths

# UPDATE PARAMETERS
params.update({'Algorithm': 'GridSearch'}); pprint(f'{params=}')
params.update({'R__max': 57})
params.update({'R_0': 0})
params.update({'eta': None})

params.update({'theta_buy_min': 10})
params.update({'theta_buy_max': 50})
params.update({'theta_sell_min': None})
params.update({'theta_sell_max': None})

params.update({'L': 2_000}) #add sample-paths
params.update({'T': 100}) #optimal theta not consistent, so raise T

params

("params={'Algorithm': 'GridSearch', 'T': 100000, 'eta': None, 'R__max': 57, "
 "'R_0': 0, 'seed': 189654913, 'theta_sell_min': None, 'theta_sell_max': None, "
 "'theta_buy_min': 10, 'theta_buy_max': 50, 'theta_inc': 1, 'L': 2}")

{'Algorithm': 'GridSearch',
 'T': 100,
 'eta': None,
 'R__max': 57,
 'R_0': 0,
 'seed': 189654913,
 'theta_sell_min': None,
 'theta_sell_max': None,
 'theta_buy_min': 10,
 'theta_buy_max': 50,
 'theta_inc': 1,
 'L': 2000}

# create a model and a policy
piNames = ['X__BuyBelow']
SVarNames = ['R_t', 'D_t']
initialState = {
    'R_t': params['R_0'],
    'D_t': dem_sim.simulate(),
}
xVarNames = ['x_t_ELA']
possibleDecisions = None
p__buy = 100
p__sell = 150
exogParams = {} # we use simulation
M = InventoryStorageModel(
    SVarNames, 
    xVarNames, 
    initialState, 
    params, 
    exogParams,
    possibleDecisions,
    p__buy,
    p__sell
)
P = InventoryStoragePolicy(M, piNames)

%%time
##########################################################################
#GridSearch #. SDAM-9.4.1
if params['Algorithm'] == 'GridSearch':
    # obtain the theta values to carry out a full grid search
    grid_search_theta_values = P.grid_search_theta_values(params)
    print(grid_search_theta_values)
    
    # cumCI_theta_I = \
    #   P.perform_grid_search(params, grid_search_theta_values[0])
    # P.plot_cumC_chart( # 1 theta
    #     cumCI_theta_I, #cumC values in here
    #     grid_search_theta_values[1]) #theta values    

    Fhat__meanI_th_I, Fhat__stdvI_th_I = \
      P.perform_grid_search_sample_paths(params, grid_search_theta_values[0])
    # P.plot_cumC_chart_sample_paths( # 1 theta
    #     Fhat__meanI_th_I, 
    #     Fhat__stdvI_th_I, 
    #     grid_search_theta_values[1], params)
##################################################################################

([(10, None), (11, None), (12, None), (13, None), (14, None), (15, None), (16, None), (17, None), (18, None), (19, None), (20, None), (21, None), (22, None), (23, None), (24, None), (25, None), (26, None), (27, None), (28, None), (29, None), (30, None), (31, None), (32, None), (33, None), (34, None), (35, None), (36, None), (37, None), (38, None), (39, None), (40, None), (41, None), (42, None), (43, None), (44, None), (45, None), (46, None), (47, None), (48, None), (49, None)], array([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]))

=== theta=(10, None) ===
Finishing theta (10, None) with cumC 63038.65. Best theta so far (10, None). Best cumC 63038.65

=== theta=(11, None) ===
Finishing theta (11, None) with cumC 63736.62. Best theta so far (11, None). Best cumC 63736.62

=== theta=(12, None) ===
Finishing theta (12, None) with cumC 64589.30. Best theta so far (12, None). Best cumC 64589.30

=== theta=(13, None) ===
Finishing theta (13, None) with cumC 65279.30. Best theta so far (13, None). Best cumC 65279.30

=== theta=(14, None) ===
Finishing theta (14, None) with cumC 66158.82. Best theta so far (14, None). Best cumC 66158.82

=== theta=(15, None) ===
Finishing theta (15, None) with cumC 67371.80. Best theta so far (15, None). Best cumC 67371.80

=== theta=(16, None) ===
Finishing theta (16, None) with cumC 68193.05. Best theta so far (16, None). Best cumC 68193.05

=== theta=(17, None) ===
Finishing theta (17, None) with cumC 69334.95. Best theta so far (17, None). Best cumC 69334.95

=== theta=(18, None) ===
Finishing theta (18, None) with cumC 70223.52. Best theta so far (18, None). Best cumC 70223.52

=== theta=(19, None) ===
Finishing theta (19, None) with cumC 70975.77. Best theta so far (19, None). Best cumC 70975.77

=== theta=(20, None) ===
Finishing theta (20, None) with cumC 71782.45. Best theta so far (20, None). Best cumC 71782.45

=== theta=(21, None) ===
Finishing theta (21, None) with cumC 72417.10. Best theta so far (21, None). Best cumC 72417.10

=== theta=(22, None) ===
Finishing theta (22, None) with cumC 73019.93. Best theta so far (22, None). Best cumC 73019.93

=== theta=(23, None) ===
Finishing theta (23, None) with cumC 73773.38. Best theta so far (23, None). Best cumC 73773.38

=== theta=(24, None) ===
Finishing theta (24, None) with cumC 74644.27. Best theta so far (24, None). Best cumC 74644.27

=== theta=(25, None) ===
Finishing theta (25, None) with cumC 75581.57. Best theta so far (25, None). Best cumC 75581.57

=== theta=(26, None) ===
Finishing theta (26, None) with cumC 76790.60. Best theta so far (26, None). Best cumC 76790.60

=== theta=(27, None) ===
Finishing theta (27, None) with cumC 77707.50. Best theta so far (27, None). Best cumC 77707.50

=== theta=(28, None) ===
Finishing theta (28, None) with cumC 78298.00. Best theta so far (28, None). Best cumC 78298.00

=== theta=(29, None) ===
Finishing theta (29, None) with cumC 78471.40. Best theta so far (29, None). Best cumC 78471.40

=== theta=(30, None) ===
Finishing theta (30, None) with cumC 78692.82. Best theta so far (30, None). Best cumC 78692.82

=== theta=(31, None) ===
Finishing theta (31, None) with cumC 78692.82. Best theta so far (30, None). Best cumC 78692.82

=== theta=(32, None) ===
Finishing theta (32, None) with cumC 78692.82. Best theta so far (30, None). Best cumC 78692.82

=== theta=(33, None) ===
Finishing theta (33, None) with cumC 78692.82. Best theta so far (30, None). Best cumC 78692.82

=== theta=(34, None) ===
Finishing theta (34, None) with cumC 78692.82. Best theta so far (30, None). Best cumC 78692.82

=== theta=(35, None) ===
Finishing theta (35, None) with cumC 78692.82. Best theta so far (30, None). Best cumC 78692.82

=== theta=(36, None) ===
Finishing theta (36, None) with cumC 78925.18. Best theta so far (36, None). Best cumC 78925.18

=== theta=(37, None) ===
Finishing theta (37, None) with cumC 78934.48. Best theta so far (37, None). Best cumC 78934.48

=== theta=(38, None) ===
Finishing theta (38, None) with cumC 79020.75. Best theta so far (38, None). Best cumC 79020.75

=== theta=(39, None) ===
Finishing theta (39, None) with cumC 79063.60. Best theta so far (39, None). Best cumC 79063.60

=== theta=(40, None) ===
Finishing theta (40, None) with cumC 79063.60. Best theta so far (39, None). Best cumC 79063.60

=== theta=(41, None) ===
Finishing theta (41, None) with cumC 79063.60. Best theta so far (39, None). Best cumC 79063.60

=== theta=(42, None) ===
Finishing theta (42, None) with cumC 79063.60. Best theta so far (39, None). Best cumC 79063.60

=== theta=(43, None) ===
Finishing theta (43, None) with cumC 79063.60. Best theta so far (39, None). Best cumC 79063.60

=== theta=(44, None) ===
Finishing theta (44, None) with cumC 79063.60. Best theta so far (39, None). Best cumC 79063.60

=== theta=(45, None) ===
Finishing theta (45, None) with cumC 79063.60. Best theta so far (39, None). Best cumC 79063.60

=== theta=(46, None) ===
Finishing theta (46, None) with cumC 79063.60. Best theta so far (39, None). Best cumC 79063.60

=== theta=(47, None) ===
Finishing theta (47, None) with cumC 79063.60. Best theta so far (39, None). Best cumC 79063.60

=== theta=(48, None) ===
Finishing theta (48, None) with cumC 79063.60. Best theta so far (39, None). Best cumC 79063.60

=== theta=(49, None) ===
Finishing theta (49, None) with cumC 79063.60. Best theta so far (39, None). Best cumC 79063.60
Finishing GridSearch in 95.60 secs
Best theta: (39, None). Best cumC: 79063.60
CPU times: user 1min 31s, sys: 218 ms, total: 1min 31s
Wall time: 1min 35s

P.plot_cumC_chart_sample_paths( # 1 theta
    Fhat__meanI_th_I, 
    Fhat__stdvI_th_I, 
    grid_search_theta_values[1], params)

4.6.2 Evaluation

M_evalu = InventoryStorageModel(
    SVarNames, 
    xVarNames, 
    initialState, 
    params, 
    exogParams,
    possibleDecisions,
    p__buy,
    p__sell
)
P_evalu = InventoryStoragePolicy(M_evalu, piNames)
params

{'Algorithm': 'GridSearch',
 'T': 100,
 'eta': None,
 'R__max': 57,
 'R_0': 0,
 'seed': 189654913,
 'theta_sell_min': None,
 'theta_sell_max': None,
 'theta_buy_min': 10,
 'theta_buy_max': 50,
 'theta_inc': 1,
 'L': 2000}

def run_policy_evalu(piInfo_evalu, piName_evalu, stop_time_evalu):
    model_copy = copy(M_evalu)
    record = []
    for t in range(stop_time_evalu):
        x_t = getattr(P_evalu, piName_evalu)(t, model_copy.S_t, piInfo_evalu, stop_time_evalu)        
        res = model_copy.step(t, x_t) # step the model forward one iteration
        record.append([res[0].R_t, res[0].D_t, res[1], res[2].x_t_ELA])
    cumC = model_copy.cumC    
    return cumC, record

# EVALUATION
piName_evalu = 'X__BuyBelow'
stop_time_evalu = 100 #180

4.6.2.1 Non-optimal policy

theta_evalu_non=(10, None)
cumC, record = run_policy_evalu(theta_evalu_non, piName_evalu, stop_time_evalu)
labels = ["R_t", "D_t", "cumC", 'x_t_ELA']
print(f'{theta_evalu_non=}')
df_non = pd.DataFrame.from_records(data=record, columns=labels); df[:10]

theta_evalu_non=(10, None)

	demand
0	6
1	6
2	9
3	10
4	10
5	8
6	14
7	12
8	12
9	13

4.6.2.2 Optimal policy

from certifi.core import where
def plot_output(df1, df2):
  legendlabels = [r'$\mathrm{opt}$', r'$\mathrm{non}$']
  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(10)
  fig.set_figwidth(13); fig.set_figheight(9)
  fig.suptitle('PERFORMANCE OF OPTIMIZED Buy-Below POLICY\nOptimal (magenta), Non-optimal (cyan)', fontsize=20)

  i = 0 #x_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['x_t_ELA'], 'm')
  axs[i].step(df2['x_t_ELA'], 'c')
  axs[i].axhline(y=0, color='k', linestyle=':')
  axs[i].set_ylabel('$x_{t,ELA}$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);

  i = 1 #D_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['D_t'], 'k')
  axs[i].set_ylabel('$D_t$'+'\n'+'$\mathrm{[units]}$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);

  i = 2 #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'], '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].step(df2['R_t'], 'c')
  axs[i].axhline(y=theta_evalu_non[0], color='c', linestyle=':')
  axs[i].text(0, theta_evalu_non[0], r'$\theta^{buy}$', size=16)

  axs[i].set_ylabel('$R_t$'+'\n'+'$\mathrm{[units]}$', 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'], 'm')
  axs[i].step(df2['cumC'], 'c')
  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{[order\ windows]}$', rotation=0, ha='right', va='center', fontweight='bold', size=ylabelsize);

  fig.legend(labels=legendlabels, loc='center', fontsize=18)

4.6.2.2.1 Optimal \(θ\) from training with a few sample-paths

# theta_evalu = (46, None)
# theta_evalu = (45, None)
theta_evalu = (43, None)
cumC, record = run_policy_evalu(theta_evalu, piName_evalu, stop_time_evalu)
labels = ["R_t", "D_t", "cumC", 'x_t_ELA']
print(f'{theta_evalu=}')
df = pd.DataFrame.from_records(data=record, columns=labels); df[:10]

theta_evalu=(43, None)

	R_t	D_t	cumC	x_t_ELA
0	57	16	-3,900.0000	57
1	41	15	-1,500.0000	0
2	42	20	-850.0000	16
3	37	15	-400.0000	15
4	42	14	-150.0000	20
5	43	11	300.0000	15
6	46	13	1,600.0000	14
7	33	15	3,700.0000	0
8	42	15	3,850.0000	24
9	42	17	4,600.0000	15

plot_output(df, df_non)

4.6.2.2.1 Optimal \(θ\) from training with many sample-paths

# theta_evalu = (37, None)
# theta_evalu = (38, None)
theta_evalu = (39, None)
cumC, record = run_policy_evalu(theta_evalu, piName_evalu, stop_time_evalu)
labels = ["R_t", "D_t", "cumC", 'x_t_ELA']
print(f'{theta_evalu=}')
df = pd.DataFrame.from_records(data=record, columns=labels); df[:10]

theta_evalu=(39, None)

	R_t	D_t	cumC	x_t_ELA
0	57	19	-3,000.0000	57
1	38	20	300.0000	0
2	37	19	1,250.0000	19
3	38	21	2,100.0000	20
4	36	17	2,900.0000	19
5	40	18	3,050.0000	21
6	22	17	5,000.0000	0
7	40	17	3,300.0000	35
8	23	13	5,400.0000	0
9	44	12	3,800.0000	34

plot_output(df, df_non)