# Model-based controlled learning of MDP policies
with an application to lost-sales inventory control

Willem van Jaarsveld
Eindhoven University of Technology

###### Abstract

Recent literature established that neural networks can represent good MDP
policies across a range of stochastic dynamic models in supply chain and
logistics. To overcome limitations of the model-free algorithms typically
employed to learn/find such neural network policies, a model-based algorithm
is proposed that incorporates variance reduction techniques. For the classical
lost sales inventory model, the algorithm learns neural network policies that
are superior to those learned using model-free algorithms, while also
outperforming heuristic benchmarks. The algorithm may be an interesting
candidate to apply to other stochastic dynamic problems in supply chain and
logistics.
Keywords: Markov Decision Processes, Neural Network Policies, Variance
Control, Inventory management

##  1 Introduction

Markov decision processes (MDPs) are a natural model for countless stochastic
dynamic decision problems in supply chain, transportation, manufacturing,
maintenance and service logistics, etc. However, optimal MDP policies can only
be computed in special cases where the state space is small. For decades,
researchers have developed approximations to find MDP policies with
satisfactory performance. Breakthroughs in Deep Reinforcement Learning (DRL)
for Atari games [[, e.g. A3C:]]mnih2016asynchronous and go, chess and shogi
[[]]silver2018general have demonstrated the capabilities of modern hardware
while contributing algorithmic ideas that have intensified the quest for such
approximations.

For operations research (OR), the most exciting aspect of these breakthroughs
has become clear through the work of [14]: The A3C algorithm is of direct
value for difficult stochastic dynamic problems in supply chain and logistics.
In particular, [14] consider a range of difficult problems in inventory
control (lost sales, dual sourcing, and one-warehouse multi-retailer) and show
that A3C yields policies that are competitive with hand-crafted heuristics
that took decades to develop. This remarkable result demonstrates that neural
networks can represent good policies for some of the most challenging
stochastic decision problems that arose in supply chain and logistics over the
past decades.

Algorithms like A3C were developed with certain benchmark applications in
mind: e.g. Atari games, action control in physics simulators, etc. Like
stochastic decision problems arising in logistics and supply chain, these
benchmarks can be formulated as an MDP and thus fall within the DRL framework,
but there is a key difference: The environments underlying these benchmarks
are nearly deterministic, while MDPs that arises in supply chain and logistics
feature transitions that are highly stochastic (e.g. because of random demand
arrivals). In the former case, series of good actions typically yield
trajectories with favorable objectives. In the latter case, this statement
only holds on average, and individual trajectories may result in high costs
even while optimal actions were taken for every state. This is relevant
because DRL algorithms like A3C rely solely on simulating state trajectories
to gather statistical evidence on the quality of actions in states. Based on
high-variance statistics, distinguishing between good and near-optimal actions
is cumbersome, and if we do not _control_ this variance, we may need to settle
for good actions. In this light, it is worth noting that the optimality gaps
reported for A3C and related algorithms for difficult problems in the
logistics domain are typically in the order of 1-3% [14, 41], and that
training is long and may require extensive tuning.

To make a further step towards the adoption of neural networks as general-
purpose policies for MDPs arising in e.g. logistics and supply chain, this
paper proposes and demonstrates an algorithm for learning good MDP policies
that aims to address some of the shortcomings of _model-free_ DRL algorithms
like A3C when dealing with highly stochastic MDPs. Like in A3C, policies are
represented as neural networks and the aim is finding neural network
parameters such that the resulting policy is close-to-optimal. In other
aspects, the proposed algorithm is rather different from model-free DRL
algorithms like A3C: It is an iterative procedure based on labelled samples,
see §2 for details and related work. To deal with the high variance
trajectories that arise in decision models for logistics and supply chain, a
generic approach for controlling variance is proposed that vastly improves the
ability to obtain high-quality samples. The approach leverages structural
information on the impact of uncertain environment factors on model dynamics
that may not be available for all MDPs. However, such information _is_
available for many problems in logistics and supply chain, which renders the
algorithm particularly attractive for such problems.

It is demonstrated that the approach, termed _model-based controlled learning_
(MCL), finds near-optimal neural network policies relatively quickly. In
particular, for the classical lost sales inventory problem [44], the algorithm
is shown to vastly outperform both the best neural-network based benchmark
(A3C) as well as a range of heuristic benchmarks that have been developed over
the past decades. Additionally, the tuning of parameters appears to be a
relatively simple task for the proposed algorithm.

The exposition in this paper is tailored to the operations research/operations
management community: apart from some familiarity with Markov Decision
Processes, no pre-knowledge (e.g. on neural networks) is assumed. The
remainder of this paper is organized as follows: §2 reviews literature, §3
contains preliminaries, and §4 discusses the MCL algorithm. Experiments for
the classical lost sales inventory problem are presented in §5, while §6
discusses limitations and future research.

##  2 Literature

Algorithms that _approximate_ the optimal MDP have been studied for decades
under several essentially equivalent terminologies: approximate dynamic
programming, neuro-dynamic programming, and reinforcement learning [[,
cf.]]bertsekas2019reinforcement. These fields are vast: Our selective review
focuses on applications in the operations research/operations management
(OR/OM) domain, while placing the proposed methodology in the perspective of
extant work.

Two main approaches in approximate dynamic programming are value function
approximation (VFA) and state aggregation [[, see
e.g.]]powell2007approximate,bertsekas2011approximate,mes2017approximate. These
approaches have been successfully applied to a wide range of problem domains.
VFA typically assumes a particular form of the value function (e.g. linear in
the features). Some recent applications include grid-level energy storage
[33], lateral transhipments [24], response-adaptive clinical trial design [1],
rehabilitation treatment planning [10] and medical evacuation dispatching
[18]. In special cases, the value function may have a known structural form
which can be exploited. E.g. for certain inventory applications, the value
function is $L\natural$ convex, leading to highly specialized ADP applications
for lost sales [37] and dual sourcing [11]. For state aggregration, likewise a
range of approaches is available, and tailoring the approach to the
application may give substantial benefits; examples include [13] and [42]. In
this stream of literature, neural networks have been used as value function
approximators [[, e.g.]]van1997neuro,jenkins2020approximate and for dynamic
state aggregation [[, e.g.]]bertsekas2018feature.

Another approach to ADP/RL are parametric policy function approximations.
These have a long history [[, for context see e.g.]]powell2019unified, and the
adoption of neural networks for this purpose was popularized through recent
advances under the banner of Deep Reinforcement Learning (DRL). [28] and [14]
were amongst the first to adopt DRL algorithms in OR/OM. [28] use DQN [27] for
the beer game, and their unique multi-agent perspective enables a
demonstration of e.g. the ability of DQN algorithms to adapt to suboptimal
behavior of other agents in the supply chain. [14] explore the feasibility of
A3C as a _general-purpose_ solution methodology for operations management
problems. They show that A3C learns policies for three difficult OM problems
in inventory control (lost sales, dual sourcing, one warehouse multi retailer)
with optimality gaps ranging from $1-3\%$: competitive with well-known
heuristics. [41] test proximal-policy optimization (PPO, [34]) for the joint
replenishment problem and show that neural network policies outperform
heuristics, and that the structure of learned neural networks policies
resembles the structure of the optimal policy. [29] adopt DQN to study
simultaneous ordering in general multi-echelon supply chains.

A3C and PPO train neural networks directly based on the cost/reward outcomes
of trajectories. There are also algorithms that train neural networks from
labelled samples, i.e. from combinations of states and actions belonging to
that state. Such algorithms are typically _model-based_ , in the sense that
they need access to a simulator of the environment for training. An example of
such an approach is _AlphaZero_ [[]]silver2018general, which focuses on two-
player deterministic games and obtains labelled samples using Monte-Carlo tree
search. Another approach to obtain labelled samples is from _roll-out_
policies [23, 12]. In either case, the algorithm may alternate between sample
collection and supervised training of a neural network. Roll-out policies have
been used in ADP/RL literature for decades [[, an early example in OR/OM
is]]secomandi2001rollout. However, it appears that _iterative_ improvement of
neural network policies from labelled samples [[,
cf.]]lagoudakis2003reinforcement has not had a very broad take-up outside the
AI community. In fact there appear to be no prior applications in inventory
control, logistics, or to other problems in the OR/OM domain. In robotics [7]
independently proposed a model-based algorithm for learning policies for
autonomous sequential repair.

The algorithm proposed in this paper employs roll-outs for collecting labelled
samples, and trains several generations of neural networks in an iterative
fashion. However, in highly stochastic settings, standard roll-outs as
proposed by [23] are ineffective. To address this, the algorithm incorporates
a variance reduction scheme that greatly improves the computational
tractability. The variance reduction that is employed is _model-based_ : It
leverages specific information about the impact of uncertainty on model
trajectories that is typically available for models arising from decision
problems in logistics, supply chains and various other domains. The author is
not aware of prior implementations of such model-based variance reduction
techniques for the iterative training of neural networks. Note that the idea
is perhaps rather obvious once high variance is identified as the bottleneck
for DRL performance: variance reduction in the context of 1-step policy
improvement is not new.

The main contribution of this paper is proposing a concrete algorithm that
combines the above ideas, and demonstrating that this model-based controlled
learning (MCL) algorithm vastly outperforms model-free algorithms like A3C for
a typical stochastic decision problem in the logistics and supply chain
domain. To this end, the classical _lost sales inventory systems_ is adopted
as test environment. This system is amongst the most well-studied _difficult_
stochastic problems in supply chain and logistics, and it is interesting both
from a practical and a theoretical perspective [9]. Breakthrough analytical
results [17, 15, 43] have demonstrated _asymptotic_ optimality for certain
simple policies. However, outside asymptotic regimes the only known approach
for finding optimal policies relies on explicitly solving the Bellman
optimality equations, which is intractable even for moderately long leadtimes.
Therefore, companies traditionally rely on heuristics, and an excellent set of
benchmark problems to test such heuristics is available from [44]; an
additional contribution of this paper is that the MCL algorithm yields very
competitive policies for this benchmark.

##  3 Preliminaries

MDP theory is reviewed in §3.1, while §3.2 discusses neural networks.

###  3.1 Markov decision processes

Consider a Markov Decision Process (MDP) [[, cf.]]puterman2014markov with a
finite action space $\mathcal{A}=\\{1,\ldots,m\\}$ and a discrete state space
$\mathcal{S}$. Let $\mathcal{A}_{\mathbf{s}}\subseteq\mathcal{A}$ denote the
actions that are applicable for state $\mathbf{s}\in\mathcal{S}$. (Often,
$\mathcal{A}_{\mathbf{s}}=\mathcal{A}$ for all $\mathbf{s}$, but in some cases
not all actions apply for all states.) Costs are accumulated over an infinite
horizon and are discounted by $\alpha\in(0,1)$. We discuss these, and other
assumptions in §6. Denote by $\mathbf{P}(\mathbf{s}^{\prime}|\mathbf{s},a)$
the probability to transition from state $\mathbf{s}$ to state
$\mathbf{s}^{\prime}$ when taking action $a\in\mathcal{A}_{\mathbf{s}}$.
Taking action $a$ in state $\mathbf{s}$ yields costs $C(\mathbf{s},a)$. Let
$\bar{c}(\mathbf{s},a)=\mathbf{E}C(\mathbf{s},a)$. A policy $\pi$ is a
function $\pi:\mathcal{S}\to\mathcal{A}$ that maps each state $\mathbf{s}$ to
a valid action $\pi(\mathbf{s})\in\mathcal{A}_{s}$. We seek policies that
minimize the total discounted costs.

The well-known _value function_ $v_{\pi}(\cdot)$ represents the expected total
discounted costs over an infinite horizon, when starting from state
$\mathbf{s}$ and following $\pi$; it satisfies

$\displaystyle\forall s\in
S:v_{\pi}(s)=\bar{c}(s,\pi(s))+\alpha\sum_{\mathbf{s}^{\prime}\in\mathcal{S}}\mathbf{P}(\mathbf{s}^{\prime}|\mathbf{s},\pi(\mathbf{s}))v_{\pi}(\mathbf{s}^{\prime})$
(1)

Similarly, denote by $q_{\pi}(\mathbf{s},a)$ the expected total discounted
costs when starting from state $\mathbf{s}$, taking action $a$, and following
$\pi$ afterwards. For $\mathbf{s}\in\mathcal{S}$ and $a\in\mathcal{A}_{s}$, it
holds that
$q_{\pi}(\mathbf{s},a)=\bar{c}(s,a)+\alpha\sum_{\mathbf{s}^{\prime}\in\mathcal{S}}\mathbf{P}(\mathbf{s}^{\prime}|\mathbf{s},a)v_{\pi}(\mathbf{s}^{\prime})$.

Given any policy $\pi(\cdot)$, define the _improved_ policy $(\pi)^{+}(\cdot)$
by
$(\pi)^{+}(\mathbf{s}):=\operatorname*{arg\,min}_{a\in\mathcal{A}_{\mathbf{s}}}q_{\pi}(\mathbf{s},a)$.
Repeatedly improving a policy leads to the optimal policy under a wide range
of assumptions [32], but policy improvement is computationally intractable for
large state spaces.

###  3.2 Neural networks

For the purposes of this paper, it is sufficient to think of neural networks
as parameterized functions from $\mathbb{R}^{\nu}$ to $\mathbb{R}^{\mu}$, for
some $\nu,\mu\in\mathbb{N}$. Denote such a (generic) function by
$N_{\theta}(\cdot)$, where $\theta$ denotes the function parameters. That is,
for each set of parameters $\theta$ and each $\mathbf{x}\in\mathbb{R}^{\nu}$,
we have $\mathbf{y}=N_{\theta}(\mathbf{x})\in\mathbb{R}^{\mu}$.

While algorithms in this paper are not restricted to any specific neural
network, experiments in this paper utilize the “vanilla” neural network: The
multi-layer perceptron (MLP). For completeness and concreteness, we next give
$N_{\theta}(\cdot)$ for the MLP [[, for details, see e.g.]]kroese2019DSML.
Like most neural networks, the MLP is composed of multiple _layers_ ; each
layer corresponds to a mathematical operation that operates on the output of
the previous layer. In the MLP, each layer $l\in\\{1,\ldots,L\\}$ combines an
affine transformation of the inputs with a non-linear _activation function_.
Accordingly, layer $l$ is parameterized by a matrix $A^{l}$ and a vector
$\mathbf{b}^{l}$, thus the parameters of the entire network are
$\theta=\\{A^{1},\mathbf{b}^{1},\ldots,A^{L},\mathbf{b}^{L}\\}$. The
activation function $f^{l}(\cdot):\mathbb{R}^{k}\to\mathbb{R}^{k}$ maps
vectors to vectors. In this paper we adopt
$f^{l}(\mathbf{x})=\max(\mathbf{x},0)$ for $l<L$, where the $\max$ is taken
component-wise, while $f^{L}(\mathbf{x})=\mathbf{x}$. More precisely, layer
operations in the neural networks are governed by:

$\displaystyle\mathbf{a}^{l}$
$\displaystyle=f^{l}(A^{l}\mathbf{a}^{l-1}+\mathbf{b}^{l})$ $\displaystyle
l\in\\{1,\ldots,L\\}$ (2)

The MLP sets $N_{\theta}(\mathbf{x}):=\mathbf{a}^{L}$, while interpreting
$\mathbf{a}^{0}:=\mathbf{x}$, i.e.
$\mathbf{a}^{1}=f^{1}(A^{1}\mathbf{x}+\mathbf{b}^{1})$.

Validity of (2) is ensured by requiring that for $l\in\\{2,\ldots,L\\}$ the
number of columns in $A^{l}$ equals the number of rows in $A^{l-1}$, etc.
Also, the choice of dimensions of $A^{l}$ and $\mathbf{b}^{l}$ governs the
dimensions of input layer $\mathbf{a}^{0}$, output layer $\mathbf{a}^{L}$, and
hidden layers $\mathbf{a}^{l}$ ($l\in\\{1,\ldots,L-1\\})$. These dimensions,
together with $L$, specify the _structure_ of the neural network.

###  3.3 Neural network policies

Neural networks can represent MDP policies. To make this explicit, note that
MDP states $\mathbf{s}$ can typically be represented as vectors in
$\mathbb{R}^{n}$ for some $n\in\mathbb{N}$. Accordingly, and for ease of
presentation, suppose that $\mathcal{S}\subset\mathbb{R}^{n}$. Also, remember
that $\mathcal{A}=\\{1,\ldots,m\\}$.

Now consider any neural network
$N_{\theta}(\cdot):\mathbb{R}^{\nu}\to\mathbb{R}^{\mu}$ such that $\nu=n$ and
$\mu=m$. For any $\theta$, such network corresponds to an MDP policy that we
shall denote by $\pi_{\theta}(\cdot)$, and that is defined as follows:

$\displaystyle\pi_{\theta}(\mathbf{s})=\operatorname*{arg\,max}_{a\in\mathcal{A}_{\mathbf{s}}}[(N_{\theta}(\mathbf{s}))_{a}]$
(3)

Here, $(\mathbf{y})_{a}$ denotes the $a$th element of the vector $\mathbf{y}$;
ties in the $\operatorname*{arg\,max}$ are broken by selecting the minimum $a$
that attains the maximum. This makes (3) well-defined since
$N_{\theta}(\mathbf{s})\in\mathbb{R}^{m}$ and
$\mathcal{A}_{\mathbf{s}}\subseteq\\{1,\ldots,m\\}$.

While defining a policy based on a neural network may seem arbitrary, it has
been demonstrated in prior research that neural networks have the flexibility
to represent good policies across domains; for a notable example in
logistics/supply chain see [14]. (Model-free algorithm such as A3C use a
slightly different approach for interpreting a neural network as a policy: The
output of the neural network is transformed into a probability distribution
over the action space, by taking the so-called soft-max of the outputs. Thus
A3C optimizes stochastic policies. )

In line with DRL literature, the neural network _structure_ is exogenous to
the model-based controlled learning (MCL) algorithm: For _given_ network
structure, MCL aims to find parameters $\theta$ _such_ that
$\pi_{\theta}(\cdot)$ is a well-performing policy.

###  3.4 Supervised learning

The MCL algorithm proposed in the next section relies on _supervised_ learning
to train neural networks. The following brief introduction to supervised
learning can be safely skipped by readers familiar with the subject.

In our context, supervised learning finds a _policy_ , i.e. a relation between
the actions $\pi(\mathbf{s})\in\mathcal{A}_{\mathbf{s}}$ taken by the policy
and the state variables $\mathbf{s}$. Thus the components of the state
$\mathbf{s}\in\mathbb{R}^{n}$ act as features. In fact, since our action space
is discrete, we are learning a classifier [23]. Suppose we obtain a set of
_data points_
$\mathcal{K}=\\{(\mathbf{s}_{k},\pi(\mathbf{s}_{k}))|k\in\\{1,\ldots,K\\}\\}$
that relate states to the action taken in those states by some policy
$\pi(\cdot)$. A vast range of algorithms have been proposed [[, see
e.g.]]goodfellow2016deep that aim to fit (or learn) neural network parameters
$\theta$ such that $\pi_{\theta}$ is _similar_ to $\pi$, e.g. $\pi$ and
$\pi_{\theta}$ tend to take the same actions in the same states, _also_ for
states that are not represented in the data points $\mathcal{K}$. (For
_general_ policies $\pi:\mathcal{S}\rightarrow\mathcal{A}$ one obviously
cannot expect $\pi$ and $\pi_{\theta}$ to be the same for states not in
$\mathcal{K}$. However, _good_ policies $\pi$ may exhibit structure that can
be captured by the neural network.)

Learning $\theta$ may involve splitting $\mathcal{K}$ in a training set and a
test set. The _training loss_ $\mathcal{L}(\theta)$ measures the “distance”
between $\pi$ and $\pi_{\theta}$ for the states in the training set. Fitting
$\theta$ is an iterative, gradient-based process: In each step, the gradient
of $\mathcal{L}(\theta)$ with respect to $\theta$ is estimated, and
subsequently $\theta$ is updated by taking a step in a direction opposite to
the gradient. An approach to avoid overfitting on the training data is to
terminate the iterative procedure when the loss on the test set, defined
analogously to the loss for the training set, no longer decreases. For
completeness, Appendix A describes the training approach adopted for the
experiments presented in §5. For later use, we will denote neural network
parameters $\theta$ trained on data set $\mathcal{K}$ by
$\theta(\mathcal{K})$.

##  4 Algorithm

###  4.1 Overview

Consider an MDP as introduced in §3.1, and some initial policy $\pi_{0}$. A
few _policy improvement_ steps (cf. §3.1) on $\pi_{0}$ typically yield a
(very) good policy, but exact policy improvement is intractable when
$|\mathcal{S}|$ is large. We shall instead make a few _approximate_ policy
improvement steps $i=1,2,\ldots$, each involving both Monte Carlo simulation
and the training of a neural network. Denote the policy resulting from step
$i$ by $\pi_{i}$. Next follows an overview of how policy $\pi_{i+1}$ is
obtained from policy $\pi_{i}$; details follow in §4.2 and 4.3.

Remember that an exact policy improvement step on $\pi_{i}$ yields policy
$(\pi_{i})^{+}$, where
$(\pi_{i})^{+}(\mathbf{s}):=\operatorname*{arg\,min}_{a\in\mathcal{A}_{\mathbf{s}}}q_{\pi_{i}}(\mathbf{s},a)$.
However, computing $q_{\pi_{i}}(\mathbf{s},a)$ is typically intractable for
$|\mathcal{S}|$ large. Note, however, that $q_{\pi}(\mathbf{s},a)$ is the
expected total discounted reward when taking action $a$ in state $\mathbf{s}$,
and subsequently following policy $\pi$. In particular, let
$Q_{\pi}(\mathbf{s},a)$ denote the discounted costs incurred in a sample path
that starts with action $a$ from state $\mathbf{s}$ and follows $\pi(\cdot)$
onwards:

$\displaystyle Q_{\pi}(\mathbf{s},a)$
$\displaystyle:=\sum_{t=1}^{\infty}\alpha^{t-1}c_{t}$ $\displaystyle\text{
with }\mathbf{s}_{0}=\mathbf{s},a_{0}=a,\forall t\geq 0:c_{t}\sim
C(\mathbf{s}_{t},a_{t}),s_{t+1}\sim\mathbf{P}(\cdot|\mathbf{s}_{t},\pi(\mathbf{s}_{t}))$
(4)

Then $\mathbf{E}Q_{\pi}(\mathbf{s},a)=q_{\pi}(\mathbf{s},a)$. Accordingly, as
detailed in §4.2, Monte-Carlo techniques can be used to obtain unbiased
estimators of the costs $q_{\pi}(\mathbf{s},a)$. Such estimates underlie a
simulation-based approximation of

$(\pi_{i})^{+}(\mathbf{s}):=\operatorname*{arg\,min}_{a\in\mathcal{A}_{\mathbf{s}}}\mathbf{E}Q_{\pi_{i}}(\mathbf{s},a)$

for any state $\mathbf{s}$, as discussed in detail in §4.2. Abusing notation,
denote the simulation-based action retrieved for state $\mathbf{s}$ in this
fashion by $\hat{\pi}^{+}_{i}(\mathbf{s})$.

Finding $\hat{\pi}_{i}^{+}(\mathbf{s})$ for a state $\mathbf{s}$ involves
extensive simulations and is computationally expensive. Therefore, setting
$\pi_{i+1}$ equal to $\hat{\pi}_{i}^{+}$ would render further approximate
policy improvement steps intractable. Instead, arriving at policy $\pi_{i+1}$
involves collecting data points
$\mathcal{K}_{i}=\\{(\mathbf{s}_{k},\hat{\pi}_{i}^{+}(\mathbf{s}_{k}))|k\in\\{1,\ldots,K\\}\\}$
and subsequently learning neural network parameters $\theta(\mathcal{K}_{i})$
from this data (cf. §3.4). The approximate policy improvement step is then
obtained as $\pi_{i+1}=\pi_{\theta(\mathcal{K}_{i})}$. In summary, the
approximate policy improvement on policy $\pi_{i}$ involves collecting sample
actions $\mathcal{K}_{i}$ from an approximate, simulation-based
_approximation_ of $(\pi_{i})^{+}$, and subsequently using that data to train
a neural network to obtain parameters $\theta(\mathcal{K}_{i})$, resulting in
policy $\pi_{i+1}=\pi_{\theta(\mathcal{K}_{i})}$. Since exact policy iteration
typically converges to close-to-optimal solutions in 3-6 steps for stochastic
models, a few approximate policy iterations steps may yield a policy with good
performance.

The remainder of this section is structured as follows. §4.2 discusses the
simulation-based approach that retrieves $\hat{\pi}_{i}^{+}(\mathbf{s})$ for
any given state $\mathbf{s}$. §4.3 discusses how to select the states that are
to be included in $\mathcal{K}_{i}$.

###  4.2 Simulation-based policy improvement

The generic approach for controlling variance while estimating
$(\pi_{i})^{+}(\mathbf{s})$ is discussed in §4.2.1, while the adopted
algorithm is discussed in §4.2.2.

####  4.2.1 Simulation estimators

Suppose we fix some $\mathbf{s}\in\mathcal{S}$ and some policy $\pi$. Consider
the following approach: 1) Obtain several independent replications of
$Q_{\pi}(\mathbf{s},a)$ for each $a\in\mathcal{A}_{\mathbf{s}}$ 2) Let
$\hat{Q}_{\pi}(\mathbf{s},a)$ be the respective average of these replications,
and 3) then $(\pi_{i})^{+}(\mathbf{s})$ can be approximated as
$\operatorname*{arg\,min}_{a\in\mathcal{A}_{s}}\hat{Q}_{\pi}(\mathbf{s},a)$.

This approximation is asymptotically exact when sufficiently many replications
of $Q_{\pi}(\mathbf{s},a)$ are collected in Step 1, but in practice the
approach does not work well when directly applied for highly stochastic MDPs
arising in operations management. Indeed, let
$a,a^{\prime}\in\mathcal{A}(\mathbf{s})$, and suppose
$(\pi_{i})^{+}(\mathbf{s})=a$ and $q_{\pi}(\mathbf{s},a)\approx
q_{\pi}(\mathbf{s},a^{\prime})$, i.e. actions $a$ and $a^{\prime}$ are both
reasonably good. Since sample path costs $Q_{\pi}(\mathbf{s},a)$ and
$Q_{\pi}(\mathbf{s},a)$ tend to have a high variance for problems in
stochastic OM, determining with reasonable accuracy whether $a$ or
$a^{\prime}$ is preferred may require _too many_ replications. The remainder
of this section discusses an approach to overcome this.

The approach simultaneously addresses another (unrelated) difficulty
pertaining to estimation of $q(\mathbf{s},a)$ via simulation:
$Q_{\pi}(\mathbf{s},a)$ corresponds to the discounted costs accumulated over
an _infinite_ horizon. However, note that

$\displaystyle\mathbf{E}\sum_{t=0}^{\infty}\alpha^{t}c_{t}=\mathbf{E}\sum_{t=0}^{T}c_{t},$
(5)

where $T$ is geometrically distributed with parameter $\alpha$
($T\sim\operatorname*{Geo}(\alpha)$). Hence, unbiased estimators of
$q_{\pi}(\mathbf{s},a)$ may be obtained by collecting the _total_ reward over
a trajectory of $T\sim\operatorname*{Geo}(\alpha)$ periods, where the initial
state and action are $\mathbf{s}$ and $a$ and where later actions are governed
by $\pi$.

#### Conceptual approach

Model-based _controlled_ learning departs from the assumption that the
randomness in transitions and costs are caused by uncertain factors or events
that influence the trajectory, and that exist _independently_ of the
trajectory, at least in the MDP model. There are many examples of models where
uncertain factors exist independently: When controlling inventory in classical
or data-driven models, demand is exogenous [44], when designing clinical
trials, patient outcomes are exogenous [1], when dispatching military medical
evacuations assets, evacuations requests are exogenous [18].

Let $\xi$ denote a (composite) random variable that comprises the uncertain
factors that may influence the trajectory (including trajectory length $T$). A
concrete example implementation of $\xi$ is presented below. The total costs
of a trajectory that encounters the uncertain factors $\xi$ will be denoted by
$Q_{\pi}(\mathbf{s},a|\xi)$, where $\mathbf{s}$ and $a$ denote the first state
and action, and $\pi$ governs subsequent actions. Trajectory randomness is
assumed to be comprised in $\xi$, thus
$\mathbf{E}_{\xi}(Q_{\pi}(\mathbf{s},a|\xi))=q_{\pi}(\mathbf{s},a)$. To reduce
the number of replications required, we fix $\xi$ while comparing different
actions $a$ for the same start state $\mathbf{s}$. (Concretely, in terms of
the inventory application discussed below, this implies comparing trajectories
for different actions using the _same_ demand sequence, instead of generating
trajectories independently.)

To see how this would help, note that estimating
$\operatorname*{arg\,min}_{a\in\mathcal{A}_{\mathbf{s}}}q_{\pi}(\mathbf{s},a)$
involves assessing whether
$q_{\pi}(\mathbf{s},a)-q_{\pi}(\mathbf{s},a^{\prime})$ is bigger or smaller
than zero, for all pairs of actions $a,a^{\prime}\in\mathcal{A}_{\mathbf{s}}$.
As an illustration, consider two actions
$a,a^{\prime}\in\mathcal{A}_{\mathbf{s}}$. Let $\xi_{i}$ for $i=1,\ldots,r$
denote independent replications of $\xi$, and let $X$ denote the resulting
estimator of $q_{\pi}(\mathbf{s},a)-q_{\pi}(\mathbf{s},a^{\prime})$. I.e.
$X=\frac{1}{r}\sum_{i=1}^{r}Q_{\pi}(\mathbf{s},a|\xi_{i})-Q_{\pi}(\mathbf{s},a^{\prime}|\xi_{i})$.
We have $\mathbf{E}[X]=q_{\pi}(\mathbf{s},a)-q_{\pi}(\mathbf{s},a^{\prime})$
and

$\operatorname{var}[X]=\frac{1}{r}\big{[}\operatorname{var}[Q_{\pi}(\mathbf{s},a|\xi)]+\operatorname{var}[Q_{\pi}(\mathbf{s},a^{\prime}|\xi)]-2\operatorname{cov}[Q_{\pi}(\mathbf{s},a|\xi),Q_{\pi}(\mathbf{s},a^{\prime}|\xi)]\big{]}.$

The covariance term, which appears only because of the pairwise comparison of
trajectory costs while fixing $\xi$, may be substantial compared to the
variance terms, especially for highly stochastic MDPs. As a consequence,
variance control may substantially reduce the number of replications $r$
required to get the variance of estimator $X$ to an appropriate level.

#### A concrete approach

For a concrete approach to obtain $\xi$ and $Q(\mathbf{s},a|\xi)$, suppose
that costs and transitions from state $\mathbf{s}_{t}$ in period $t$ are
completely determined by the action $a_{t}$ and by some _exogenous_ random
variable $W_{t}$, where $W_{t}$ are i.i.d. random variables with common
distribution $W$ [[, cf.]]bertsekas1997differential. (E.g. $W_{t}$ could be
the per-period random demand.) More precisely, suppose that $W$ takes on
values in some set $\mathcal{W}$, and that there exist functions
$f:\mathcal{S}\times\mathcal{A}\times\mathcal{W}\to\mathcal{S}$ and
$g:\mathcal{S}\times\mathcal{A}\times\mathcal{W}\to\mathbb{R}$ such that
taking action $a$ in state $\mathbf{s}$ results in a transition to random
state $f(\mathbf{s},a,W)$ and random costs $g(\mathbf{s},a,W)$. We may relate
$f$ and $g$ to the standard notation of §3.1:

$\displaystyle\mathbf{E}_{W}[g(\mathbf{s},a,W)]=\bar{c}(\mathbf{s},a)$
$\displaystyle\mathbf{s}\in\mathcal{S},a\in\mathcal{A}_{\mathbf{s}}$ (6)
$\displaystyle\mathbf{P}_{W}[f(\mathbf{s},a,W)=\mathbf{s}^{\prime}]=\mathbf{P}(\mathbf{s}^{\prime}|\mathbf{s},a)$
$\displaystyle\mathbf{s},\mathbf{s}^{\prime}\in\mathcal{S};a\in\mathcal{A}_{\mathbf{s}}$
(7)

An example next illustrates these concepts.

Example: For the classical lost sales inventory system [[, cf.]]zipkin2008old
with leadtime $\tau>1$, states are represented as vectors in
$\mathbf{s}\in\mathbb{R}^{\tau}$. In particular, let $(\mathbf{s})_{1}$
represent on hand inventory, and $(\mathbf{s})_{2}$ to $(\mathbf{s})_{\tau}$
the pipeline inventory due to arrive in $1,\ldots,\tau-1$ periods,
respectively. Let $D$ denote per-period demand, then $W:=D$ and
$\mathcal{W}:=\mathbb{R}$. $f(\mathbf{s},a,d)=\mathbf{s}^{\prime}$ corresponds
to the next state when taking action $a$ in state $\mathbf{s}$ and facing
demand $d\in\mathbb{R}$. We have
$(\mathbf{s}^{\prime})_{1}=((\mathbf{s})_{1}-d)^{+}+\mathbf{s}_{2}$. Since $a$
corresponds to the order due to arrive in $\tau$ time periods, we have
$(\mathbf{s}^{\prime})_{\tau}=a$ and for $i\in\\{2,\ldots,\tau-1\\}$ that
$(\mathbf{s}^{\prime})_{i}=(\mathbf{s})_{i+1}$. Therefore,
$f((\mathbf{s}_{1},\ldots,\mathbf{s}_{\tau}),a,d):=((\mathbf{s}_{1}-D)^{+}+\mathbf{s}_{2},\mathbf{s}_{3},\ldots,\mathbf{s}_{\tau},a)$.
The costs $c$ when facing demand $d$ in state $\mathbf{s}$ are
$g(\mathbf{s},a,d)=h(\mathbf{s}_{1}-d)^{+}+p(d-\mathbf{s}_{1})^{+}$, where $h$
are the unit cost of holding inventory and $p$ are the unit cost penalty of
lost sales.

Now, let $\xi$ correspond to a composite random variable containing the random
duration of the trajectory and, for every period of the trajectory, an
independent replication of $W$: I.e.

$\xi=\big{(}T(\xi),\forall t\in\\{1,\ldots,T(\xi)\\}:W_{t}(\xi)\big{)}$

Here, $T(\xi)\sim\operatorname*{Geo}(\alpha)$ denotes the duration of the
trajectory, and $W_{t}(\xi)$ for $t\in\\{1,\ldots,T(\xi)\\}$ the random
variables that will govern state transitions and costs (cf. (6) and (7)) for
period $t$ of the trajectory. The random variables $W_{t}(\xi)$ are
independent, and $W_{t}(\xi)\sim W$. When starting in state
$\mathbf{s}=\mathbf{s}_{1}$ and taking action $a=a_{1}$, the accumulated costs
equal:

$\displaystyle Q_{\pi}(s_{1},a_{1}|\xi):=\sum_{t=1}^{T(\xi)}c_{t},$
$\displaystyle\forall
t\in\\{1,\ldots,T(\xi)\\}:c_{t}=g(\mathbf{s}_{t},a_{t},W_{t}(\xi)),\mathbf{s}_{t+1}=f(\mathbf{s}_{t},a_{t},W_{t}(\xi)),a_{t+1}=\pi(\mathbf{s}_{t+1}).$
(8)

By (5-7), we have
$\mathbf{E}_{\xi}[Q_{\pi}(\mathbf{s},a|\xi)]=q_{\pi}(\mathbf{s},a)$.

####  4.2.2 A bandit problem

For given policy $\pi$ and state $\mathbf{s}$, it is discussed how to obtain
the approximation $\hat{\pi}^{+}(\mathbf{s})$ for $(\pi)^{+}(\mathbf{s})$
using simulation optimization based on (8). To avoid burdensome notation,
dependence on $\pi$ and $\mathbf{s}$ will be suppressed, e.g.
$Q(a|\xi):=Q_{\pi}(\mathbf{s},a|\xi)$ and $q(a):=q_{\pi}(\mathbf{s},a)$. The
approach draws independent samples of $\xi$: denote sample $i$ by $\xi_{i}$.
From these samples, we obtain estimators of $q(a)$ and $q(a)-q(a^{\prime})$.
It will be convenient to denote a set containing the first $n$ such estimators
by

$\displaystyle\mathcal{Z}_{n}(a)$
$\displaystyle:=\big{\\{}Q(a|\xi_{i})\big{|}i\in\\{1,\ldots,n\\}\big{\\}}$ (9)
$\displaystyle\mathcal{Z}_{n}(a,a^{\prime})$
$\displaystyle:=\big{\\{}Q(a|\xi_{i})-Q(a^{\prime}|\xi_{i})\big{|}i\in\\{1,\ldots,n\\}\big{\\}}$
(10)

For any such set $\mathcal{Z}$, denote its mean by $\mu[\mathcal{Z}]$ and the
standard error corresponding to this mean by $\sigma_{\mu}[\mathcal{Z}]$.
E.g., $\mathcal{Z}_{n}(a,a^{\prime})$ yields an estimate
$\mu[\mathcal{Z}_{n}(a,a^{\prime})]\pm\sigma_{\mu}[\mathcal{Z}_{n}(a,a^{\prime})]$
for $q(a|\xi_{i})-q(a^{\prime}|\xi_{i})$.

After computing $Q(a|\xi_{i})$ for a number of samples $\xi_{i}$ for all
actions $a\in\mathcal{A}_{\mathbf{s}}$, it may still be unclear which action
is in $\operatorname*{arg\,min}_{a}q(a)$, but it may already be clear that
some actions are not in the $\operatorname*{arg\,min}$, and we would like to
avoid wasting computational resources for those actions. Optimal allocation of
resources in such settings is a widely studied problem, e.g. in literature on
_bandits_ and on _ranking and selection_ [[, see
e.g.]]berry1985bandit,kim2007recent. An approach is proposed in line with
general ideas in this literature, while accounting for our variance control
scheme.

This approach for approximating
$\hat{\pi}^{+}(\mathbf{s})=\operatorname*{arg\,min}_{a}q(a)$ is governed by
parameters $\underaccent{\bar}{n},\bar{n}\in\mathbb{N}$ with
$1<\underaccent{\bar}{n}\leq\bar{n}$, and by some $\epsilon>0$. For the first
$\underaccent{\bar}{n}$ samples $\xi$, $Q(a|\xi_{i})$ is computed for _all_
$a\in\mathcal{A}_{\mathbf{s}}$. For samples $i>\underaccent{\bar}{n}$,
$Q(a|\xi_{i})$ is only computed for actions $a\in\mathcal{A}_{i}$, where
$\mathcal{A}_{i}$ denotes the actions $a$ for which a (paired) Student’s
t-test indicates that the probability that $(\pi)^{+}(\mathbf{s})=a$ exceeds
$\epsilon$. The algorithm terminates when only a single action remains in
$\mathcal{A}_{i}$, or when $\bar{n}$ samples have been drawn. More precisely,
computing $\hat{\pi}^{+}(\mathbf{s})$ involves the following procedure:

  1. 1.

Initialization: For $i\in\\{1,\ldots,\underaccent{\bar}{n}\\}$, draw an
independent sample $\xi_{i}\sim\xi$, and compute $Q(a,\xi_{i})$ for
$a\in\mathcal{A}_{\mathbf{s}}$. Set $i=\underaccent{\bar}{n}$, and let
$\mathcal{A}_{\underaccent{\bar}{n}}=\mathcal{A}_{\mathbf{s}}$.

  2. 2.

Let
$a^{*}_{i}\in\operatorname*{arg\,min}_{a\in\mathcal{A}_{i}}\mu(\mathcal{Z}_{i}(a))$,
and let $\mathcal{A}_{i+1}$ consist of actions $a\in\mathcal{A}_{i}$ for which

$\displaystyle\mu[\mathcal{Z}_{i}(a,a^{*}_{i})]\leq\sigma_{\mu}[\mathcal{Z}_{i}(a,a^{*}_{i})]\Phi^{-1}(1-\epsilon),$
(11)

  3. 3.

If $\mathcal{A}_{i}$ contains only a single action $a$, or if $i=\bar{n}$, set
$\hat{\pi}^{+}(\mathbf{s})=a^{*}_{i}$ and terminate.

  4. 4.

Otherwise, draw an independent sample $\xi_{i+1}$ and compute $Q(a,\xi_{i+1})$
for $a\in\mathcal{A}_{i+1}$. Set $i\leftarrow i+1$ and go to step $2$.

Based on the sample $\mathcal{Z}_{i}(a,a^{*}_{i})$, our best estimate of
$q(a)-q(a^{*}_{i})$ equals
$\mu[\mathcal{Z}_{i}(a,a^{*}_{i})]\pm\sigma_{\mu}[\mathcal{Z}_{i}(a,a^{*}_{i})]$.
Hence, based on normality assumptions for these statistics in isolation, the
probability that $q(a)-q(a^{*}_{i})<0$ can be approximated as
$1-\Phi(\mu/\sigma_{\mu})$, which does not exceed $\epsilon$ unless $a$
satisfies (11). (Note that $a^{\prime}=a^{*}_{i}$ always satisfies (11).)
Appendix B discusses efficient implementations of this procedure.

###  4.3 Choice and collection of samples

Recall that obtaining $\pi_{i+1}$ from $\pi_{i}$ involves collecting a dataset
$\mathcal{K}_{i}=\\{(\mathbf{s}_{k},\hat{\pi}_{i}^{+}(\mathbf{s}_{k}))|k\in\\{1,\ldots,K\\}\\}$.
$K\in\mathbb{N}$ is a parameter of the algorithm. In §4.2.2 the approach
underlying $\hat{\pi}_{i}^{+}$ was discussed; it remains to discuss the
selection of the states $\mathbf{s}_{k}$.

Samples $(\mathbf{s}_{k},\hat{\pi}_{i}^{+}(\mathbf{s}_{k}))$ are collected in
parallel: Each worker is tasked to collect a number of samples such that work
is evenly distributed and such that in total $K$ samples are collected.
Appropriate selection of states $\mathbf{s}_{k}$ for inclusion in
$\mathcal{K}$ is important. An approach related to the one proposed by [23]
was found to have rather robust performance. In particular, suppose some
(possibly random) starting state $\mathbf{s}\in\mathcal{S}$ is available. The
approach selects states that are encountered when following, in each period,
the simulation-based policy $\hat{\pi}_{i}^{+}(\cdot)$ with probability
$1-\beta$, and taking a random action with probability $\beta$, where
$\beta\in[0,1]$ is an algorithm parameter. More precisely, and supposing
$K^{\prime}$ samples are assigned:

  1. 1.

Initialization: Let $\mathbf{s}_{1}=\mathbf{s}$, $k=1$, and
$\mathcal{K}=\emptyset$.

  2. 2.

Compute $a_{k}=\hat{\pi}_{i}^{+}(\mathbf{s}_{k})$, and add sample
$(\mathbf{s}_{k},a_{k})$ to $\mathcal{K}$.

  3. 3.

If $k=K^{\prime}$, terminate.

  4. 4.

Otherwise, with probability $\beta$ draw action $a$ uniformly from
$\mathcal{A}_{\mathbf{s}_{k}}$, and with probability $1-\beta$ set $a=a_{k}$.
Draw $\mathbf{s}_{k+1}\sim\mathbf{P}(\cdot|\mathbf{s}_{k},a)$, set
$k\leftarrow k+1$, and go to step $2$.

##  5 Experiments

The MCL algorithm is tested for the _classical lost sales inventory system_ ,
as motivated in §2. A brief description of this model is in §4.2, for details
refer to [44].

We consider instances of the testbed proposed by [44], which allow for exact
evaluation. Our main benchmark is the A3C implementation of [14], which report
results for a subset of the lost sales instances considered here. To place the
results in perspective while keeping the discussion focused, in addition the
following benchmarks from [44] are reported:

  * •

The _base-stock policy_ is a reasonable heuristic that might be used in
practice. We will report results for the best base-stock policy.

  * •

The Myopic-2 policy is a relatively complex policy that is the best-performing
benchmark policy in [44].

To _find_ a neural network policy, MCL is applied to the discounted model with
$\alpha=0.975$. In line with [44], solutions are evaluated and compared based
on exact evaluation of the _average_ cost per time unit over an infinite
horizon, by explicitly solving the average-cost Bellman equations. For
implementing the MCL algorithm, the bounds in [44] are used to arrive at a
finite action space $\mathcal{A}$ and allowed actions per state
$\mathcal{A}_{\mathbf{s}}$. To initiate the algorithm, $\pi_{0}$ always
selects the largest order in $\mathcal{A}_{\mathbf{s}}$. (Similar results were
obtained when $\pi_{0}$ selects actions randomly from
$\mathcal{A}_{\mathbf{s}}$.)

We next detail the choice of hyperparameters for our experiments with MCL. We
train four generation of neural networks for each problem instance. Each
neural network has the MLP architecture of §3.2 with three hidden layers with
dimensions $128$, $64$, and $64$, respectively, and is trained using a
minibatch size of $64$. Training uses $K=4000$ samples obtained following the
procedures in §4.2.2 and 4.3 with $\underaccent{\bar}{n}=500$, $\bar{n}=4000$,
$\epsilon=0.02$ and $\beta=0.05$. We report the performance of the best-
performing neural network, which is typically (but not invariably) the
fourth/last generation.

Table 1 reports the optimality cap for the various benchmarks, computed as
$(C(\pi)-C(\pi^{*}))/C(\pi^{*})\times 100\%$, with $C(\pi)$ the average cost
per time unit of the policy and $C(\pi^{*})$ the average cost per time unit of
the optimal policy. Gaps for _best base-stock_ and _myopic-2_ are computed
from [44], while the gaps for A3C are from [14]. [14] did not obtain results
for all problems in the benchmark, perhaps because training the policies with
A3C was rather time-consuming. Gaps of benchmarks are reported with rounding
to one decimal place, following literature. For MCL, rounding aims to retain a
sense of the rough magnitude of the optimality gap.

Optimality gaps for MCL are very small but non-zero, demonstrating that the
learned neural network is able to capture the structure of the optimal policy
to the extent that it can almost match its performance. There does not appear
to be a very clear pattern in optimality gaps.

Learning the neural networks took a bit under $3$ hours for all 24 instances
of Table 1 together, where instances are processed sequentially on $5$ nodes
of a supercomputer (see appendix B for hardware and implementation details).
On a single modern high-performance computer, computation times would be about
2 to 4 times longer. Computation times for individual instances ranged from
$5$ to $11$ minutes, and were mainly dependent on $p$, which makes sense since
computation times are dominated by sample collection and higher values of $p$
imply a higher variance of $Q(s,a|\xi)$, which in turn implies that more
samples are needed to arrive at the desired accuracy.

Some comments on hyperparameter tuning follows. The choice for the presented
experiments was based on earlier experience with related MDP models, with the
objective of adopting a single set of hyperparameters that gives good
performance across all cases. No tuning was involved except for the parameter
$\alpha$, which was initially set at $0.99$: This resulted in similar
performance but at the expense of computation times that are roughly 2 times
slower. In general, for average-cost MDPs in an infinite horizon, $\alpha$
must be chosen close enough to $1$ such that the resulting policy does not
exhibit overly myopic behavior, but not too close to $1$ to avoid needlessly
long computation times. Another key parameter is the minibatch size (Appendix
A), that must be chosen small enough to ensure sufficient generalization of
the neural network, even if this comes at the expense of longer training
times. Results are rather insensitive to the precise choice of other
hyperparameters. When choosing $\underaccent{\bar}{n}$, a consideration may be
that it must be high enough to ensure approximate normality of the estimators
involved in order to avoid premature pruning of actions based on skewed
statistical evidence.

Table 1: Optimality gap for selected methods for instances from the [44] testbed. A $-$ indicates that no result is available for the benchmark.  | Poisson | Geometric
---|---|---
Lead-time $\tau$ |  Lead-time $\tau$
Lost sales penalty | Policy | 2 | 3 | 4 | 2 | 3 | 4
$p=4$ | Best base-stock | 5.5% | 8.2% | 9.9% | 4.5% | 6.4% | 7.8%
Myopic-2 | 0.2% | 0.8% | 1.9% | 0.5% | 1.2% | 1.7%
A3C | 3.2% | 3.0% | 6.7% | $-$ | $-$ | $-$
MCL | 0.0003% | 0.001% | 0.03% | 0.01% | 0.01% | 0.01%
$p=9$ | Best base-stock | 3.7% | 5.1% | 6.4% | 3.1% | 4.6% | 5.8%
Myopic-2 | 0.2% | 0.6% | 1.2% | 0.5% | 1.1% | 1.9%
A3C | 4.8% | 3.1% | 3.4% | $-$ | $-$ | $-$
MCL | 0.001% | 0.004% | 0.02% | 0.01% | 0.01% | 0.01%
$p=19$ | Best base-stock | 2.3% | 2.9% | 3.9% | 2.0% | 3.0% | 3.9%
Myopic-2 | 0.1% | 0.4% | 0.8% | 0.3% | 0.8% | 1.4%
MCL | 0.001% | 0.01% | 0.04% | 0.007% | 0.03% | 0.01%
$p=39$ | Best base-stock | 0.9% | 1.8% | 2.5% | 1.3% | 2.0% | 2.6%
Myopic-2 | 0.1% | 0.3% | 0.4% | 0.2% | 0.5% | 0.9%
MCL | 0.002% | 0.02%% | 0.097% | 0.02% | 0.04% | 0.06%

##  6 Discussion and conclusion

In this paper, the MCL algorithm was proposed and tested. Inspired by [14],
further evidence was provided demonstrating that the potential of neural
network based policies as a _general-purpose solution methodology_ for dynamic
decision problems arising in logistics and supply chain is substantial. In
particular, §5 demonstrated that neural networks can vastly outperform
heuristics that are carefully designed for a specific problem. Given the
generic nature of ideas discussed in this paper, it seems likely that these
results carry over to other relevant and difficult problems in supply chain
and logistics. Thus further research into neural network policies is
appropriate.

Unlike prior applications of neural-network based policies for stochastic
problems in operations management, the MCL algorithm is based on roll-outs and
it incorporates model-based variance reduction techniques. Apart from variance
control, the following properties of MCL may partially explain its apparent
superiority to algorithms like A3C for the focal problems of this paper: 1)
Supervised training with stable targets; 2) Avoids value approximations; 3)
well-suited to efficient hardware implementation [[, cf.]]silver2018general.
It seems to be reasonable to expect that the benefits of MCL compared to
model-free algorithms like A3C will extend to other high-variance problems,
especially when uncertain factors that influence trajectories can be
separately identified and controlled for.

The exposition in this paper focuses on infinite-horizon discounted MDPs with
a discrete, finite action set and a discrete state space. The model-based
learning approach extends to finite-horizon MDPs by defining
$Q_{\pi}(\mathbf{s},a)$ as the costs when taking action $a$ in state
$\mathbf{s}$, and following $\pi(\cdot)$ until the end of the horizon. Finite-
horizon cases can be treated in this fashion for discounted ($\alpha<1$) and
undiscounted ($\alpha=1$) MDPs. The undiscounted _infinite-horizon_ case
cannot be directly handled using the approach, but policies for discounted
MDPs typically transfer well to the corresponding undiscounted case (cf. §5).
The assumption of a discrete state space $\mathcal{S}$ eases the burden of
notation, but it not essential to the approach.

It is essential that the MDP has a discrete and finite action space. _Truly_
continuous action spaces are rare in operations management, so assuming
discrete action spaces may not be a serious limitation in our field. A more
pressing challenge arises when each period involves _multiple_ decisions, e.g.
in multi-SKU inventory problems: The Cartesian product of action spaces may
quickly grow very large, which is problematic since $|\mathcal{A}|=m$ is the
output dimension of the neural network. However, considering that MCL supports
state-dependent action spaces, an appropriate redefinition of the MDP may
still enable effective treatment of such models.

The proposed algorithm controls variance by constructing a random variable
$\xi$ and a deterministic function $Q_{\pi}(\mathbf{s},a|\xi)$ such that
$\mathbf{E}_{\xi}(Q_{\pi}(\mathbf{s},a|\xi))=q_{\pi}(\mathbf{s},a)$. The
concrete approach discussed in this paper assumed the existence of exogenous
random variables $W$ [[, cf.]]bertsekas1997differential that govern state
transitions, and as a consequence $\xi$ corresponds to a sequence of i.i.d.
replications of $W$. The assumption of exogenous and i.i.d. $W$ already allows
for the modelling of numerous problems besides the lost sales inventory
problem that was explored in detail in this paper, e.g. lost sales inventory
problems with a fixed setup [8], multi-echelon inventory control [21], fire-
brigade dispatching [39], etc.

However, the framework also supports variance control when the assumption of
i.i.d. random variables per time period is inappropriate. An extension to
inventory models with time-dependent or correlated demand is relatively
straightforward, once it is noted that the construction of $\xi$ may be state
and time dependent. Some other interesting cases are problems where the true
distribution of random variables is only learned while taking actions, such as
response-adaptive clinical trial design [[]]ahuja2020approximation. It would
be interesting to investigate whether appropriate constructions of $\xi$ exist
for such cases.

## Acknowledgements

This work was carried out on the Dutch national e-infrastructure with the
support of SURF Cooperative, and the author is deeply indebted to the
Cartesius helpdesk team for their invaluable support in resolving various
technical difficulties. The author thanks presenters and participants in the
Eindhoven seminar on reinforcement learning for stimulating discussions and
feedback, and in particular Paulo de Oliveira da Costa for helpful discussions
on supervised training of neural networks.

## References

  * [1] Vishal Ahuja and John R Birge  “An Approximation Approach for Response-Adaptive Clinical Trial Design”  In _INFORMS Journal on Computing_ INFORMS, 2020
  * [2] Donald A Berry and Bert Fristedt  “Bandit problems”  London: ChapmanHall, 1985
  * [3] Dimitri P Bertsekas  “Differential training of rollout policies”, Citeseer, 1997
  * [4] Dimitri P Bertsekas  “Approximate policy iteration: A survey and some new methods”  In _Journal of Control Theory and Applications_ 9.3 Springer, 2011, pp. 310–335
  * [5] Dimitri P Bertsekas  “Feature-based aggregation and deep reinforcement learning: A survey and some new implementations”  In _IEEE/CAA Journal of Automatica Sinica_ 6.1 IEEE, 2018, pp. 1–31
  * [6] Dimitri P Bertsekas  “Reinforcement learning and optimal control”  Athena Scientific Belmont, MA, 2019
  * [7] Sushmita Bhattacharya et al.  “Multiagent Rollout and Policy Iteration for POMDP with Application to Multi-Robot Repair Problems”  In _arXiv preprint arXiv:2011.04222_ , 2020
  * [8] Marco Bijvank, Sandjai Bhulai and Woonghee Tim Huh  “Parametric replenishment policies for inventory systems with lost sales and fixed order cost”  In _European Journal of Operational Research_ 241.2 Elsevier, 2015, pp. 381–390
  * [9] Marco Bijvank and Iris FA Vis  “Lost-sales inventory theory: A review”  In _European Journal of Operational Research_ 215.1 Elsevier, 2011, pp. 1–13
  * [10] Ingeborg A Bikker, Martijn RK Mes, Antoine Sauré and Richard J Boucherie  “Online capacity planning for rehabilitation treatments: An approximate dynamic programming approach”  In _Probability in the Engineering and Informational Sciences_ 34.3 Cambridge University Press, 2020, pp. 381–405
  * [11] Wenbo Chen and Huixiao Yang  “A heuristic based on quadratic approximation for dual sourcing problem with general lead times and supply capacity uncertainty”  In _IISE Transactions_ 51.9 Taylor & Francis, 2019, pp. 943–956
  * [12] Christos Dimitrakakis and Michail G Lagoudakis  “Rollout sampling approximate policy iteration”  In _Machine Learning_ 72.3 Springer, 2008, pp. 157–171
  * [13] Jiarui Fang, Lei Zhao, Jan C Fransoo and Tom Van Woensel  “Sourcing strategies in supply risk management: An approximate dynamic programming approach”  In _Computers & Operations Research_ 40.5 Elsevier, 2013, pp. 1371–1382
  * [14] Joren Gijsbrechts, Robert N Boute, Jan A Van Mieghem and Dennis Zhang  “Can deep reinforcement learning improve inventory management? performance on dual sourcing, lost sales and multi-echelon problems”  In _Performance on Dual Sourcing, Lost Sales and Multi-Echelon Problems (July 29, 2019)_ , 2019
  * [15] David A Goldberg et al.  “Asymptotic optimality of constant-order policies for lost sales inventory models with large lead times”  In _Mathematics of Operations Research_ 41.3 INFORMS, 2016, pp. 898–913
  * [16] Ian Goodfellow, Yoshua Bengio, Aaron Courville and Yoshua Bengio  “Deep learning”  MIT press Cambridge, 2016
  * [17] Woonghee Tim Huh, Ganesh Janakiraman, John A Muckstadt and Paat Rusmevichientong  “Asymptotic optimality of order-up-to policies in lost sales inventory systems”  In _Management Science_ 55.3 INFORMS, 2009, pp. 404–420
  * [18] Phillip R Jenkins, Matthew J Robbins and Brian J Lunday  “Approximate dynamic programming for military medical evacuation dispatching policies”  In _INFORMS Journal on Computing_ INFORMS, 2020
  * [19] Seong-Hee Kim and Barry L Nelson  “Recent advances in ranking and selection”  In _2007 Winter Simulation Conference_ , 2007, pp. 162–172  IEEE
  * [20] Diederik P Kingma and Jimmy Ba  “Adam: A method for stochastic optimization”  In _arXiv preprint arXiv:1412.6980_ , 2014
  * [21] Ton Kok et al.  “A typology and literature review on stochastic multi-echelon inventory models”  In _European Journal of Operational Research_ 269.3 Elsevier, 2018, pp. 955–983
  * [22] D.P. Kroese, Z.I. Botev, T. Taimre and R. Vaisman  “Data Science and Machine Learning: Mathematical and Statistical Methods”, Chapman & Hall/CRC machine learning & pattern recognition  Boca Raton: CRC Press, 2019  URL: https://people.smp.uq.edu.au/DirkKroese/DSML/
  * [23] Michail G Lagoudakis and Ronald Parr  “Reinforcement learning as classification: Leveraging modern classifiers”  In _Proceedings of the 20th International Conference on Machine Learning (ICML-03)_ , 2003, pp. 424–431
  * [24] Joern Meissner and Olga V Senicheva  “Approximate dynamic programming for lateral transshipment problems in multi-location inventory systems”  In _European Journal of Operational Research_ 265.1 Elsevier, 2018, pp. 49–64
  * [25] Martijn RK Mes and Arturo Pérez Rivera  “Approximate dynamic programming by practical examples”  In _Markov Decision Processes in Practice_ Springer, 2017, pp. 63–101
  * [26] Volodymyr Mnih et al.  “Asynchronous methods for deep reinforcement learning”  In _International conference on machine learning_ , 2016, pp. 1928–1937
  * [27] Volodymyr Mnih et al.  “Human-level control through deep reinforcement learning”  In _nature_ 518.7540 Nature Publishing Group, 2015, pp. 529–533
  * [28] Afshin Oroojlooyjadid, MohammadReza Nazari, Lawrence Snyder and Martin Takáč  “A Deep Q-Network for the Beer Game: A Deep Reinforcement Learning algorithm to Solve Inventory Optimization Problems”  In _arXiv preprint arXiv:1708.05924_ , 2017
  * [29] Mohammad Pirhooshyaran and Lawrence V Snyder  “Simultaneous Decision Making for Stochastic Multi-echelon Inventory Optimization with Deep Neural Networks as Decision Makers”  In _arXiv preprint arXiv:2006.05608_ , 2020
  * [30] Warren B Powell  “Approximate Dynamic Programming: Solving the curses of dimensionality”  John Wiley & Sons, 2007
  * [31] Warren B Powell  “A unified framework for stochastic optimization”  In _European Journal of Operational Research_ 275.3 Elsevier, 2019, pp. 795–821
  * [32] Martin L Puterman  “Markov decision processes: discrete stochastic dynamic programming”  John Wiley & Sons, 2014
  * [33] Daniel F Salas and Warren B Powell  “Benchmarking a scalable approximate dynamic programming algorithm for stochastic control of grid-level energy storage”  In _INFORMS Journal on Computing_ 30.1 INFORMS, 2018, pp. 106–123
  * [34] John Schulman et al.  “Proximal policy optimization algorithms”  In _arXiv preprint arXiv:1707.06347_ , 2017
  * [35] Nicola Secomandi  “A rollout policy for the vehicle routing problem with stochastic demands”  In _Operations Research_ 49.5 INFORMS, 2001, pp. 796–802
  * [36] David Silver et al.  “A general reinforcement learning algorithm that masters chess, shogi, and Go through self-play”  In _Science_ 362.6419 American Association for the Advancement of Science, 2018, pp. 1140–1144
  * [37] Peng Sun, Kai Wang and Paul Zipkin  “Quadratic approximation of cost functions in lost sales and perishable inventory control problems”  In _Fuqua School of Business, Duke University, Durham, NC_ , 2014
  * [38] “SurfSara” Accessed: 2020-05-06, https://userinfo.surfsara.nl/systems/cartesius/description, 2020
  * [39] Joep van den Bogaert  “The data-driven fire department” PDEng thesis  Technische Universiteit Eindhoven, 2019
  * [40] Benjamin Van Roy, Dimitri P Bertsekas, Yuchun Lee and John N Tsitsiklis  “A neuro-dynamic programming approach to retailer inventory management”  In _Proceedings of the 36th IEEE Conference on Decision and Control_ 4, 1997, pp. 4052–4057  IEEE
  * [41] Nathalie Vanvuchelen, Joren Gijsbrechts and Robert Boute  “Use of Proximal Policy Optimization for the Joint Replenishment Problem”  In _Computers in Industry_ 119, 2020, pp. 103239
  * [42] Michael A Voelkel, Anna-Lena Sachs and Ulrich W Thonemann  “An aggregation-based approximate dynamic programming approach for the periodic review model with random yield”  In _European Journal of Operational Research_ 281.2 Elsevier, 2020, pp. 286–298
  * [43] Linwei Xin and David A Goldberg  “Optimality gap of constant-order policies decays exponentially in the lead time for lost sales models”  In _Operations Research_ 64.6 INFORMS, 2016, pp. 1556–1565
  * [44] Paul Zipkin  “Old and new methods for lost-sales inventory systems”  In _Operations Research_ 56.5 INFORMS, 2008, pp. 1256–1263

##  Appendix A Training neural network classifiers

The proposed algorithm does not critically depend on specific approaches for
training the neural network. For the experiments, a standard combination of
best practices was adopted from the body of literature [[, see
e.g.]]mnih2015human. For completeness and to make the paper more accessible
for readers without a background in neural networks, a summary follows of how
$\theta_{i}$ is obtained from the samples in $\mathcal{K}_{i}$. First, a
random subset of $95\%$ of the samples in
$\mathcal{K}_{i}=\\{(\mathbf{s}_{k},a_{k}=\hat{\pi}_{i}^{+}(\mathbf{s}_{k}))|k\in\\{1,\ldots,K\\}\\}$
are designated as the training set $\bar{\mathcal{K}}$, while the remaining
samples are the test set $\tilde{\mathcal{K}}$. The loss for a sample is
defined as the cross-entropy loss between 1) the softmax/softargmax of the
output of the neural network for state $\mathbf{s}_{k}$, where actions that
are not applicable for the state are masked out and 2) the “correct” action
$a$. To arrive at a precise definition, remember that
$\mathcal{A}_{\mathbf{s}}\subseteq\\{1,\ldots,m\\}$ and that
$N_{\theta}(\cdot)$ maps states $\mathbf{s}\in\mathcal{S}$ to vectors in
$\mathbb{R}^{m}$. Then the loss for sample $(\mathbf{s},a)$ is defined as:

$\displaystyle\operatorname{loss}(\mathbf{s},a):=\log\frac{\exp[(N_{\theta}(\mathbf{s}))_{a}]}{\sum_{a^{\prime}\in\mathcal{A}_{\mathbf{s}}}\exp[(N_{\theta}(\mathbf{s}))_{a^{\prime}}]}$
(12)

(This loss decreases when the entry $(N_{\theta}(\mathbf{s}))_{a}$
corresponding to the right label $a$ increases and when entries corresponding
to other actions decrease. Thus decreasing the loss typically results in more
samples for which $\pi_{\theta}(\mathbf{s})=a$, see (3).)

The average loss over the samples in the training set is minimized using
stochastic gradient descent: Each _epoch_ , first the samples in the test set
are divided into minibatches of some fixed size. (This minibatch size is a
hyperparameter of the algorithm.) Subsequently, iterating over the
minibatches, the gradient of the average loss of samples in the minibatch with
respect to $\theta$ is computed, and $\theta$ is updated in a direction
opposite to the gradient [[, cf.], using standard parameters]kingma2014adam.

After every $5$ _epochs_ , the loss on the test set is computed, and the
network parameters $\theta$ are check-pointed. The algorithm terminates when
the best test loss does not improve for $20$ consecutive epochs, and outputs
the check-pointed network parameters corresponding to the best test loss until
then.

##  Appendix B Implementation and hardware

The MDL algorithm was implemented in c++17 with PyTorch (c++ front-end) for
neural network training and inference. OpenMPI supports the collection of
samples over multiple nodes on a cluster computer; on each node, the
collection of samples is in parallel over the hardware threads available on
that node. Neural network training is performed on one GPU-enabled node.

The implementation follows closely the overall setup in §4, but it _emulates_
the bandit algorithm in §4.2.2 rather than exactly implementing it: A better
use of hardware is ensured by inferencing $\pi_{\theta}(\mathbf{s})$ for $250$
vectors $\mathbf{s}$ in parallel (on _each_ hardware thread). To achieve this,
the algorithm speculatively initiates trajectories $Q(a,\xi_{i})$ such that
there are at all times $250$ trajectories for different combinations of action
$a$ and sample $\xi_{i}$ for which simulation is ongoing. (The outcome of some
of these initiated trajectories may eventually turn out not to be needed,
depending on the outcome of other pending trajectories.) This _implementation
detail_ still dramatically reduces computation time without affecting the
outcome of the algorithm.

Computations were performed on $1$ GPU node and $4$ thin nodes of the
Cartesius cluster: thin and GPU nodes respectively have Intel E5-2690v3 and
E5-2450v2 processors, while GPU nodes are additionally equipped with two Tesla
K40m GPGPUs [[, cf.]]cartesius.