# Going faster to see further: GPU-accelerated value iteration and simulation
for perishable inventory control using JAX

Joseph Farringtona∗, Kezhi Lia, Wai Keong Wonga,b, Martin Utleyc
aInstitute of Health Informatics, University College London, 222 Euston Road,
London, NW1 2DA, UK
bNIHR University College London Hospitals Biomedical Research Centre,
University College London and University
College London Hospitals NHS Foundation Trust, 149 Tottenham Court Road,
London, W1T 7DN, UK
cClinical Operational Research Unit, University College London, 4 Taviton
Street, London, NW1 2PG, UK

∗Corresponding author<EMAIL_ADDRESS>

###### Abstract

Value iteration can find the optimal replenishment policy for a perishable
inventory problem, but is computationally demanding due to the large state
spaces that are required to represent the age profile of stock. The parallel
processing capabilities of modern GPUs can reduce the wall time required to
run value iteration by updating many states simultaneously. The adoption of
GPU-accelerated approaches has been limited in operational research relative
to other fields like machine learning, in which new software frameworks have
made GPU programming widely accessible. We used the Python library JAX to
implement value iteration and simulators of the underlying Markov decision
processes in a high-level API, and relied on this library’s function
transformations and compiler to efficiently utilize GPU hardware. Our method
can extend use of value iteration to settings that were previously considered
infeasible or impractical. We demonstrate this on example scenarios from three
recent studies which include problems with over 16 million states and
additional problem features, such as substitution between products, that
increase computational complexity. We compare the performance of the optimal
replenishment policies to heuristic policies, fitted using simulation
optimization in JAX which allowed the parallel evaluation of multiple
candidate policy parameters on thousands of simulated years. The heuristic
policies gave a maximum optimality gap of 2.49%. Our general approach may be
applicable to a wide range of problems in operational research that would
benefit from large-scale parallel computation on consumer-grade GPU hardware.

Keywords: value iteration; simulation optimization; reinforcement learning;
perishable inventory;
parallel algorithm

##  1 Introduction

Perishable items, such as fresh food and blood products, “undergo change in
storage so that in time they may become partially or entirely unfit for
consumption” [43]. This means that wastage must be considered alongside the
impact of shortages and stock-holding levels when making replenishment
decisions. Wastage may be a concern for economic, sustainability, or ethical
reasons. One target of the United Nations Sustainable Development Goals is to
reduce wastage at the retail and consumer levels of the food supply chain by
half, from an estimated 17% of total food production [58]. In the blood supply
chain, platelets can only be stored for between three and seven days leading
to high reported wastage rates of 10-20% [20]. Better policies for perishable
inventory control could help to reduce wastage and make the best possible use
of limited resources.

Early theoretical work demonstrated that optimal policies for perishable
inventory replenishment could be found using dynamic programming [42, 22].
Value iteration is a dynamic programming approach that can be used to find the
optimal policy when the problem is framed as a Markov decision process [7].
The optimal policy depends on the age profile of the inventory, not just on
the total number of units in stock. This approach is therefore limited by the
“curse of dimensionality”: the computational requirements grow exponentially
with the maximum useful life of the product [44]. [43] observed that, at the
time, this made dynamic programming approaches impractical for problems where
the maximum useful life of the product was more than two periods. More
recently, despite advances in computational power, researchers have stated
that value iteration remains infeasible or impractical when the maximum useful
life of the product is longer than two or three periods or when additional
complexities (e.g. substitution between products or random remaining useful
life on arrival) are introduced [15, 27, 40]

The prevalent view in the literature about the scale of problems for which
value iteration is feasible appears to neglect the recent developments in
graphical processing unit (GPU) hardware, and in software libraries that make
it possible for researchers and practitioners to take advantage of GPU
capabilities without detailed knowledge of GPU hardware and GPU-specific
programming approaches. GPUs were developed for rendering computer graphics,
which requires the same operations to be efficiently applied to many inputs in
parallel. Compared to a central processing unit (CPU), GPUs therefore have
many more, albeit individually less powerful, cores. GPU-acceleration refers
to offloading computationally intensive tasks that benefit from large-scale
parallelization from the CPU to the GPU. The first mainstream software
framework to support general computing tasks on GPUs using a common general
purpose programming language was Nvidia’s CUDA platform, which launched in
2007. One of the areas in which GPUs have since had a major impact is the
field of deep learning. This impact has led to, and in turn been supported by
[30], the development of higher-level software libraries including TensorFlow
[2], PyTorch [46] and JAX [10]. These libraries provide comparatively simple
Python application programming interfaces (APIs) to support easy
experimentation and developer productivity, while utilising highly optimized
CUDA code “under-the-hood” to exploit the parallel processing capabilities of
GPUs.

In this work, we implemented value iteration using JAX. JAX provides
composable transformations of Python functions, which make it easy to apply
functions in parallel over arrays of input data, and performs just-in-time
compilation (JIT) to run workloads efficiently on hardware accelerators
including GPUs. This is ideal for value iteration: each iteration requires
many independent updates which can be performed in parallel, and the up-front
computational cost of JIT can be amortised over many repeats of the compiled
operation for each iteration. By decreasing the wall time required to run
value iteration, we increase the size of problems for which the optimal policy
can be calculated in practice: by going faster, we can see further. These
policies may themselves be used to guide decision making. They can also
support research into new heuristics and approximate approaches, including
reinforcement learning, by providing performance benchmarks on much larger
problems than has previously been possible.

Inspired by the recent development of GPU-based simulators in reinforcement
learning [21, 39, 37, 9], we also implemented simulators for perishable
inventory problems using JAX, which enabled us to run large numbers of
simulations in parallel.

We consider perishable inventory scenarios from three recent studies, where
running value iteration for certain settings was described as computationally
infeasible or impractical. For most of these settings, we have been able to
find the optimal policy using a consumer-grade GPU and report the wall time
required for each experiment. We compare the performance of the policies found
using value iteration with the performance of heuristic policies with
parameters fitted using simulation optimization.

The main contributions of this work are in demonstrating that:

  * •

value iteration can be used to find optimal policies, for scenarios for which
it has recently been described as computationally infeasible or impractical,
on consumer-grade GPU hardware (summarised in Table 1);

  * •

this performance can be achieved without in-depth knowledge of GPU-specific
programming approaches and frameworks, using the Python library JAX;

  * •

simulation optimization for perishable inventory control can also be
effectively run in parallel using JAX, particularly for larger problems where
policies that perform well can be identified in a fraction of the time
required to run value iteration.

For one large problem with over 16 million states, a CPU-based MATLAB
implementation of value iteration did not converge within a week in a prior
study. Using our method, value iteration converges in under 3.5 hours on a
consumer-grade GPU and, without any code changes, in less than 30 minutes
using four data-centre grade GPUs. Our simulation optimization method is able
to evaluate 50 possible sets of parameters for a heuristic policy, each on
4,000 years of simulated data, in parallel in under 15 seconds. The largest
optimality gap we observed for the heuristic policies fit using our simulation
optimization method was 2.49%.

Our open-source code is available at
https://github.com/joefarrington/viso_jax. The repository includes a Google
Colab notebook that enables interested readers to reproduce our experiments
using free cloud-based GPUs provided by Google.

|  Maximum useful life $m$
---|---
| Problem features |  2 |  3 |  4 |  5 |  8
A |  Lead time $>$ 1 |  ⚫ |  ✱ |  ✱ |  ✱ |  ❍
|  |  |  |  |  |
B | Substitution between products |  ⚫ |  ✱ |  ❍ |  ❍ |  ❍
|  |  |  |  |  |
C | Not all arrivals fresh, periodic demand |  ❍ |  ⚫ |  ❍ |  ✱ |  ✖

| Key
---|---
⚫ | Value iteration feasible for all experiments in the original study.
✱ | Our method extends value iteration to experiments that were considered
| infeasible or impractical in the original study.
✖ | All experiments infeasible in the original study and with our method.
❍ | Setting not considered in the original study.

Table 1: Summary of our contribution extending value iteration to larger
problems with a longer
maximum useful life.

##  2 Related work

[31] demonstrated that value iteration could be effectively run in parallel on
a GPU soon after the introduction of CUDA. Subsequent research has evaluated
the performance of GPU-accelerated value iteration on problems from economics
and finance [5, 1, 18, 34, 35] and route-finding and navigation [12, 29, 50,
13]. We have only identified a single study that applied this approach to an
inventory control problem: [45] implemented a custom value iteration algorithm
in CUDA to find replenishment policies for a subset of perishable inventory
problems originally described by [27].

Previous studies have reported impressive reductions in wall time achieved by
running value iteration on GPU. The GPU-accelerated method in [45] was up to
11.7$\times$ faster than a sequential CPU-based method written in C. Despite
this, the approach has not been widely adopted. The majority of the
aforementioned studies focus on implementing value iteration using CUDA or
OpenCL, a multi-platform alternative to CUDA, and comparing the performance of
a GPU-accelerated method with a CPU-based method, instead seeking to use GPU-
acceleration to solve problems for which value iteration is otherwise
impractical or infeasible. One of the main barriers to entry for other
researchers may be the perceived difficulty of GPU programming. Writing
efficient code using CUDA or OpenCL requires careful consideration of memory
access, balancing resource usage when mapping parallel processes to the
hardware, and interaction between the CPU and GPU [28]. There have been
efforts to make GPU-accelerated value iteration more accessible. [31] created
a solver framework using his CUDA implementation of value iteration as a back-
end, but this does not appear to be publicly available. More recently, [34]
created a toolkit in MATLAB to solve infinite horizon value iteration problems
which automatically uses a GPU when available and appropriate. This toolkit
requires the state-transition matrix to be provided as an input, which is not
practical for some of the larger problems we consider due to memory
limitations.

Our approach is broadly similar to that of [18] and [51] who used machine
learning frameworks to implement GPU-accelerated value iteration for economics
models. A key observation made by [18] is that their TensorFlow implementation
is an order of magnitude faster than their custom CUDA C++ implementation on a
GPU. This demonstrates that the comparative accessibility provided by a
machine learning framework need not come at the cost of poorer performance.
TensorFlow, like PyTorch and JAX, translates the high level API instructions
into highly-optimized CUDA code. Experts in GPU programming may be able
achieve better performance than a machine learning framework by working at the
level of CUDA or OpenCL. Researchers without that expertise are likely to both
save development time and achieve performance benefits by working at a higher
level of abstraction using a machine learning framework and relying on it to
make best use of the available hardware.

Due to the computational challenges of using dynamic programming methods to
find policies for perishable inventory management, research has focused on
approximate solutions. One straightforward way to make a dynamic programming
approach more computationally tractable is to reduce the size of the state
space by aggregating stock items into batches. The solution to the down-sized
dynamic program can then by factored up to give an approximate solution to the
the original problem (e.g. [8, 26]). An alternative approach is to use a
heuristic policy with a small number of parameters, such as a base-stock
policy. Research in this area has concentrated on both identifying suitable
structures for heuristic policies (see [41] for an early example, and [25] for
a recent example), and finding suitable parameters for those policies in
specific situations - commonly using stochastic mixed integer linear
programming (e.g. [16, 24, 49]) or simulation optimization (e.g. [14, 17]).
Recently, reinforcement learning methods have also been used to find
approximate polices for managing perishable inventory [32, 55, 15, 3].

Of these other approaches to the problem, we focus on simulation optimization
in addition to value iteration because GPU-accelerated simulation is feasible
using available software libraries but not yet widely adopted. Applied
research on specific mixed integer linear programs often relies on commercial
solver software such as IBM ILOG CPLEX Optimization Studio or Gurobi Optimizer
which do not currently support GPU-acceleration. Adapting mixed integer
programming solution strategies to suit the architecture of GPUs is the
subject of active research [48]. The use of GPU-acceleration supported by
machine learning frameworks is already widespread in reinforcement learning
research.

Simulation optimization can be used to solve optimization problems where the
objective function cannot be computed exactly, but can be estimated using
simulation. Sampling error in the objective function can be reduced by running
simulations for a longer period, or by running additional simulations. The
relevance of parallel computing to simulation optimization is well recognised
[6, 23], but we have identified few simulation optimization studies using GPUs
to run multiple simulations in parallel. In inventory management, [53]
exhaustively evaluated the possible order quantities for a newsvendor problem
using parallel simulations on GPU. More recently, [38] used simulation
optimization to solve a chemical process monitoring problem using GPU-
acceleration for both the simulation and the metaheuristic search process that
proposed candidate solutions. Similar to the value iteration studies discussed
above, both of these projects used custom CUDA code which may explain the
limited subsequent adoption despite the established reductions in wall time
relative to CPU baselines. Following recent work in the reinforcement learning
community [21, 36] we implemented our simulators using the Python library
gymnax [37], which enabled us to write our simulation operations using JAX and
readily evaluate each of numerous policies on thousands of simulated years in
parallel on GPU.

##  3 Methods

###  3.1 Scenarios

We considered three scenarios, all of which are periodic review, single-
echelon perishable inventory problems with a fixed, known delivery lead time
$L$. The three scenarios were selected as recent examples from the perishable
inventory literature in which value iteration was reported as infeasible or
impractical for at least some experimental settings, and which include
elements relevant to our wider work investigating the potential of
reinforcement learning methods to support blood product inventory management.
Scenario A, from [15], is a straightforward perishable inventory replenishment
problem but for some experimental settings the lead time, $L$, is greater than
one period and therefore we need to consider inventory in transit when placing
an order. Scenario B is the two product scenario described by [27] which adds
the complexity of substitution between perishable products. Substitution is an
important aspect of managing blood product inventory, where compatibility
between the blood groups of the donor and the recipient is critical. Scenario
C, from [40], models the management of platelets in a hospital blood bank and
adds two complicating factors: periodic patterns of demand, and uncertainty in
the remaining useful life of products on arrival, which may depend on the
order quantity. In every scenario demand is stochastic, unmet demand is
assumed to be lost, and units in stock with a remaining useful life of one
period are assumed to expire at the end of the day. Except in Scenario C, the
products have a fixed, known useful life $m$ and are all assumed to arrive
fresh. We summarise the key differences between the scenarios in Table 2.

For readers who are unfamiliar with inventory management problems and the
associated terminology we recommend Chapters 3 and 4 of [52] for a general
introduction and [11] and [44] for more focused coverage of perishable
inventory control.

| Problem features |  | Reward function components |
---|---|---|---|---
| Source |  Products |  Lead time $>1$ |  Substitution |  Not all arrivals fresh |  Periodic demand |  |  |  |  Variable ordering |  Fixed ordering |  Wastage |  Shortage |  Holding |  Revenue |
A | [15] |  1 |  ✓ |  |  |  |  |  |  |  ✓ |  |  ✓ |  ✓ |  ✓ |  |
|  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
B | [27] |  2 |  |  ✓ |  |  |  |  |  |  ✓ |  |  |  |  |  ✓ |
|  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
C | [40] |  1 |  |  |  ✓ |  ✓ |  |  |  |  |  ✓ |  ✓ |  ✓ |  ✓ |  |
Table 2: Summary of the key differences between our three scenarios.

[15], [27] and [40] each used different notation to describe their work. In an
effort to aid the reader in understanding the similarities and differences
between the scenarios we have adopted a single notation which we apply to all
three scenarios. In Appendices A.1, B.1 and C.1 we present the key equations
describing each scenario in our notation and provide a table summarising our
notation in Appendix D.

All of the scenarios are defined as Markov decision processes (MDPs). An MDP
is a formal description of a sequential decision problem in which, at a
discrete series of points in time, an agent observes the state of its
environment $S_{t}$ and selects an action $A_{t}$. At the next point in time,
the agent will receive a reward signal $R_{t+1}$, observe the updated state of
its environment $S_{t+1}$ and must select its next action $A_{t+1}$. An MDP
can be defined in terms of a set of states $s\in\mathbb{S}$, a set of actions
$a\in\mathbb{A}$, a set of a rewards $r\in\mathbb{\Psi}$, a function defining
the dynamics of the MDP (Equation 1), and a discount factor $\gamma\in[0,1]$
[56]. The discount factor controls the relative contribution of future rewards
and immediate rewards. The decision process is Markovian because the dynamics
of the system obey the Markov property: state transitions and rewards at time
$t$ are conditionally independent of the sequence of state-action pairs
$(S_{0},A_{0})$ to $(S_{t-1},A_{t-1})$ given $(S_{t},A_{t})$.

$p(s^{\prime},r|s,a)=\text{Prob}\left(S_{t}=s^{\prime},R_{t}=r|S_{t-1}=s,A_{t-1}=a\right)$
(1)

Within this framework, MDP agents select their actions by following a policy.
In this work, we only consider deterministic policies, $a=\pi(s)$. The
objective is to find a policy that maximises the expected return, the
discounted sum of future rewards, when interacting with the environment. In an
infinite horizon problem the return at timestep $t$ is
$G_{t}=\sum_{k=0}^{\infty}\gamma^{k}R_{t+k+1}$.

###  3.2 Value iteration

Following the treatment of [56], the value of a state under a policy $\pi$,
$V^{\pi}(s)$, is the expected return when starting in state $s$ and following
policy $\pi$. For a finite MDP, one in which the sets of states, actions and
rewards are finite, we can define an optimal policy, $\pi^{*}$, as a policy
for which $V^{\pi^{*}}(s)\geq V^{\pi^{\prime}}(s)$ for every state
$s\in\mathbb{S}$, for any policy $\pi^{\prime}$. There may be more than one
optimal policy, but they all share the same optimal value function. Value
functions satisfy recursive relationships, called Bellman equations, between
the value at the current state and the immediate reward plus the discounted
value at the next state. The Bellman equation for the optimal policy, the
Bellman optimality equation, is:

$V^{\pi^{*}}(s)=\max_{a\in\mathbb{A}}\sum_{s^{\prime}\in\mathbb{S},r\in\mathbb{\Psi}}p(s^{\prime},r|s,a)\left[r+\gamma
V^{\pi^{*}}(s^{\prime})\right]$ (2)

Value iteration is a dynamic programming algorithm, which uses the Bellman
optimality equation as an update operation to estimate the optimal value
function:

$V_{i+1}(s)=\max_{a\in\mathbb{A}}\sum_{s^{\prime}\in\mathbb{S},r\in\mathbb{\Psi}}p(s^{\prime},r|s,a)\left[r+\gamma
V_{i}(s^{\prime})\right]$ (3)

For a finite MDP this operation will, in the limit of infinite iterations,
converge to the optimal value function. The optimal policy can be extracted
from the value function using a one-step ahead search:

$\pi(s)=\operatorname*{arg\,max}_{a\in\mathbb{A}}\sum_{s^{\prime}\in\mathbb{S},r\in\mathbb{\Psi}}p(s^{\prime},r|s,a)\left[r+\gamma
V(s^{\prime})\right]$ (4)

Similar to the approach of [27], we used a deterministic transition function,
$(s^{\prime},r)=T(s,a,\omega)$, where $\omega\in\mathbb{\Omega}$ is a possible
realisation of the stochastic element(s) of the transition. For a specific
state-action pair, $(s,a)$, and a specific random outcome, $\omega$, the next
state and the reward can be calculated deterministically. In the most
straightforward example we consider, Scenario A, the only uncertainty is in
the daily demand, and therefore $\mathbb{\Omega}$ is the set of possible
values that demand may take in any period. Under this formulation, the value
iteration update equation can be rewritten as:

$V_{i+1}(s)=\max_{a\in\mathbb{A}}\sum_{\omega\in\mathbb{\Omega}}P(\omega|s,a)\left[r_{\omega}+\gamma
V_{i}(s^{\prime}_{\omega})\right],\text{where
}(r_{\omega},s^{\prime}_{\omega})=T(s,a,\omega)$ (5)

and we extract the optimal policy using the equation:

$\pi(s)=\operatorname*{arg\,max}_{a\in\mathbb{A}}\sum_{\omega\in\mathbb{\Omega}}P(\omega|s,a)\left[r_{\omega}+\gamma
V(s^{\prime}_{\omega})\right],\text{where
}(r_{\omega},s^{\prime}_{\omega})=T(s,a,\omega)$ (6)

where
$P(\omega|s,a)=\text{Prob}\left(\Omega_{t}=\omega|S_{t}=s,A_{t}=a\right)$ is
the probability of random outcome $\omega$ having observed state $S_{t}=s$ and
then taken action $A_{t}=a$. $\Omega_{t}$ represents the stochastic elements
of the transition that occur between the observation of state $S_{t}$ and the
observation of state $S_{t+1}$.

Since we cannot run an infinite number of iterations, we use a convergence
test to determine when to stop value iteration and extract the policy. We are
interested in the policy, and not the value function itself, and therefore in
certain cases we can reduce the number of iterations required by stopping when
further updates to the value function will not change the policy. We describe
the convergence test used for each scenario in the corresponding section
below.

We implemented value iteration using a custom Python class, VIRunner, which
defines the common functionality required to run value iteration and extract
the optimal policy using the approach in Equations 5 and 6. The base VIRunner
class includes eight placeholders for methods which must be defined for a
specific scenario. For each scenario we defined a subclass of VIRunner,
replacing the placeholder methods with custom functions that:

  * •

return a list of all possible states as tuples;

  * •

return an array that maps from a state to its index in the list of all
possible states;

  * •

return an array of all possible actions;

  * •

return an array of all possible random outcomes;

  * •

return the immediate reward and next state following the deterministic
transition function given a state, action and random outcome;

  * •

return an array with the probability of each random outcome given a state-
action pair;

  * •

return an initial estimate of the value function; and

  * •

test for convergence of the value iteration procedure.

The VIRunner class could be easily adapted to solve new problems by creating a
new subclass and replacing the placeholder entries for these eight methods.

A naive implementation of value iteration following Equation 5 would require a
nested for-loop over every state, every action, and every random outcome. In a
single iteration, these updates are independent and therefore can be performed
in parallel. JAX provides two main composable function transformations that
facilitate running functions in parallel: vectorizing map (vmap) and parallel
map (pmap). Both of these transformations create new functions that map the
original function over specified axes of the input, enabling the original
function to applied to a large number of inputs in parallel. The key
difference is that vmap provides vectorization, so the operations happen on
the same device (e.g. the same GPU), while pmap supports single-program,
multiple-data parallelism and runs the operation for different inputs on
separate (but identical) devices.

An important feature of vmap and pmap is that they are composable and
therefore can be readily nested. For our basic value iteration update, set out
in Algorithm 1, we nest vmap operations over states, actions, and random
outcomes instead of using nested loops. This is only feasible if there is
sufficient GPU memory to update all of the states simultaneously. For larger
instances we grouped the states into batches for which the update can be
performed simultaneously and performed the update for one batch of states at a
time. To enable multiple identical devices to be used where available, we
automatically detected the the number of available of devices and used pmap to
map our update function for multiple batches of states over the available
devices. Each batch must contain the same number of states, and each device
must receive the same number of batches, to efficiently loop over batches of
states and use pmap. We therefore padded the array of states so that it could
be reshaped to an array with dimensions (number of devices, number of batches,
maximum batch size, number of elements in state). Each device received an
array with dimensions (number of batches, maximum batch size, number of
elements in state), performed a loop over the leading dimension, and
calculated the update one batch of states at a time. The same process was used
to extract the policy in parallel at the end of value iteration.

Algorithm 1 Value iteration using vmap

Initialise array of all states $s$: $\mathbb{S}$

Initialise array of all actions $a$: $\mathbb{A}$

Initialise array of all random outcomes $\omega$: $\mathbb{\Omega}$

Initialise initial estimate of value function: $V_{0}(s)\quad\forall
s\in\mathbb{S}$

Initialise discount factor: $\gamma$

Initialise iteration counter: $i=0$

Define deterministic transition function which returns next state and reward:
$T(s,a,\omega)$

Perform value iteration

while not converged do

$i\leftarrow i+1$

vmap over $s\in\mathbb{S}$

vmap over $a\in\mathbb{A}$

vmap over $\omega\in\mathbb{\Omega}$

$(s^{\prime}_{\omega},r_{\omega})\leftarrow T(s,a,\omega)$

$Q_{i}(s,a)\leftarrow\sum_{\omega}P(\omega|s,a)\left[r_{\omega}+\gamma
V_{i-1}(s^{\prime}_{\omega})\right]$

$V_{i}(s)\leftarrow\max_{a}Q_{i}(s,a)$

Test for convergence

end while

Extract the policy, $\pi(s)\approx\pi^{*}(s)$

vmap over $s\in\mathbb{S}$

vmap over $a\in\mathbb{A}$

vmap over $\omega\in\mathbb{\Omega}$

$(s^{\prime}_{\omega},r_{\omega})\leftarrow T(s,a,\omega)$

$Q_{i+1}(s,a)\leftarrow\sum_{\omega}P(\omega|s,a)\left[r_{\omega}+\gamma
V_{i}(s^{\prime}_{\omega})\right]$

$\pi(s)\leftarrow\operatorname*{arg\,max}_{a}Q_{i+1}(s,a)$

Functions transformed by pmap are automatically JIT compiled with XLA, a
domain-specific compiler for linear algebra. JAX traces the function the first
time it is run, and the traced function is compiled using XLA into optimized
code for the available devices.

We used double-precision (64-bit) numbers when running value iteration,
instead of the single-precision (32-bit) numbers that JAX uses by default
because, in preliminary experiments, we found that convergence was not always
stable.

We report the wall time required to run each value iteration experiment. The
reported times include JIT compilation time, writing checkpoints and writing
final outputs, including the policy, because we believe this represents a
realistic use case.

###  3.3 Simulation of the Markov decision processes

We created a simulator to represent each scenario. The simulators have two
purposes: firstly, to fit parameters for heuristic replenishment policies
using simulation optimization and, secondly, to evaluate the performance of
the policies produced by value iteration and simulation optimization based on
the return and three key performance indicators (KPIs): service level, wastage
and holding. The service level is the percentage of demand that was met over a
simulated rollout, wastage is the proportion of units received that expired
over a simulated rollout and holding is the mean number of units in stock at
the end of each day during a simulated rollout.

Each simulator is a reinforcement learning environment written using the
Python library gymnax [37], which is based on JAX. This provides a standard
interface for working with MDPs, while enabling many simulations to be run in
parallel on a GPU using vmap. For a single policy, we can use vmap to
implement our simulation rollout over multiple random seeds. We can
simultaneously evaluate multiple sets of parameters for the same policy on a
shared set of random seeds by nesting vmapped functions. We note that it would
be straightforward to use reinforcement learning software libraries to learn
policies for these scenarios using these environments.

We selected different heuristic policies for the different scenarios from the
literature, considering which (if any) heuristic was used in the original
study and the structure of each problem. We describe the heuristic policy used
for each scenario in the corresponding section below. All of the heuristic
policies use one or both of an order-up-to level parameter S and reorder point
parameter s. The order quantity is the difference between the current stock on
hand and in transit (potentially subject to some modification, as in Scenario
B) and the order-up-to level S. If there were no stock on hand or in transit
the heuristic policy would order S units, and therefore S corresponds to the
largest order that would be placed following the heuristic policy. If the
heuristic policy also has a reorder point parameter s then an order is only
placed when the current stock on hand and in transit is less than or equal to
the reorder point s [52].

We used the Python library Optuna [4] to suggest parameters for the heuristic
policies. For heuristic policies with a single parameter, we evaluated all
feasible values simultaneously in parallel using Optuna’s grid sampler. When
there was more than one parameter we instead used a genetic algorithm,
Optuna’s NSGAII sampler, to search the parameter space. For each suggested set
of parameters, we ran 4,000 rollouts, each 365 days long following a warm-up
period of 100 days. When using the grid search sampler, we took as the best
parameter value the one with the highest mean return after the single parallel
run. When using the NSGAII sampler we ran 50 sets of parameters in parallel,
representing a single generation for the genetic algorithm, and ranked them
based on the mean return. We terminated the NSGAII search procedure when the
best combination of parameters had not changed for five generations, or when
100 generations had been completed.

For each scenario we compare the performance of the value iteration policy and
best heuristic policy identified using simulation optimization on 10,000
simulated rollouts, each 365 days long following a warm-up period of 100 days.
We report the mean and standard deviation of the return and, in the
appendices, the service level, wastage and stock holding over these rollouts.
For each rollout, the return is the discounted sum of rewards from the end of
the warm-up period until the end of the simulation. The components of the
reward function for Scenarios A, B and C are summarised in Table 2 and the
reward functions are set out in Appendices A.1, B.1 and C.1 respectively. The
standard deviation of the return and the KPIs shows the effect of the
stochasticity (due to random demand, random willingness to accept substitution
and/or random useful life on arrival) in each scenario.

###  3.4 Reproducibility

There are two key reproducibility considerations for this work: firstly,
accurately implementing the scenarios described in previous studies and,
secondly, ensuring that others are able to reproduce our own experiments.

We compared outputs from our value iteration and simulation optimization
methods to outputs from the original studies, and these checks are included as
automated tests in our publicly available GitHub repository. [15] made their
code available on GitHub and fully specified the optimal and heuristic
policies for two experiments in their paper which we used to test our
implementation of Scenario A. For Scenario B, we compared the best parameters
for heuristic policies and mean daily reward values to those reported in [27],
and performed additional comparisons to the output of a MATLAB implementation
of their value iteration method that the authors kindly made available to us.
[40] plotted value iteration policies for a subset of his experiments, and he
kindly provided us with the underlying data for those plots so that we could
confirm the policies from our implementation of Scenario C matched those he
had reported.

Our code is available on GitHub, and is based on open-source software
libraries. Our GitHub repository includes a Google Colab notebook that can be
used to reproduce our experiments using a free, cloud-based GPU, avoiding
local hardware requirements or configuration challenges. The type of GPU
allocated to a session in Colab is not guaranteed and there are service limits
that restrict the maximum continuous running time. Experiments may be
restarted from a checkpoint if a session terminates before the experiment is
completed.

###  3.5 Hardware

All experiments were conducted on a desktop computer running Ubuntu 20.04 LTS
via Windows Subsystem for Linux on Windows 11 with an AMD Ryzen 9 5900X
processor, 64GB RAM, and an Nvidia GeForce RTX 3060 GPU. The Nvidia GeForce
RTX 3060 is a consumer-grade GPU that, at the time of writing in February
2023, can be purchased for less than £400 in the United Kingdom [19].

To demonstrate the potential benefits of more powerful data-centre grade GPU
devices, and how our approach can be easily scaled to utilise multiple GPUs,
we additionally ran value iteration for one large problem case of Scenario B
using one, two or four Nvidia A100 40GB GPUs.

##  4 Scenario A: lead time may be greater than one period

###  4.1 Problem description

[15] described a single-product, single-echelon, periodic review perishable
inventory replenishment problem and investigated whether using heuristic
replenishment policies to shape the reward function can improve the
performance of reinforcement learning methods.

At the start of each day $t$ the agent observes the state $S_{t}$, the current
inventory in stock (split by remaining useful life) and in transit (split by
period ordered), and places a replenishment order
$A_{t}\in\\{0,1,...,A_{\max}\\}$. Demand for day $t$, $D_{t}$, is sampled from
a truncated gamma distribution and rounded to the nearest integer. Demand is
filled from available stock following either a first-in first-out (FIFO) or
last-in first-out (LIFO) issuing policy. At the end of the day, the state is
updated to reflect the ageing of stock and the reward, $R_{t+1}$, is
calculated. The reward function comprises four components: a holding cost per
unit in stock at the end of the period ($C_{h}$), a variable ordering cost per
unit ($C_{v}$), a shortage cost per unit of unmet demand ($C_{s}$) and a
wastage cost per unit that perishes at the end of the period ($C_{w}$). The
order placed on day $t-(L-1)$ is received immediately prior to the start of
day $t+1$, and is included in the stock element of the state $S_{t+1}$.

The stochastic element in the transition is the daily demand $D$,
$\mathbb{\Omega}=\\{0,1,...,\infty\\}$, in the problem described by [15]. The
state transition and the reward are deterministic given a state-action pair
and the realisation of the daily demand. Daily demand is modelled by a gamma
distribution with mean $\mu$ and coefficient of variation
$\frac{\mu}{\sigma}$. We truncated the demand distribution at
$D_{\max}\gg\mu+5\sigma$, such that $\mathbb{\Omega}=\\{0,1,...,D_{\max}\\}$,
for the purposes of implementation.

The initial value function $V_{0}(s)$ was initialised at zero for every state.
[15] did not specify a particular convergence test for their value iteration
experiments. The problem is not periodic and includes a discount factor, and
we therefore we used a standard convergence test for the value function [56]
as set out in Appendix A.1.

[15] considered products with a maximum useful life $m$ of two, three, four or
five periods, and evaluated eight different experimental settings for each
value of $m$. For a product with $m=2$, they found the optimal policy using
value iteration, and used this as a benchmark for their deep reinforcement
learning policies. For larger values of $m$, they instead used a heuristic
policy as the benchmark on grounds of computational feasibility. The
experiments for each value of $m$ evaluate different combinations of lead time
$L$, wastage cost $C_{w}$, and issuing policy. We demonstrate that, using JAX
and a consumer-grade GPU, it is feasible to obtain the optimal policy for all
of the experimental settings, up to and including a maximum useful life $m$ of
five periods, and report the wall time required to run value iteration for
each experiment.

We compare the policy from value iteration with a standard base-stock policy,
parameterised by order-up-to level S, such that the order quantity on day $t$,
given total current stock (on hand and in transit) $I_{t}$ is:

$A_{t}=\left[\texttt{S}-I_{t}\right]^{+}$ (7)

We evaluated the mean return for each value of
$\texttt{S}\in\\{0,...,A_{\max}\\}$ using the Optuna grid sampler. We compare
the base-stock policy that achieves the highest mean return, characterised by
parameter S${}_{\text{best}}$, to the value iteration policy.

See Appendix A.1 for additional information about Scenario A.

###  4.2 Results

In Table 3 we present the wall time (WT) in seconds required to run value
iteration and simulation optimization for each experimental setting. We also
present the mean and standard deviation of the return obtained when using
value iteration and best heuristic policies on 10,000 simulated rollouts, each
365 days long following a warm-up period of 100 days.

The wall times reported in Table 3 show that, using our approach, the largest
cases, with $m=5$ and $L=2$, can be solved using value iteration in under 20
minutes. The running time for simulation optimization is approximately
constant at two seconds over the different problem sizes. This is consistent
with the fact that the parameter space for the heuristic base-stock policy is
the same for each experimental setting, because the maximum order quantity
$A_{\max}$ does not change.

As we would expect, we observe higher mean returns under a FIFO issuing policy
than under a LIFO issuing policy and as $m$ increases due to lower wastage.
The optimality gap is consistently higher for experiments with a longer lead
time, suggesting that the age profile of the stock is more important when lead
times are longer.

See Appendix A.2 for the best parameters for the heuristic policy and KPIs for
each experiment.

|  |  | Value | Simulation
---|---|---|---|---
|  |  | iteration | optimization
$m$ | Exp | $L$ | $C_{w}$ | Issuing | $|\mathbb{S}|$ | $|\mathbb{A}|$ | $|\mathbb{\Omega}|$ | WT (s) | Return | WT (s) | Return | Optimality
|  |  |  | policy |  |  |  |  |  |  |  | gap (%)
2 | 1 | 1 | 7 | LIFO | 121 | 11 | 101 | 5 |  -1,553 $\pm$ 61 | 2 |  -1,565 $\pm$ 62 | 0.80
| 2 | 1 | 7 | FIFO | 121 | 11 | 101 | 4 |  -1,457 $\pm$ 59 | 2 |  -1,474 $\pm$ 56 | 1.20
| 3 | 1 | 10 | LIFO | 121 | 11 | 101 | 5 |  -1,571 $\pm$ 61 | 2 |  -1,581 $\pm$ 62 | 0.64
| 4 | 1 | 10 | FIFO | 121 | 11 | 101 | 5 |  -1,463 $\pm$ 60 | 2 |  -1,485 $\pm$ 61 | 1.46
| 5 | 2 | 7 | LIFO | 1,331 | 11 | 101 | 5 |  -1,551 $\pm$ 62 | 2 |  -1,590 $\pm$ 64 | 2.49
| 6 | 2 | 7 | FIFO | 1,331 | 11 | 101 | 5 |  -1,461 $\pm$ 58 | 2 |  -1,495 $\pm$ 60 | 2.31
| 7 | 2 | 10 | LIFO | 1,331 | 11 | 101 | 6 |  -1,569 $\pm$ 61 | 2 |  -1,606 $\pm$ 64 | 2.35
| 8 | 2 | 10 | FIFO | 1,331 | 11 | 101 | 5 |  -1,469 $\pm$ 59 | 2 |  -1,504 $\pm$ 60 | 2.41
3 | 1 | 1 | 7 | LIFO | 1,331 | 11 | 101 | 5 |  -1,490 $\pm$ 58 | 2 |  -1,500 $\pm$ 59 | 0.71
| 2 | 1 | 7 | FIFO | 1,331 | 11 | 101 | 5 |  -1,424 $\pm$ 56 | 2 |  -1,435 $\pm$ 52 | 0.74
| 3 | 1 | 10 | LIFO | 1,331 | 11 | 101 | 5 |  -1,498 $\pm$ 61 | 2 |  -1,512 $\pm$ 58 | 0.90
| 4 | 1 | 10 | FIFO | 1,331 | 11 | 101 | 5 |  -1,425 $\pm$ 55 | 2 |  -1,436 $\pm$ 52 | 0.82
| 5 | 2 | 7 | LIFO | 14,641 | 11 | 101 | 13 |  -1,513 $\pm$ 61 | 2 |  -1,533 $\pm$ 61 | 1.32
| 6 | 2 | 7 | FIFO | 14,641 | 11 | 101 | 13 |  -1,435 $\pm$ 56 | 2 |  -1,456 $\pm$ 58 | 1.42
| 7 | 2 | 10 | LIFO | 14,641 | 11 | 101 | 13 |  -1,526 $\pm$ 60 | 2 |  -1,544 $\pm$ 61 | 1.16
| 8 | 2 | 10 | FIFO | 14,641 | 11 | 101 | 13 |  -1,437 $\pm$ 56 | 2 |  -1,457 $\pm$ 58 | 1.42
4 | 1 | 1 | 7 | LIFO | 14,641 | 11 | 101 | 14 |  -1,459 $\pm$ 56 | 2 |  -1,476 $\pm$ 54 | 1.15
| 2 | 1 | 7 | FIFO | 14,641 | 11 | 101 | 14 |  -1,422 $\pm$ 56 | 2 |  -1,430 $\pm$ 52 | 0.54
| 3 | 1 | 10 | LIFO | 14,641 | 11 | 101 | 14 |  -1,465 $\pm$ 56 | 2 |  -1,481 $\pm$ 60 | 1.08
| 4 | 1 | 10 | FIFO | 14,641 | 11 | 101 | 14 |  -1,422 $\pm$ 56 | 2 |  -1,430 $\pm$ 52 | 0.54
| 5 | 2 | 7 | LIFO | 161,051 | 11 | 101 | 111 |  -1,480 $\pm$ 59 | 2 |  -1,496 $\pm$ 59 | 1.07
| 6 | 2 | 7 | FIFO | 161,051 | 11 | 101 | 110 |  -1,432 $\pm$ 55 | 2 |  -1,453 $\pm$ 58 | 1.44
| 7 | 2 | 10 | LIFO | 161,051 | 11 | 101 | 110 |  -1,489 $\pm$ 59 | 2 |  -1,505 $\pm$ 58 | 1.07
| 8 | 2 | 10 | FIFO | 161,051 | 11 | 101 | 109 |  -1,432 $\pm$ 55 | 2 |  -1,453 $\pm$ 58 | 1.44
5 | 1 | 1 | 7 | LIFO | 161,051 | 11 | 101 | 114 |  -1,443 $\pm$ 55 | 2 |  -1,454 $\pm$ 55 | 0.73
| 2 | 1 | 7 | FIFO | 161,051 | 11 | 101 | 113 |  -1,422 $\pm$ 56 | 2 |  -1,430 $\pm$ 52 | 0.54
| 3 | 1 | 10 | LIFO | 161,051 | 11 | 101 | 114 |  -1,446 $\pm$ 56 | 2 |  -1,460 $\pm$ 55 | 0.94
| 4 | 1 | 10 | FIFO | 161,051 | 11 | 101 | 114 |  -1,422 $\pm$ 56 | 2 |  -1,430 $\pm$ 52 | 0.54
| 5 | 2 | 7 | LIFO | 1,771,561 | 11 | 101 | 1,191 |  -1,463 $\pm$ 58 | 2 |  -1,480 $\pm$ 60 | 1.22
| 6 | 2 | 7 | FIFO | 1,771,561 | 11 | 101 | 1,185 |  -1,432 $\pm$ 55 | 2 |  -1,453 $\pm$ 58 | 1.44
| 7 | 2 | 10 | LIFO | 1,771,561 | 11 | 101 | 1,188 |  -1,467 $\pm$ 58 | 2 |  -1,484 $\pm$ 59 | 1.15
| 8 | 2 | 10 | FIFO | 1,771,561 | 11 | 101 | 1,190 |  -1,432 $\pm$ 55 | 2 |  -1,453 $\pm$ 58 | 1.44
Table 3: Our results on Scenario A for all of the experimental settings from
[15]. The longest wall times, for value iteration when $m=5$ and $L=2$, are
approximately 20 minutes. Value iteration was considered intractable for
experiments where $m>2$ in the original study.

##  5 Scenario B: substitution between two perishable products

###  5.1 Problem description

[27] applied value iteration and simulation optimization to fit replenishment
policies for two perishable inventory problems: a single-product scenario that
is similar to Scenario A and a scenario with two products and the potential
for substitution which we consider here as Scenario B. In Scenario B we manage
two perishable products, product A and product B, with the same fixed, known
useful life $m$. Some customers who want product B are willing to accept
product A instead if product B is out of stock. The lead time $L=1$ and
therefore there is no in transit component to the state.

At the start of each day $t$, the agent observes state $S_{t}$, the current
inventory of each product in stock split by remaining useful life, and places
a replenishment order. The action consists of two elements, one order for each
product: $A_{t}=[A^{a}_{t},A^{b}_{t}]$ where
$A^{a}_{t}\in\\{0,1,...,A^{a}_{\max}\\}$ and
$A^{b}_{t}\in\\{0,1,...,A^{b}_{\max}\\}$. Demand for day $t$ is sampled from
independent Poisson distributions for each product, parameterised respectively
by mean demand $\mu^{a}$ and $\mu^{b}$, and is initially filled for each
product independently using a FIFO issuing policy. Some customers with unmet
demand for product B may be willing to accept product A instead. The
substitution demand is sampled from a binomial distribution, with a
probability of accepting substitution $\rho$ and a number of trials equal to
the unmet demand for product B. After demand for product A has been filled as
far as possible, demand for product B willing to accept product A is filled by
any remaining units of product A using a FIFO issuing policy. At the end of
the day, the state is updated to reflect the ageing of stock, and the reward,
$R_{t+1}$ is calculated. The reward function comprises revenue per unit sold
($C_{r}^{a},C_{r}^{b}$) and variable order cost ($C_{v}^{a},C_{v}^{b}$) for
each product. The order placed on day $t$ is received immediately prior to the
start of day $t+1$ and is included in the stock element of state $S_{t+1}$.

The daily demand and willingness to accept substitution are both stochastic.
We capture the effect of both by considering the stochastic element in the
transition to be the number of units issued for each product type: $H^{a}$ and
$H^{b}$. The state transition and the reward are deterministic given a state-
action pair and the number of units issued of product A and of product B. The
set of possible realisations of the stochastic elements is:

$\displaystyle\mathbb{\Omega}=\\{(h^{a},h^{b})\\}\quad\quad$ $\displaystyle
h^{a}\in\\{0,1,...,H^{a}_{\max}=mA^{a}_{\max}\\}$ (8) $\displaystyle
h^{b}\in\\{0,1,...,H^{b}_{\max}=mA^{b}_{\max}\\}$

The initial value function, $V_{0}(s)$, is set to the expected sales revenue
for state $s$ with $I^{a}$ units of product A and $I^{b}$ units of product B
in stock. This is an infinite horizon problem with no discount factor, and
therefore we used the convergence test specified in [27], which stops value
iteration when the value of each state is changing by approximately the same
amount on each iteration. If the value of every state is changing by the same
amount there will be no further changes to the best action for each state,
indicating a stable estimate of the optimal policy.

[27] considered products with a maximum useful life $m$ of two and three
periods and evaluated two experimental settings for each value of $m$. For
experiments 1 and 2 with $m=3$, where the maximum order quantities were set
based on the newsvendor model, they reported that it was not possible to
complete value iteration within one week. They therefore repeated the two
experiments for $m=3$ with lower values of $A^{a}_{\max}$ and $A^{b}_{\max}$.
With this adjustment, one of the cases could be completed within 80 hours,
while the other could still not be solved within a week. We demonstrate that
using JAX and a consumer-grade GPU it is feasible to obtain the optimal policy
for all of these settings and report the wall time required to run value
iteration for each experiment. Additionally, to investigate how our method can
benefit from more powerful GPUs, and how it scales to multiple GPUs, we report
the wall times for running the largest problem on one, two and four Nvidia
A100 40GB GPUs.

Separately, we consider the experimental settings used by [45] to evaluate
their GPU-accelerated method, in which value iteration was always run for 100
iterations instead of to convergence. The different experiments evaluate mean
daily demands between five and seven with maximum order quantities based on
the newsvendor model.

In each case, we compare the policy from value iteration with the modified
base-stock policy used by [27], based on the work of [25], which has an order-
up-to level parameter for each product: $\texttt{S}^{a}$ and $\texttt{S}^{b}$.
The order quantity for each product is determined considering only the on hand
inventory of that product and includes an adjustment for expected waste. The
order quantity on day $t$, given total stock on hand $I^{a}_{t}$ and
$I^{a}_{t}$ and stock that expires at the end of the current period
$X^{a}_{1,t}$ and $X^{b}_{1,t}$, is:

$A_{t}=\left[A^{a}_{t},A^{b}_{t}\right]=\biggl{[}\Bigl{[}\texttt{S}^{a}-I^{a}_{t}+\bigl{[}X^{a}_{1,t}-\mu^{a}\bigr{]}^{+}\Bigr{]}^{+},\Bigl{[}\texttt{S}^{b}-I^{b}_{t}+\bigl{[}X^{b}_{1,t}-\mu^{b}\bigr{]}^{+}\Bigr{]}^{+}\biggr{]}$
(9)

There are two parameters, and we used Optuna’s NSGAII sampler to search the
parameter space
$\texttt{S}^{a}\in\\{0,1,...,\texttt{S}^{a}_{\max}=2A^{a}_{\max}\\}$ and
$\texttt{S}^{b}\in\\{0,1,...,\texttt{S}^{b}_{\max}=2A^{b}_{\max}\\}$. We
considered values of the order-up-to level up to twice the maximum order
quantity used for value iteration because [27] reported best values of S that
were higher than the values of $A_{\max}$ specified for value iteration for
some of their experiments. We compare the modified base-stock policy that
achieved the highest mean return, characterised by the pair of parameters
$\left(\texttt{S}^{a},\texttt{S}^{b}\right)_{\text{best}}$, to the value
iteration policy.

See Appendix B.1 for additional information about Scenario B.

###  5.2 Results

We present results for the experimental settings for the two product scenario
from [27] in Table 4. The wall times in Table 4 show that, using our method,
value iteration can be used to find the optimal policy for all four settings
of the two product scenario with $m=3$ in under 3.2 hours. [27] reported that,
for $m=3$, value iteration did not converge within a week for experiments 1, 2
and 3 using a MATLAB implementation and experiment 4 converged in 80 hours.
Our implementation of experiment 4, running on a consumer-grade GPU, converges
in just over two minutes: more than 2000${\times}$ faster.

We present results for the four experimental settings, P1 to P4, from [45] in
Table 5. The wall times for our approach are at least six times faster than
those reported by [45] for all four settings. We cannot conclude on the
relative performance of our method and the GPU-accelerated method from [45]
without running both implementations on the same hardware and accounting for
the difference between up-front and JIT compilation. However, the results
suggest that our method is at least competitive with a custom CUDA
implementation of value iteration for the two product case while requiring
less specialist knowledge of GPU programming.

Simulation optimization scales well to larger problems, with wall times less
than one minute for all of the experimental settings. The optimality gap is
never greater than 1%, and reduces as both mean demand and the maximum useful
life increase. This suggests that there is a limited advantage to making
ordering decisions based on the stock of both products, compared to making
independent decisions for each product using a simple heuristic policy, under
the reward function and substitution process proposed by [27].

Figure 1 illustrates the clear benefits of both more powerful GPUs, and of
using multiple GPUs. Using a single Nvidia A100 40GB GPU, experiment 1 when
$m=3$ can be run in 4,838s: 2.4$\times$ faster than the Nvidia RTX 3060 in our
local machine. The A100 40GB has more GPU RAM and more CUDA cores than the RTX
3060 [57], 40GB vs 16GB and 6,912 vs 3,584 respectively, which means that it
can update the value function for a larger number of states simultaneously.
Using two A100 40GB GPUs is 1.8$\times$ faster than one, and using four A100
40GB GPUs is 2.8$\times$ faster than one, demonstrating how the wall time can
further reduced and how larger problems can be solved with additional
computational resources using exactly the same code.

See Appendix B.2 for the best parameters for the heuristic policy and KPIs for
each experiment.

|  |  | Value | Simulation
---|---|---|---|---
|  |  | iteration | optimization
$m$ | Exp | $\mu^{a}$ | $\mu^{b}$ | $A^{a}_{\max}$ | $A^{b}_{\max}$ | $|\mathbb{S}|$ | $|\mathbb{A}|$ | $|\mathbb{\Omega}|$ | WT (s) | Return | WT (s) | Return | Optimality
|  |  |  |  |  |  |  |  |  |  |  |  | gap (%)
2 | 1 | 5 | 5 | 10 | 10 | 14,641 | 121 | 441 | 5 |  1,644 $\pm$ 33 | 24 |  1,632 $\pm$ 34 | 0.70
| 2 | 7 | 3 | 14 | 6 | 11,025 | 105 | 377 | 4 |  1,650 $\pm$ 33 | 23 |  1,639 $\pm$ 34 | 0.67
3 | 1 | 5 | 5 | 15 | 15 | 16,777,216 | 256 | 2,116 | 11,496 |  1,761 $\pm$ 32 | 33 |  1,758 $\pm$ 32 | 0.16
| 2 | 7 | 3 | 21 | 9 | 10,648,000 | 220 | 1,792 | 4,013 |  1,762 $\pm$ 32 | 44 |  1,759 $\pm$ 32 | 0.18
| 3 | 5 | 5 | 13 | 13 | 7,529,536 | 196 | 1,600 | 3,058 |  1,761 $\pm$ 32 | 32 |  1,758 $\pm$ 32 | 0.16
| 4 | 7 | 3 | 20 | 4 | 1,157,625 | 105 | 793 | 134 |  1,762 $\pm$ 32 | 43 |  1,759 $\pm$ 32 | 0.17
Table 4: Our results on Scenario B for all of the experimental settings from [27]. The longest wall time, for value iteration on experiment 1 when $m=3$, is approximately $3.2$ hours. Value iteration could not be completed within a week for experiments 1, 2 and 3, and required 80 hours for experiment 4, when $m=3$ in [27]. |  |  | Value | Simulation
---|---|---|---|---
|  |  | iteration | optimization
$m$ | Exp | $\mu^{a}$ | $\mu^{b}$ | $A^{a}_{\max}$ | $A^{b}_{\max}$ | $|\mathbb{S}|$ | $|\mathbb{A}|$ | $|\mathbb{\Omega}|$ | WT (s) | Return | WT (s) | Return | Optimality
|  |  |  |  |  |  |  |  |  |  |  |  | gap (%)
2 | P1 | 5 | 5 | 10 | 10 | 14,641 | 121 | 441 | 11 |  1,644 $\pm$ 33 | 25 |  1,632 $\pm$ 34 | 0.70
| P2 | 5 | 6 | 10 | 12 | 20,449 | 143 | 525 | 18 |  1,826 $\pm$ 35 | 31 |  1,816 $\pm$ 34 | 0.58
| P3 | 6 | 6 | 12 | 12 | 28,561 | 169 | 625 | 27 |  2,011 $\pm$ 36 | 31 |  2,000 $\pm$ 37 | 0.55
| P4 | 7 | 7 | 13 | 13 | 38,416 | 196 | 729 | 56 |  2,379 $\pm$ 39 | 29 |  2,368 $\pm$ 40 | 0.46
Table 5: Our results on Scenario B for all of the experimental settings used
by [45] to test their GPU-accelerated approach. Value iteration was run for
100 iterations for each experiment instead of to convergence. Aside from this,
experiment P1 is the same as experiment 1 with $m=2$ in Table 4. The longest
wall time, for value iteration on experiment P4, is approximately one minute.
Value iteration was tractable for all of these settings in the original study
but wall times using our method are at least six times faster than those
reported by [45]. This improvement may be at least partially attributable to
hardware differences. Figure 1: Wall times required to run value iteration for
experiment 1 with $m=3$ for Scenario B using different GPUs. The Nvidia
GeForce RTX 3060 is a consumer-grade GPU. The Nvidia A100 40GB is a data-
centre grade GPU. JAX enables our method to be run on multiple identical GPUs
without any code changes.

##  6 Scenario C: periodic demand and uncertain useful life on
arrival

###  6.1 Problem description

[40] described a perishable inventory problem that models the management of
platelets in a hospital blood bank. There are two problem features not
included in Scenarios A or B. Firstly, the demand is periodic, with an
independent demand distribution for each day of the week. Secondly, the
remaining useful life of products on arrival is uncertain, and this
uncertainty may be exogenous or endogenous. The lead time $L$ is assumed to be
zero and therefore there is no in transit component to the state.

At the start of each day $t$ the agent observes state $S_{t}$, which specifies
the day of the week and the current inventory in stock split by remaining
useful life, and places a replenishment order
$A_{t}\in\\{0,1,...,A_{\max}\\}$. This order is assumed to arrive instantly.
The remaining useful life of the units on arrival is sampled from a
multinomial distribution, the parameters of which may depend on the order
quantity $A_{t}$. Demand for day $t$, $D_{t}$ is sampled from a truncated
negative binomial distribution and filled from available stock using an
oldest-unit first-out (OUFO) policy. At the end of the day, the state is
updated to reflect the ageing of stock and the reward, $R_{t+1}$ is
calculated. The reward function comprises four components: a holding cost per
unit in stock at the end of the period ($C_{h}$), a shortage cost per unit of
unmet demand ($C_{s}$), a wastage cost per unit that perishes at the end of
the period ($C_{w}$) and a fixed ordering cost ($C_{f}$). Unlike Scenario A
which also includes a holding cost, the holding cost is charged on units that
expire at the end of the period.

To reduce the number of possible states, we consider a limited case of this
problem in which there is a maximum capacity of $A_{\max}$ for stock of each
age. If, when an order is received, the sum of units in stock and units
received with $k$ days of remaining useful life is greater than $A_{\max}$
then we assume the excess units are not accepted at delivery. The stock level
with $k$ days of remaining useful life is therefore at most $A_{\max}$ when
demand is sampled. This constraint is chosen to be consistent with the
calculation of the total number of states in [40], but there are alternative
ways to apply the constraint (e.g. by discarding excess units at the end of
each day along with wastage) and these may have different optimal policies.

The stochastic elements in the transition are the daily demand, $D$, and the
age profile of the units received to fill the order placed at the start of the
day: $\underline{\text{Y}}=\left[Y_{m},Y_{m-1},...,Y_{1}\right]$. The state
transition and the reward are deterministic given a state-action pair, the
daily demand, and the age profile of the units received. The set of possible
realisations of the stochastic elements is:

$\displaystyle\mathbb{\Omega}=\\{(d,\underline{\text{y}})\\}\quad$
$\displaystyle d\in\\{0,1,...,D_{\max}\\}$ (10) $\displaystyle
y_{i}\in\\{0,1,...,A_{\max}\\},\quad\forall i\in\\{1,2,...,m\\}$
$\displaystyle\sum_{i=1}^{m}y_{i}\leq A_{\max}$

The initial value function $V_{0}(s)$ was initialised at zero for every state.
[40] did not specify a particular convergence test for his value iteration
experiments. The problem is periodic, with a discount factor, and therefore we
use a convergence test based on those described in [54] which stops value
iteration when the undiscounted change in the value function over a period (in
this case, seven days) is approximately the same for every state. As in
Scenario B, when the change in value is the approximately the same for every
state there will be no further changes to the best action for every state, and
hence, the estimated optimal policy is stable.

[40] considered products with a maximum useful life of three, five or eight
periods and stated that, due to the large state space, value iteration was
intractable for this problem when $m\geq 5$. We were able to run value
iteration when $m=5$, but not when $m=8$. For each value of $m$, [40]
investigated five different settings for the distribution of useful life on
arrival (one where the uncertainty was exogenous, and four where the
uncertainty was endogenous). For each of these five settings, he evaluated six
combinations of $C_{f}$ and $C_{w}$. Our objective was to demonstrate the
feasibility of our approach and therefore, given the large number of
experiments and long wall times when $m=5$, we ran two experiments for each
value of $m$: one where the uncertainty in useful life on arrival was
exogenous and one where it was endogenous. For $m=5$, we selected the settings
from [40] that are based on real, observed data from a network of hospitals in
Ontario, Canada instead of the additional settings created for sensitivity
analysis. We report the wall time required to run value iteration for each
experiment.

We compare the policy from value iteration with an $(\texttt{s},\texttt{S})$
policy. [40] did not fit heuristic policies, but suggested (s, S) as an
example of a suitable heuristic policy for future work: the addition of a
fixed ordering cost to the reward function means that it may be beneficial to
include the reorder point parameter s to avoid uneconomically small orders. We
fit one pair of s and S for each day of the week, a total of 14 parameters.
The order quantity on day $t$, given that the day of the week is $\tau$ and
the total current stock on hand is $I_{t}$ is:

$A_{t}=\begin{cases}\left[\texttt{S}^{\tau}-I_{t}\right]^{+}&\text{if
$I_{t}\leq\texttt{s}^{\tau}$}\\\ 0&\text{if
$I_{t}>\texttt{s}^{\tau}$}\end{cases}$ (11)

where $(\texttt{s}^{\tau},\texttt{S}^{\tau})$ is the pair of parameters for
day of the week $\tau$.

We used Optuna’s NSGAII sampler to search for combinations of
$\texttt{s}^{\tau}\in\\{0,1,...,\texttt{s}_{\max}=A_{\max}\\}$ and
$\texttt{S}^{\tau}\in\\{0,1,...,\texttt{S}_{\max}=A_{\max}\\}\;\forall\tau\in\\{0,1,..,6\\}$.
This heuristic policy has a hard constraint that
$\texttt{s}^{\tau}<\texttt{S}^{\tau}\;\forall\tau\in\\{0,1,..,6\\}$. Optuna
does not support using hard constraints to restrict the search space, so we
enforced the constraint by only allowing a non-zero order to be placed if the
constraint was met. We compare the heuristic policy that achieved the highest
mean return, characterised by parameters
$\left(\left(\texttt{s}^{0},\texttt{S}^{0}\right),...,\left(\texttt{s}^{6},\texttt{S}^{6}\right)\right)_{\text{best}}$,
to the value iteration policy.

See Appendix C.1 for additional information about Scenario C.

###  6.2 Results

In Table 6 we present the results for the experimental settings from [40] that
we have selected, two for each value of $m$. Using our method, it is possible
to find the optimal policy using value iteration for $m=3$ and $m=5$ while
accounting for uncertainty in useful life on arrival. The experiments where
$m=5$ represent a real world problem: [40] fit the parameters for the demand
distribution and distribution of useful life on arrival to observed data from
a network of hospitals in Ontario, Canada. This is an important application of
our value iteration method, demonstrating that it can be used to find optimal
policies for problems of a realistic size. The alternative experimental
settings evaluated by [40] but not repeated here have the same numbers of
states, actions and possible random outcomes and therefore we would expect the
wall times to be of a similar order as corresponding experiments reported in
Table 6.

We were unable to complete value iteration when $m=8$. This problem has over
12.6 billion possible states, even with the restriction that we placed on the
maximum stock holding of each age, and over 65 million possible random
outcomes. It is not feasible to store the state array in the memory of our
local machine, let alone run value iteration. However, we were able to fit a
heuristic policy using simulation optimization in less than 20 minutes.

The simulation optimization experiments for this scenario take longer than
those of the other scenarios, between five and 20 minutes. This is due to the
large number of possible combinations of parameters, because our heuristic
policy require seven pairs of parameters $\left(\texttt{s},\texttt{S}\right)$,
one for each weekday. Optuna does not support restricting the search space
based on the constraint that
$\texttt{s}^{\tau}<\texttt{S}^{\tau}\;\forall\tau\in\\{0,1,..,6\\}$ and
therefore the size of the search space for each experiment is
$(A_{\max}+1)^{14}=3.2\times 10^{18}$, compared to only 11 possible parameters
for the base-stock policy used for Scenario A and fewer than 1,000 possible
combinations of parameters for even the largest scenarios from Scenario B. The
heuristic policies perform well, with a maximum optimality gap of 1.22%. See
Appendix C.2 for the best parameters for the heuristic policy and KPIs for
each experiment.

|  |  | Value | Simulation
---|---|---|---|---
|  |  | iteration | optimization
$m$ | Exp | Uncertainty | $|\mathbb{S}|$ | $|\mathbb{A}|$ | $|\mathbb{\Omega}|$ | WT (s) | Mean return | WT (s) | Mean return | Optimality
|  | in useful life |  |  |  |  |  |  |  | gap (%)
3 | 1 | Exogenous | 3,087 | 21 | 37,191 | 15 |  -410 $\pm$ 62 | 507 |  -411 $\pm$ 63 | 0.26
| 2 | Endogenous | 3,087 | 21 | 37,191 | 17 |  -349 $\pm$ 53 | 305 |  -352 $\pm$ 55 | 1.04
5 | 1 | Exogenous | 1,361,367 | 21 | 1,115,730 | 178,078 |  -312 $\pm$ 46 | 514 |  -313 $\pm$ 50 | 0.34
| 2 | Endogenous | 1,361,367 | 21 | 1,115,730 | 178,023 |  -312 $\pm$ 47 | 393 |  -315 $\pm$ 46 | 1.22
8 | 1 | Exogenous | 12,607,619,787 | 21 | 65,270,205 | — | — | 618 |  -293 $\pm$ 42 | —
| 2 | Endogenous | 12,607,619,787 | 21 | 65,270,205 | — | — | 972 |  -297 $\pm$ 43 | —
Table 6: Our results on Scenario C for a subset of the experimental settings
from [40]: two examples for each value of $m$. The longest wall time, for
value iteration on experiment 1 when $m=5$, is approximately 49.5 hours. Value
iteration was considered intractable for the experiments where $m\geq 5$ in
the original study. We were able to use value iteration when $m=5$, but not
when $m=8$.

In Figure 2 we draw together the results from Scenario C with those from the
preceding scenarios, and plot the optimality gap between the heuristic policy
that achieved the highest mean return and the value iteration policy against
the wall time required for value iteration.

Figure 2: The optimality gap between the heuristic policy that achieved the
highest mean return and the value iteration policy plotted against the wall
time required to run value iteration for the experiments from Scenarios A, B
and C.

##  7 Discussion

We have found JAX to provide an effective way to expand the scale of
perishable inventory problems for which value iteration is tractable, using
only consumer-grade hardware. Expanding the range of problems to which value
iteration can be applied is not just a matter of considering settings with
greater demand, more products, or products with a longer useful life. It also
enables us to incorporate complexity that might otherwise be neglected due to
its effect on the computational tractability such as substitution and
endogenous uncertainty in the useful life of products on arrival. An “optimal”
policy fit using value iteration is optimal for the situation as modelled, but
may not perform well in practice if the model neglects challenging aspects of
the real problem. For example, [40] reported large optimality gaps, with an
average of 51%, when policies obtained under the assumption that all stock
arrived fresh were applied to a scenario with endogenous uncertainty in useful
life. One avenue for future work is to consider scenarios that combine the
more challenging elements: a multi-period lead time, substitution between
multiple products, and uncertainty in the useful life on arrival which may all
be relevant to managing blood product inventory in reality.

Our simulation optimization approach scales well to larger problems and large
policy parameter spaces, and performs well relative to the optimal policies.
One benefit of increasing the size of problems for which value iteration is
tractable is being able to better understand how the relative of performance
of heuristic and approximate policies scales with properties that influence
the problem size. This may help with the development of new heuristics, and
determining the utility of reinforcement learning and other approximate
methods. The optimality gap in our experiments was never larger than 2.5%, and
in the experiments from Scenario A and Scenario B the optimality gap decreased
as the demand and/or maximum useful life of the product increased. This is
encouraging because it suggests that in some circumstances where the problem
size remains too large for value iteration there may actually be little to
gain by using the optimal policy over one of these heuristic policies.

In our simulation optimization experiments we only used GPUs to run the
simulated rollouts. The heuristic search methods for proposing the next sets
of candidate parameters are CPU-based. We did not find Optuna’s NSGAII sampler
to be a bottleneck, but during preliminary experiments we found that some
alternative methods took longer to propose the next set of candidate
parameters than was required to evaluate them on simulated rollouts. In future
work, the optimization process suggesting parameters could also be run on GPU
similar to the work of [38] and recent work using evolutionary strategies on
GPUs to search for neural network parameters [36]. The gymnax-based simulators
would also be well suited to ranking and selection methods because it would be
straightforward to run a small number of rollouts for a large number of
possible parameters in parallel and then, at a second stage, run a large
number of rollouts for the most competitive parameters in parallel to obtain
more accurate estimates of their performance.

One of the main contributions of this work is to demonstrate an accessible way
of using GPUs to accelerate value iteration and simulation optimization and
therefore solve larger problems that are closer to those faced in reality. On
the software side, we implemented our approach using the relatively high-level
JAX API and relied on the XLA compiler to efficiently utilise GPU hardware. On
the hardware side, we primarily report results on a consumer-grade GPU, and
make available a Google Colab notebook so that our experiments can be
reproduced at no cost using cloud-based computational resources. However, a
significant strength of JAX is support for easily distributing a workload over
multiple identical GPU devices using the pmap function transformation and we
discuss in Section 5 how additional devices can be used to further reduce the
wall time and potentially make even larger problems tractable. Modern cloud
computing platforms provide on-demand access to data-centre grade GPUs,
including the A100 40GB GPU we used to run the scaling experiments in Section
5. At the time of writing in February 2023, a single A100 40GB GPU is
available on-demand for $3.67 per hour and four A100 40GB GPUs are available
for $14.68 per hour through Google Cloud Platform [33]. This may provide a
cost-effective way for research teams without access to local high-performance
computing resources to investigate problems that are too large for freely
available or consumer-grade GPU hardware.

For some cases it may be possible to further reduce the wall time for value
iteration on the same hardware by using single-precision numbers instead of
double-precision numbers. We experienced instabilities in convergence during
preliminary experiments with a large number of iterations which we resolved by
changing from JAX’s default single-precision to double-precision. The
additional precision comes with a performance cost. Experiment 1 when $m=5$
for Scenario C has the longest wall time at double-precision: 49.5 hours. This
can be reduced to 23.5 hours when using single-precision, over twice as fast
at single-precision, and the final policy is the same using both approaches. A
similar reduction in wall time can be obtained for experiment 2 with $m=5$ for
Scenario C, but in this case the final policies are not identical. The best
order quantities for just six out of 1,115,730 states differ from the policy
found at double-precision, and each by one unit. This suggests that there may
be a significant potential benefit to using single-precision numbers for GPU-
accelerated value iteration if the user is willing to tolerate potentially
larger approximation errors.

Future hardware development will make value iteration feasible for even larger
problems. In addition to future generations of GPUs, one promising direction
is field programmable gate arrays (FPGAs): integrated circuits that can be
reprogrammed to customise the hardware to implement a specific algorithm,
including value iteration [47]. Customising the hardware currently requires
specialist knowledge but, just as machine learning frameworks and higher-level
tools have made GPU programming more accessible, [47] suggests that FPGA
compilers able to translate high level code into customised circuit designs
may facilitate wider adoption.

We have focused on perishable inventory management in this study, but our
computational approach has much wider applicability. For each scenario we
created a custom subclass of our base value iteration runner class and a
custom subclass of the gymnax reinforcement learning environment as our
simulator, each with methods to implement the scenario-specific logic. The
same approach could be followed for any problem that can be modelled as an
MDP. More broadly, we believe that JAX (and other software libraries
originally developed to support deep learning including PyTorch and
Tensorflow) offers an efficient way for researchers to run large workloads in
parallel on relatively affordable GPU hardware which may support research on a
range of operational research problems.

##  8 Conclusion

JAX and similar software libraries provide a way for researchers without
extensive experience of GPU programming to take advantage of the parallel
processing capabilities of modern GPUs. In this study we have shown how a JAX-
based approach can expand the range of perishable inventory management
problems for which value iteration is tractable, using only consumer-grade
hardware. We also created GPU-accelerated simulators for each scenario, in the
form of JAX-based reinforcement learning environments, and demonstrated how
these can be used to quickly fit the parameters of heuristic policies by
simultaneously evaluating many sets of policy parameters on thousands of
simulated rollouts in parallel. By reducing the wall time required to run
value iteration and simulation optimization, these methods can support
research into larger problems, both in terms of scale and the incorporation of
aspects of reality that increase the computational complexity. The ability to
find optimal policies using value iteration may provide a valuable benchmark
for the evaluation of new heuristic and approximate methods, helping efforts
to make the best use of scarce resources and reduce wastage of perishable
inventory. This work is focused on perishable inventory management but we
believe that our methods, and the underlying principle of using software
developed by the machine learning community to parallelize workloads on GPU,
may be applicable to many problems in operational research and have made our
code publicly available to support future work.

## Author roles

Joseph Farrington: Conceptualization, Methodology, Software, Investigation,
Validation, Writing - Original Draft Kezhi Li: Supervision, Writing - Review &
Editing Wai Keong Wong: Supervision Martin Utley: Investigation, Supervision,
Writing - Review & Editing

## Acknowledgements

The authors are grateful to Professor Eligius Hendrix and Dr Mahdi Mirjalili
for providing additional material that enabled us to test our implementation
of the scenarios from their work. Any errors or differences are our
responsibility alone. The authors would also like to thank Dr Thomas Monks for
sharing his expertise on simulation optimization during the preliminary stages
of this work.

## Funding Statement

JF is supported by UKRI training grant EP/S021612/1, the CDT in AI-enabled
Healthcare Systems. This study was supported by the Clinical and Research
Informatics Unit at the National Institute for Health and Care Research
University College London Hospitals Biomedical Research Centre.

The authors acknowledge the use of the UCL Myriad High Performance Computing
Facility (Myriad@UCL), and associated support services, in the completion of
this work.

The sponsors of the research did not have a role in the study design, in the
collection, analysis and interpretation of the data, in the writing of this
report, or in the decision to submit this article for publication.

For the purpose of open access, the authors have applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

## Competing interests

The authors have no competing interests to declare.

## References

  * [1] Ammar Aamer, LuhPutu Eka Yani and IMade Alan Priyatna “Data analytics in the supply chain management: Review of machine learning applications in demand forecasting” In _Operations and Supply Chain Management: An International Journal_ 14.1, 2020, pp. 1–13 DOI: http://doi.org/10.31387/oscm0440281
  * [2] Martín Abadi et al. “TensorFlow: a system for large-scale machine learning” In _Proceedings of the 12th USENIX Conference on Operating Systems Design and Implementation_ , 2016, pp. 265–283 DOI: 10.48550/arXiv.1605.08695
  * [3] Ehsan Ahmadi et al. “Intelligent inventory management approaches for perishable pharmaceutical products in a healthcare supply chain” In _Computers & Operations Research_ 147, 2022, pp. 105968 DOI: 10.1016/j.cor.2022.105968
  * [4] Takuya Akiba et al. “Optuna: A next-generation hyperparameter optimization framework” In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, KDD ’19 New York, NY, USA: Association for Computing Machinery, 2019, pp. 2623–2631 DOI: 10.1145/3292500.3330701
  * [5] Eric M. Aldrich, Jesús Fernández-Villaverde, A. Ronald Gallant and Juan F. Rubio-Ramírez “Tapping the supercomputer under your desk: Solving dynamic equilibrium models with graphics processors” In _Journal of Economic Dynamics and Control_ 35.3, 2011, pp. 386–393 DOI: 10.1016/j.jedc.2010.10.001
  * [6] Satyajith Amaran, Nikolaos V. Sahinidis, Bikram Sharda and Scott J. Bury “Simulation optimization: a review of algorithms and applications” In _Annals of Operations Research_ 240.1, 2016, pp. 351–380 DOI: 10.1007/s10479-015-2019-x
  * [7] R. Bellman “Dynamic Programming” Princeton, NJ, USA: Princeton University Press, Princeton, NJ, USA, 1957
  * [8] John T. Blake et al. “Optimizing the platelet supply chain in Nova Scotia” In _Proceedings of the 29th meeting of the European Working Group on Operational Research Applied to Health Services (ORAHS 2003)_ , 2003, pp. 47–66 URL: http://orahs.di.unito.it/docs/2003-ORAHS-proceedings.pdf#page=47
  * [9] Clément Bonnet et al. “Jumanji: Industry-driven hardware-accelerated RL environments”, 2022 URL: https://github.com/instadeepai/jumanji
  * [10] James Bradbury et al. “JAX: composable transformations of Python+NumPy programs”, 2018 URL: http://github.com/google/jax
  * [11] Vaibhav Chaudhary, Rakhee Kulshrestha and Srikanta Routroy “State-of-the-art literature review on inventory models for perishable products” In _Journal of Advances in Management Research_ 15.3, 2018, pp. 306–346 DOI: http://dx.doi.org/10.1108/JAMR-09-2017-0091
  * [12] Peng Chen and Lu Lu “Markov decision process parallel value iteration algorithm on GPU” In _Proceedings of 2013 International Conference on Information Science and Computer Applications_ Atlantis Press, 2013, pp. 299–304 DOI: 10.2991/isca-13.2013.51
  * [13] Denisa-Andreea Constantinescu, Angeles Navarro, Juan-Antonio Fernández-Madrigal and Rafael Asenjo “Performance evaluation of decision making under uncertainty for low power heterogeneous platforms” In _Journal of Parallel and Distributed Computing_ 137, 2020, pp. 119–133 DOI: 10.1016/j.jpdc.2019.11.009
  * [14] Doraid Dalalah, Omar Bataineh and Khaled A. Alkhaledi “Platelets inventory management: a rolling horizon sim–opt approach for an age-differentiated demand” In _Journal of Simulation_ 13.3, 2019, pp. 209–225 DOI: 10.1080/17477778.2018.1497461
  * [15] Bram J. De Moor, Joren Gijsbrechts and Robert N. Boute “Reward shaping to improve the performance of deep reinforcement learning in perishable inventory management” In _European Journal of Operational Research_ 301.2, 2022, pp. 535–545 DOI: 10.1016/j.ejor.2021.10.045
  * [16] Mary Dillon, Fabricio Oliveira and Babak Abbasi “A two-stage stochastic programming model for inventory management in the blood supply chain” In _International Journal of Production Economics_ 187.May 2016, 2017, pp. 27–41 DOI: 10.1016/j.ijpe.2017.02.006
  * [17] Qinglin Duan and T. Liao “A new age-based replenishment policy for supply chain inventory optimization of highly perishable products” In _International Journal of Production Economics_ 145.2, 2013, pp. 658–671 DOI: 10.1016/j.ijpe.2013.05.020
  * [18] Victor Duarte, Diogo Duarte, Julia Fonseca and Alexis Montecinos “Benchmarking machine-learning software and hardware for quantitative economics” In _Journal of Economic Dynamics and Control_ 111, 2020, pp. 103796 DOI: 10.1016/j.jedc.2019.103796
  * [19] Ebuyer “NVIDIA GeForce RTX 3060 Graphics Card” URL: https://www.ebuyer.com/store/Components/cat/Graphics-Cards-Nvidia/subcat/GeForce-RTX-3060?q=nvidia+3060
  * [20] Andrew WJ Flint et al. “Is platelet expiring out of date? A systematic review” In _Transfusion Medicine Reviews_ 34.1, 2020, pp. 42–50 DOI: 10.1016/j.tmrv.2019.08.006
  * [21] Daniel Freeman et al. “Brax - A differentiable physics engine for large scale rigid body simulation” In _Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks_ 1, 2021 DOI: 10.48550/arXiv.2106.13281
  * [22] Brant E. Fries “Optimal ordering policy for a perishable commodity with fixed lifetime” In _Operations Research_ 23.1, 1975, pp. 46–61 DOI: 10.1287/opre.23.1.46
  * [23] Michael C. Fu et al. “Simulation optimization: A panel on the state of the art in research and practice” In _Proceedings of the Winter Simulation Conference 2014_ , 2014, pp. 3696–3706 DOI: 10.1109/WSC.2014.7020198
  * [24] Serkan Gunpinar and Grisselle Centeno “Stochastic integer programming models for reducing wastages and shortages of blood products at hospitals” In _Computers and Operations Research_ 54, 2015, pp. 129–141 DOI: 10.1016/j.cor.2014.08.017
  * [25] René Haijema and Stefan Minner “Improved ordering of perishables: The value of stock-age information” In _International Journal of Production Economics_ 209, The Proceedings of the 19th International Symposium on Inventories, 2019, pp. 316–324 DOI: 10.1016/j.ijpe.2018.03.008
  * [26] René Haijema, Jan Wal and Nico M. Dijk “Blood platelet production: optimization by dynamic programming and simulation” In _Computers & Operations Research_ 34.3, Logistics of Health Care Management, 2007, pp. 760–779 DOI: 10.1016/j.cor.2005.03.023
  * [27] Eligius MT Hendrix et al. “On computing optimal policies in perishable inventory control using value iteration” In _Computational and Mathematical Methods_ 1.4, 2019, pp. e1027 DOI: https://doi.org/10.1002/cmm4.1027
  * [28] Pieter Hijma et al. “Optimization techniques for GPU programming” In _ACM Computing Surveys_ , 2022 DOI: 10.1145/3570638
  * [29] Tsutomu Inamoto et al. “An implementation of dynamic programming for many-core computers” In _SICE Annual Conference 2011_ , 2011, pp. 961–966 URL: https://ieeexplore.ieee.org/abstract/document/6060648
  * [30] Won Jeon et al. “Chapter Six - Deep learning with GPUs” In _Advances in Computers_ 122, Hardware Accelerator Systems for Artificial Intelligence and Machine Learning Elsevier, Cambridge, MA, USA, 2021, pp. 167–215 DOI: 10.1016/bs.adcom.2020.11.003
  * [31] Arsæll Pór Jóhannsson “GPU-based Markov decision process solver”, 2009 URL: https://en.ru.is/media/skjol-td/MSThesis_ArsaellThorJohannsson.pdf
  * [32] Ahmet Kara and Ibrahim Dogan “Reinforcement learning approaches for specifying ordering policies of perishable inventory systems” In _Expert Systems with Applications_ 91, 2018, pp. 150–158 DOI: 10.1016/j.eswa.2017.08.046
  * [33] Sergey Karayev and Charles Frye “Cloud GPUs” In _Full Stack Deep Learning_ , 2023 URL: https://fullstackdeeplearning.com/cloud-gpus/
  * [34] Robert Kirkby “A toolkit for value function iteration” In _Computational Economics_ 49.1, 2017, pp. 1–15 DOI: 10.1007/s10614-015-9544-1
  * [35] Robert Kirkby “Quantitative macroeconomics: Lessons learned from fourteen replications” In _Computational Economics_ , 2022 DOI: 10.1007/s10614-022-10234-w
  * [36] Robert Tjarko Lange “evosax: JAX-based evolution strategies” arXiv:2212.04180 [cs] arXiv, 2022 DOI: 10.48550/arXiv.2212.04180
  * [37] Robert Tjarko Lange “gymnax: A JAX-based reinforcement learning environment library”, 2022 URL: http://github.com/RobertTLange/gymnax
  * [38] Mai Chan Lau and Rajagopalan Srinivasan “A hybrid CPU-Graphics Processing Unit (GPU) approach for computationally efficient simulation-optimization” In _Computers & Chemical Engineering_ 87, 2016, pp. 49–62 DOI: 10.1016/j.compchemeng.2016.01.001
  * [39] Viktor Makoviychuk et al. “Isaac Gym: High Performance GPU-Based Physics Simulation For Robot Learning” arXiv:2108.10470 [cs] arXiv, 2021 DOI: 10.48550/arXiv.2108.10470
  * [40] Mahdi Mirjalili “Data-driven modelling and control of hospital blood inventory”, 2022 URL: https://tspace.library.utoronto.ca/bitstream/1807/124976/1/Mirjalili_Mahdi_202211_PhD_thesis.pdf
  * [41] Steven Nahmias “A comparison of alternative approximations for ordering perishable inventory” In _INFOR_ 13.2, 1975, pp. 175–184 DOI: 10.1080/03155986.1975.11731604
  * [42] Steven Nahmias “Optimal ordering policies for perishable inventory—II” In _Operations Research_ 23.4, 1975, pp. 735–749 DOI: 10.1287/opre.23.4.735
  * [43] Steven Nahmias “Perishable inventory theory: a review” In _Operations Research_ 30.4, 1982, pp. 680–708 DOI: 10.1287/opre.30.4.680
  * [44] Steven Nahmias “Perishable Inventory Systems”, International Series in Operations Research & Management Science Springer, New York, NY, USA, 2011 DOI: 10.1007/978-1-4419-7999-5
  * [45] G. Ortega, Eligius MT Hendrix and I. García “A CUDA approach to compute perishable inventory control policies using value iteration” In _The Journal of Supercomputing_ 75.3, 2019, pp. 1580–1593 DOI: 10.1007/s11227-018-2692-z
  * [46] Adam Paszke et al. “PyTorch: An Imperative Style, High-Performance Deep Learning Library” In _Advances in Neural Information Processing Systems_ 32 Curran Associates, Inc., 2019 DOI: https://doi.org/10.48550/arXiv.1912.01703
  * [47] Alessandro Peri “A hardware approach to value function iteration” In _Journal of Economic Dynamics and Control_ 114, 2020, pp. 103894 DOI: 10.1016/j.jedc.2020.103894
  * [48] Kalyan Perumalla and Maksudul Alam “Design considerations for GPU-based mixed integer programming on parallel computing platforms” In _50th International Conference on Parallel Processing Workshop_ Lemont IL USA: ACM, 2021, pp. 1–7 DOI: 10.1145/3458744.3473366
  * [49] Suchithra Rajendran and A. Ravindran “Platelet ordering policies at hospitals using stochastic integer programming model and heuristic approaches to reduce wastage” In _Computers & Industrial Engineering_ 110, 2017, pp. 151–164 DOI: 10.1016/j.cie.2017.05.021
  * [50] Sergio Ruiz and Benjamín Hernández “A parallel solver for Markov decision process in crowd simulations” In _Proceedings of the Fourteenth Mexican International Conference on Artificial Intelligence (MICAI)_ , 2015, pp. 107–116 DOI: 10.1109/MICAI.2015.23
  * [51] Thomas J Sargent and John Stachurski “Dynamic programming on the GPU via JAX - QuantEcon Notes”, 2022 URL: https://notes.quantecon.org/submission/622ed4daf57192000f918c61/comments
  * [52] Lawrence V Snyder and Zuo-Jun Shen “Fundamentals of supply chain theory” Wiley, Hoboken, NJ, USA, 2019
  * [53] Ganid Srimool, Putchong Uthayopas and Juta Pichitlamkhen “Speeding up a large logistics optimization problems using GPU technology” In _The 8th Electrical Engineering/Electronics, Computer, Telecommunications and Information Technology (ECTI) Association of Thailand - Conference 2011_ , 2011, pp. 450–454 DOI: 10.1109/ECTICON.2011.5947872
  * [54] Shiaw Y. Su and Rolf A. Deininger “Generalization of White’s method of successive approximations to periodic Markovian decision processes” In _Operations Research_ 20.2, 1972, pp. 318–326 DOI: 10.1287/opre.20.2.318
  * [55] Ruoying Sun, Peng Sun, Jiahui Li and Gang Zhao “Inventory cost control model for fresh product retailers based on DQN” In _Proceedings of the 2019 IEEE International Conference on Big Data (Big Data)_ , 2019, pp. 5321–5325 DOI: 10.1109/BigData47090.2019.9006424
  * [56] Richard S. Sutton and Andrew G. Barto “Reinforcement learning: an introduction” Cambridge, MA, USA: The MIT Press, Cambridge, MA, USA, 2018
  * [57] TechnicalCity “A100 PCIe 40 GB vs GeForce RTX 3060” URL: https://technical.city/en/video/GeForce-RTX-3060-vs-A100-PCIe-40-GB
  * [58] United Nations “The Sustainable Development Goals Report 2022”, 2022 URL: https://unstats.un.org/sdgs/report/2022/The-Sustainable-Development-Goals-Report-2022.pdf

##  Appendix A Additional information for Scenario A

###  A.1 Scenario description

In this appendix we recast the problem formulated by [15] into a consistent
notation used for all three of the scenarios.

The state of the system, $S_{t}$, comprises two components: the orders in
transit $\underline{\text{O}_{t}}$ and the units in stock
$\underline{\text{X}_{t}}$:

$\displaystyle\underline{\text{O}_{t}}$
$\displaystyle=\left[O_{L-1,t}=A_{t-1},O_{L-2,t},...,O_{1,t}\right]$ (12)
$\displaystyle\underline{\text{X}_{t}}$
$\displaystyle=\left[X_{m,t}=O_{1,t-1},X_{m-1,t},...,X_{1,t}\right]$ (13)

for a total of $(m+L-1)$ elements, with lead time $L\geq 1$. The total number
of possible states is therefore $(A_{\max}+1)^{m+L-1}$. In the state
$S_{t}=\left[\underline{\text{O}_{t}},\underline{\text{X}_{t}}\right]$, the
entries are ordered by ascending age: the first element is the order placed on
day $t-1$ and the last element is the stock that will expire at the end of the
current day. The total number of units in stock at the start of day $t$ is
$X_{t}=\sum_{i=1}^{m}X_{i,t}$, the total number of units in transit at the
start of day $t$ is $O_{t}=\sum_{i=1}^{L-1}O_{i,t}$. The total number of units
in stock or in transit at the start of day $t$ is $I_{t}=X_{t}+O_{t}$. If
$L=1$ there is no in transit component to the state, and the first element of
$\underline{\text{X}_{t}}$ is $A_{t-1}$. In Table 7 we present the parameter
values that are the same for all of the experiments for Scenario A.

|  $D_{\max}$ |  $A_{\max}$ |  $C_{v}$ |  $C_{s}$ |  $C_{h}$ |  $\mu$ |  $\frac{\mu}{\sigma}$ |  $\gamma$ |  $\epsilon$
---|---|---|---|---|---|---|---|---|---
Value |  100 |  10 |  3 |  5 |  1 |  4 |  0.5 |  0.99 |  $1\times 10^{-4}$
Table 7: Parameter values that are consistent for all of the experiments for
Scenario A.

Daily demand, the stochastic element in the transition, is modelled using a
truncated gamma distribution. It does not depend on the state or the action.
The demand for the product is discrete and the gamma distribution is
continuous, so the probability that the daily demand is equal to
$d\in\\{0,1,...,D_{\max}\\}$ is:

$\displaystyle\text{Prob}(\Omega=\omega|S=s,A=a)$ $\displaystyle=P(\Omega=d)$
(14) $\displaystyle=P(D=d)$
$\displaystyle=\begin{cases}F(d+\frac{1}{2};\mu,\frac{\mu}{\sigma})-F(d-\frac{1}{2};\mu,\frac{\mu}{\sigma}),&\text{if
$d\in\\{0,1,...,D_{\max}-1$\\}}\\\
1-F(D_{\max}-\frac{1}{2};\mu,\frac{\mu}{\sigma}),&\text{if
$d=D_{\max}$}\end{cases}$

where $F(x;\mu,\frac{\mu}{\sigma})$ is the cumulative distribution function of
the gamma distribution parameterised by mean $\mu$ and coefficient of
variation $\frac{\mu}{\sigma}$, and $F(x;\mu,\frac{\mu}{\sigma})=0$ when
$x\leq 0$.

The reward function comprises four components: a holding cost per unit in
stock at the end of the period ($C_{h}$), a variable ordering cost per unit
($C_{v}$), a shortage cost per unit of unmet demand ($C_{s}$) and a wastage
cost per unit that perishes at the end of the period ($C_{w}$). The single-
step reward after taking action $A_{t}$ in state $S_{t}$ with
$\Omega_{t}=\left(D_{t}\right)$ is:

$R_{t+1}=-C_{v}A_{t}-C_{h}\left[X_{t}-D_{t}-W_{t}\right]^{+}-C_{s}\left[D_{t}-X_{t}\right]^{+}-C_{w}W_{t}$
(15)

where $W_{t}$ is the number of units that expire at the end of period $t$.

Equations 16 and 17 set out how the number of expired units, $W_{t}$, is
calculated and how the elements of $\underline{\text{X}_{t}}$ are updated when
following a FIFO issuing policy and a LIFO issuing policy, respectively.

$\displaystyle W_{t}$ $\displaystyle=[X_{1,t}-D_{t}]^{+}$ (16) $\displaystyle
X_{j,t+1}$
$\displaystyle=\left[X_{j+1,t}-\left[D_{t}-\sum_{k=1}^{j}X_{k,t}\right]^{+}\right]^{+}\quad\forall
j\in\\{1,2,...,m-1\\}$ $\displaystyle X_{m,t+1}$
$\displaystyle=O_{1,t}=A_{t-L+1}$

$\displaystyle W_{t}$
$\displaystyle=\left[X_{1,t}-\left[D_{t}-\sum^{m}_{k=2}X_{k,t}\right]^{+}\right]^{+}$
(17) $\displaystyle X_{j,t+1}$
$\displaystyle=\left[X_{j+1,t}-\left[D_{t}-\sum_{k=j+2}^{m}X_{k,t}\right]^{+}\right]^{+}\quad\forall
j\in\\{1,2,...,m-1\\}$ $\displaystyle X_{m,t+1}$
$\displaystyle=O_{1,t}=A_{t-L+1}$

The scenario is an infinite horizon MDP with a discount factor and no
periodicity in the state space. We therefore used a standard convergence test
for the value function [56], evaluating:

$\max_{s\in\mathbb{S}}|V_{i}(s)-V_{i-1}(s)|<\epsilon$ (18)

after each iteration. The inequality tests for the convergence of the values
themselves, and requires more iterations than the convergence tests used for
the other scenarios which are testing for convergence of the change in value
for each state. The test compares the current estimate of the value function
with the estimate from the immediately preceding iteration and does not
require previous checkpoints for evaluation. Therefore, to save storage space
and writing time, we saved a checkpoint every 100 iterations.

###  A.2 Additional results

We present additional results for Scenario A in Table 8: the order-up-to level
parameter $\texttt{S}_{\text{best}}$ fit using simulation optimization and the
mean and standard deviation of three KPIs calculated over 10,000 evaluation
rollouts for each policy.

| Service level (%) | Wastage (%) | Holding (units)
---|---|---|---
$m$ | Exp | $\texttt{S}_{\text{best}}$ | VI | SO | VI | SO | VI | SO
2 | 1 | 5 |  61.0 $\pm$ 1.4 |  58.6 $\pm$ 1.3 |  2.4 $\pm$ 0.6 |  2.2 $\pm$ 0.6 |  0.2 $\pm$ 0.0 |  0.2 $\pm$ 0.0
| 2 | 7 |  72.7 $\pm$ 1.6 |  76.6 $\pm$ 1.5 |  0.7 $\pm$ 0.4 |  1.5 $\pm$ 0.5 |  0.5 $\pm$ 0.1 |  0.8 $\pm$ 0.1
| 3 | 5 |  61.0 $\pm$ 1.4 |  58.6 $\pm$ 1.3 |  2.4 $\pm$ 0.6 |  2.2 $\pm$ 0.6 |  0.2 $\pm$ 0.0 |  0.2 $\pm$ 0.0
| 4 | 6 |  71.7 $\pm$ 1.6 |  68.6 $\pm$ 1.5 |  0.7 $\pm$ 0.3 |  0.7 $\pm$ 0.3 |  0.5 $\pm$ 0.1 |  0.5 $\pm$ 0.0
| 5 | 7 |  61.0 $\pm$ 1.4 |  55.4 $\pm$ 1.3 |  2.4 $\pm$ 0.6 |  2.4 $\pm$ 0.7 |  0.2 $\pm$ 0.0 |  0.2 $\pm$ 0.0
| 6 | 9 |  73.5 $\pm$ 1.7 |  69.4 $\pm$ 1.5 |  0.9 $\pm$ 0.4 |  1.1 $\pm$ 0.4 |  0.6 $\pm$ 0.1 |  0.6 $\pm$ 0.1
| 7 | 7 |  61.0 $\pm$ 1.4 |  55.4 $\pm$ 1.3 |  2.4 $\pm$ 0.6 |  2.4 $\pm$ 0.7 |  0.2 $\pm$ 0.0 |  0.2 $\pm$ 0.0
| 8 | 9 |  72.3 $\pm$ 1.6 |  69.4 $\pm$ 1.5 |  0.8 $\pm$ 0.4 |  1.1 $\pm$ 0.4 |  0.6 $\pm$ 0.1 |  0.6 $\pm$ 0.1
3 | 1 | 6 |  69.5 $\pm$ 1.5 |  68.3 $\pm$ 1.4 |  1.3 $\pm$ 0.4 |  1.4 $\pm$ 0.4 |  0.5 $\pm$ 0.1 |  0.5 $\pm$ 0.0
| 2 | 8 |  79.3 $\pm$ 1.5 |  83.3 $\pm$ 1.4 |  0.1 $\pm$ 0.1 |  0.2 $\pm$ 0.2 |  0.9 $\pm$ 0.1 |  1.3 $\pm$ 0.1
| 3 | 6 |  65.2 $\pm$ 1.4 |  68.3 $\pm$ 1.4 |  0.7 $\pm$ 0.3 |  1.4 $\pm$ 0.4 |  0.4 $\pm$ 0.0 |  0.5 $\pm$ 0.0
| 4 | 8 |  79.3 $\pm$ 1.5 |  83.3 $\pm$ 1.4 |  0.1 $\pm$ 0.1 |  0.2 $\pm$ 0.2 |  0.9 $\pm$ 0.1 |  1.3 $\pm$ 0.1
| 5 | 8 |  65.6 $\pm$ 1.6 |  62.6 $\pm$ 1.4 |  1.7 $\pm$ 0.5 |  1.4 $\pm$ 0.5 |  0.3 $\pm$ 0.0 |  0.4 $\pm$ 0.0
| 6 | 10 |  78.1 $\pm$ 1.6 |  75.5 $\pm$ 1.5 |  0.1 $\pm$ 0.1 |  0.1 $\pm$ 0.1 |  0.9 $\pm$ 0.1 |  0.9 $\pm$ 0.1
| 7 | 8 |  65.6 $\pm$ 1.6 |  62.6 $\pm$ 1.4 |  1.7 $\pm$ 0.5 |  1.4 $\pm$ 0.5 |  0.3 $\pm$ 0.0 |  0.4 $\pm$ 0.0
| 8 | 10 |  77.9 $\pm$ 1.6 |  75.5 $\pm$ 1.5 |  0.1 $\pm$ 0.1 |  0.1 $\pm$ 0.1 |  0.9 $\pm$ 0.1 |  0.9 $\pm$ 0.1
4 | 1 | 7 |  74.4 $\pm$ 1.4 |  76.4 $\pm$ 1.5 |  0.7 $\pm$ 0.3 |  1.5 $\pm$ 0.4 |  0.7 $\pm$ 0.1 |  0.8 $\pm$ 0.1
| 2 | 8 |  79.3 $\pm$ 1.5 |  83.3 $\pm$ 1.4 |  0.0 $\pm$ 0.0 |  0.0 $\pm$ 0.0 |  0.9 $\pm$ 0.1 |  1.3 $\pm$ 0.1
| 3 | 6 |  73.7 $\pm$ 1.4 |  68.5 $\pm$ 1.4 |  0.6 $\pm$ 0.3 |  0.5 $\pm$ 0.3 |  0.7 $\pm$ 0.1 |  0.5 $\pm$ 0.1
| 4 | 8 |  79.3 $\pm$ 1.5 |  83.3 $\pm$ 1.4 |  0.0 $\pm$ 0.0 |  0.0 $\pm$ 0.0 |  0.9 $\pm$ 0.1 |  1.3 $\pm$ 0.1
| 5 | 9 |  69.5 $\pm$ 1.5 |  69.3 $\pm$ 1.5 |  1.0 $\pm$ 0.4 |  1.0 $\pm$ 0.4 |  0.5 $\pm$ 0.1 |  0.6 $\pm$ 0.1
| 6 | 10 |  78.9 $\pm$ 1.7 |  75.5 $\pm$ 1.5 |  0.0 $\pm$ 0.0 |  0.0 $\pm$ 0.0 |  1.0 $\pm$ 0.1 |  0.9 $\pm$ 0.1
| 7 | 9 |  68.7 $\pm$ 1.5 |  69.3 $\pm$ 1.5 |  0.9 $\pm$ 0.4 |  1.0 $\pm$ 0.4 |  0.5 $\pm$ 0.1 |  0.6 $\pm$ 0.1
| 8 | 10 |  78.9 $\pm$ 1.7 |  75.5 $\pm$ 1.5 |  0.0 $\pm$ 0.0 |  0.0 $\pm$ 0.0 |  1.0 $\pm$ 0.1 |  0.9 $\pm$ 0.1
5 | 1 | 7 |  76.3 $\pm$ 1.5 |  76.6 $\pm$ 1.5 |  0.4 $\pm$ 0.2 |  0.7 $\pm$ 0.3 |  0.8 $\pm$ 0.1 |  0.8 $\pm$ 0.1
| 2 | 8 |  79.3 $\pm$ 1.5 |  83.3 $\pm$ 1.4 |  0.0 $\pm$ 0.0 |  0.0 $\pm$ 0.0 |  0.9 $\pm$ 0.1 |  1.3 $\pm$ 0.1
| 3 | 7 |  75.6 $\pm$ 1.4 |  76.6 $\pm$ 1.5 |  0.3 $\pm$ 0.2 |  0.7 $\pm$ 0.3 |  0.8 $\pm$ 0.1 |  0.8 $\pm$ 0.1
| 4 | 8 |  79.3 $\pm$ 1.5 |  83.3 $\pm$ 1.4 |  0.0 $\pm$ 0.0 |  0.0 $\pm$ 0.0 |  0.9 $\pm$ 0.1 |  1.3 $\pm$ 0.1
| 5 | 9 |  71.9 $\pm$ 1.6 |  69.5 $\pm$ 1.5 |  0.6 $\pm$ 0.3 |  0.4 $\pm$ 0.3 |  0.6 $\pm$ 0.1 |  0.6 $\pm$ 0.1
| 6 | 10 |  78.9 $\pm$ 1.7 |  75.5 $\pm$ 1.5 |  0.0 $\pm$ 0.0 |  0.0 $\pm$ 0.0 |  1.0 $\pm$ 0.1 |  0.9 $\pm$ 0.1
| 7 | 9 |  71.5 $\pm$ 1.6 |  69.5 $\pm$ 1.5 |  0.5 $\pm$ 0.3 |  0.4 $\pm$ 0.3 |  0.6 $\pm$ 0.1 |  0.6 $\pm$ 0.1
| 8 | 10 |  78.9 $\pm$ 1.7 |  75.5 $\pm$ 1.5 |  0.0 $\pm$ 0.0 |  0.0 $\pm$ 0.0 |  1.0 $\pm$ 0.1 |  0.9 $\pm$ 0.1
Table 8: The best order-up-to level $\texttt{S}_{\text{best}}$, fit using
simulation optimization, and KPIs for policies fit using value iteration (VI)
and simulation optimization (SO) for all of the experimental settings for
Scenario A from [15].

##  Appendix B Additional information for Scenario B

###  B.1 Scenario description

In this appendix we recast the problem formulated by [27] into a consistent
notation used for all three of the scenarios.

The state of the environment, $S_{t}$ comprises two components, one for each
product type. In the combined state
$S_{t}=\left[\underline{\text{X}^{a}_{t}},\underline{\text{X}^{b}_{t}}\right]$,
the elements in each component are ordered by ascending age:

$\displaystyle\underline{\text{X}^{a}_{t}}$
$\displaystyle=[X^{a}_{m,t}=A^{a}_{t-1},X^{a}_{m-1,t},...,X^{a}_{1,t}]$ (19)
$\displaystyle\underline{\text{X}^{b}_{t}}$
$\displaystyle=[X^{b}_{m,t}=A^{b}_{t-1},X^{b}_{m-1,t},...,X^{b}_{1,t}]$ (20)

for a total number of $2m$ elements. The total number of possible states is
therefore $\left(A_{\max}^{a}+1\right)^{m}+\left(A_{\max}^{b}+1\right)^{m}$.
The total number of units in stock at the start of period $t$ is
$I^{a}_{t}=X^{a}_{t}=\sum_{i=1}^{m}X^{a}_{i,t}$ for product A and
$I^{b}_{t}=X^{b}_{t}=\sum_{i=1}^{m}X^{b}_{i,t}$ for product B. In Table 9 we
present the parameter values that are the same for all of the experiments for
Scenario B.

|  $C^{a}_{v}$ |  $C^{b}_{v}$ |  $C^{a}_{r}$ |  $C^{b}_{r}$ |  $\rho$ |  $\gamma$ |  $\epsilon$
---|---|---|---|---|---|---|---
Value |  0.5 |  0.5 |  1.0 |  1.0 |  0.5 |  1.0 |  $1\times 10^{-4}$
Table 9: Parameter values that are consistent for all of the experiments for
Scenario B.

The stochastic element of the transition is the number of products of each
type issued, ($H^{a},H^{b}$). The number of units of product B issued only
depends on the demand for product B and the total stock of product B, but the
number of units of product A that are issued depends on the demand for product
A, the total stock of product A, and any excess demand for product B for which
the customer is willing to accept product A.

Let the demand for product A be $D^{a}$, the demand for product B be $D^{b}$,
the excess demand for product B where the customer is willing to accept
product A be $D^{u}$ and the total demand for product A including any
substitution be $D^{z}=D^{a}+D^{u}$. To calculate the probability of a
combination $\omega=(h^{a},h^{b})$ given a particular state $s$, we consider
five possible cases:

$\displaystyle\text{Prob}\left(\Omega=\omega|S=s,A=a\right)$
$\displaystyle=P\left(\Omega=\left(h^{a},h^{b}\right)|S=s\right)$ (21)
$\displaystyle=P\left(H^{a}=h^{a},H^{b}=h^{b}|S=s\right)$
$\displaystyle=\begin{cases}0,&\text{if }h^{a}>I^{a}\text{ or }h^{b}>I^{b}\\\
P(D^{a}=h^{a})P(D^{b}=h^{b}),\phantom{onetwothreefourfivesixse}&\text{if
}h^{a}<I^{a}\text{ and }h^{b}<I^{b}\\\ P(D^{a}\geq
I^{a})P(D^{b}=h^{b}),&\text{if }h^{a}=I^{a}\text{ and }h^{b}<I^{b}\\\
P(D^{z}=h^{a}|S=s)P(D^{b}\geq I^{b}),&\text{if }h^{a}<I^{a}\text{ and
}h^{b}=I^{b}\\\ P(D^{z}\geq I^{a})P(D^{b}\geq I^{b}),&\text{if
}h^{a}=I^{a}\text{ and }h^{b}=I^{b}\\\ \end{cases}$
$\displaystyle=\begin{cases}0,&\text{if }h^{a}>I^{a}\text{ or }h^{b}>I^{b}\\\
P(h^{a};\mu^{a})P(h^{b};\mu^{b}),&\text{if }h^{a}<I^{a}\text{ and
}h^{b}<I^{b}\\\
\left[1-F\left(I^{a}-1;\mu^{a}\right)\right]P(h^{b};\mu^{b}),&\text{if
}h^{a}=I^{a}\text{ and }h^{b}<I^{b}\\\
P(D^{z}=h^{a}|S=s)\left[1-\left(F(I^{b}-1;\mu^{b}\right)\right],&\text{if
}h^{a}<I^{a}\text{ and }h^{b}=I^{b}\\\
\left[1-\sum_{d=0}^{I^{a}-1}P(D^{z}=d|S=s)\right]\left[1-\left(F(I^{b}-1;\mu^{b}\right)\right],&\text{if
}h^{a}=I^{a}\text{ and }h^{b}=I^{b}\\\ \end{cases}$

For the fourth and fifth cases there may be substitution, and therefore we
need to consider the distribution of the total demand for product A and the
distribution of the demand for substitution:

$\displaystyle P(D^{z}=d^{z}|S=s)$ $\displaystyle=P(D^{z}=d^{z}|I^{b}=y)$
$\displaystyle=\sum_{k=0}^{d^{z}}P(D^{a}=k)P(D^{u}=d^{z}-k|I^{b}=y,D^{b}\geq
y)$ (22)
$\displaystyle=\sum_{k=0}^{d^{z}}P(k;\mu^{a})P(D^{u}=d^{z}-k|I^{b}=y,D^{b}\geq
y)$

$\displaystyle P(D^{u}=d^{u}|I^{b}=y,D^{b}\geq y)$
$\displaystyle=\begin{cases}\sum_{c=0}^{\infty}P(D^{b}=c+y)(1-\rho)^{c},&\text{if
$d^{u}=0$}\\\ \sum_{c=d^{u}}^{\infty}P(D^{b}=c+y)P(d^{u};c,\rho),&\text{if
$d^{u}>0$}\\\ \end{cases}$ (23)
$\displaystyle=\begin{cases}\sum_{c=0}^{\infty}P(c+y;\mu^{b})(1-\rho)^{c},&\phantom{ac}\text{if
$d^{u}=0$}\\\
\sum_{c=d^{u}}^{\infty}P(c+y;\mu^{b})P(d^{u};c,\rho),&\phantom{ac}\text{if
$d^{u}>0$}\end{cases}$

where $P(x;c,\rho)$ is a binomial probability mass function representing the
probability that there are $x$ units of excess demand for product B willing to
accept product A out of a total of $c$ units of excess demand for product B
and the probability of being willing to accept the substitution is $\rho$.
$P(x;\mu^{a})$ and $P(x;\mu^{b})$ are the probability mass functions of
independent Poisson distributions for the daily demand of product A and B
parameterised by mean daily demands $\mu^{a}$ and $\mu^{b}$ respectively, and
$F(x;\mu^{a})$ and $F(x;\mu^{b})$ are the corresponding cumulative
distribution functions.

We calculated the values of $P(D^{u}=d^{u}|I^{b}=y,D^{b}\geq y)$ and
$P(D^{z}=d^{z}|S=s)$, for $d^{u}\in\\{0,1,...,D_{\max}\\}$ and
$d^{z}\in\\{0,1,...,D_{\max}\\}$, where
$D_{\max}=\left(\left(m\max(A^{a}_{\max},A^{b}_{\max})\right)+2\right)$, at
the start of value iteration following the MATLAB implementation of [27].

The reward function comprises two components which can be different for each
product: a variable ordering cost per unit ($C^{a}_{v}$, $C^{b}_{v}$) and
revenue per unit sold ($C^{a}_{r}$, $C^{b}_{r}$). The single step reward after
taking action $A_{t}$ in state $S_{t}$ with
${\Omega}_{t}=\left(H^{a}_{t},H^{b}_{t}\right)$ is:

$R_{t+1}=-\left(C^{a}_{v}A^{a}_{t}+C^{a}_{v}A^{b}_{t}\right)+\left(C^{a}_{r}H^{a}_{t}+C^{b}_{r}H^{b}_{t}\right)$
(24)

Equation 25 shows how the elements of $\underline{\text{X}^{a}_{t}}$ and
$\underline{\text{X}^{b}_{t}}$ are updated following a FIFO issuing policy.

$\displaystyle X^{a}_{j,t+1}$
$\displaystyle=X^{a}_{j+1,t}-\left[H^{a}_{t}-\sum_{k=1}^{j}X^{a}_{k,t}\right]^{+}\quad\forall
j\in\\{1,2,...,m-1\\}$ (25) $\displaystyle X^{a}_{m,t+1}$
$\displaystyle=A^{a}_{t}$ $\displaystyle X^{b}_{j,t+1}$
$\displaystyle=X^{b}_{j+1,t}-\left[H^{b}_{t}-\sum_{k=1}^{j}X^{b}_{k,t}\right]^{+}\quad\forall
j\in\\{1,2,...,m-1\\}$ $\displaystyle X^{b}_{m,t+1}$ $\displaystyle=A^{b}_{t}$

The maximum order quantities for value iteration, $A^{a}_{\max}$ and
$A^{b}_{\max}$ are calculated independently for each product following the
newsvendor model [52]:

$A^{k}_{\max}=\Bigg{\lceil}F^{-1}\left(\frac{C^{k}_{r}-C^{k}_{v}}{C^{k}_{r}};m\mu^{k}\right)\Bigg{\rceil}^{+},\quad\forall
k\in\\{a,b\\}$ (26)

where $F(x;m\mu^{k})$ is the cumulative distribution function of a Poisson
distribution parameterised by $m\mu^{k}$ and $\mu^{k}$ is the mean daily
demand for product $k$.

The maximum order quantities reported for experiments P1 to P4 in Table 3 of
[45] are not consistent with the number of states reported in that table. We
have assumed that they instead represent the number of actions (one higher
than the maximum order quantity, due to the possibility of ordering zero
units) as this is consistent with the number of states reported, with Table 1
of [45] and, for experiment P1, with the corresponding experiment in [27].

The initial estimate for the value function is the expected one-step ahead
sales revenue:

$V_{0}(s)=\sum^{I^{a}_{t}}_{h^{a}=0}\sum^{I^{b}_{t}}_{h^{b}=0}P(H^{a}=h^{a},H^{b}=h^{b}|S=s)(h^{a}C^{a}_{r}+h^{b}C^{b}_{r})$
(27)

We used the same convergence test as [27], evaluating:

$\max_{s\in\mathbb{S}}\left[V_{i}(s)-V_{i-1}(s)\right]-\min_{s\in\mathbb{S}}\left[V_{i}(s)-V_{i-1}(s)\right]<\epsilon$
(28)

after each iteration. The inequality tests for the convergence of the change
in value for each state. When the value of each state is changing by the same
amount, the best action for each state will not change and therefore the
estimate of the optimal policy is stable. We saved a checkpoint after every
iteration.

###  B.2 Additional results

We present additional results for the experimental settings from [27] in Table
10 and for the experimental settings from [45] in Table 11. In each table we
present the best combination of order-up-to level parameters
$\left(\texttt{S}^{a},\texttt{S}^{b}\right)_{\text{best}}$ fit using
simulation optimization and the mean and standard deviation of three KPIs
calculated over 10,000 evaluation rollouts for each policy. Demand for product
B was considered to be satisfied for the purposes of calculating the service
level if filled by product A when substitution was acceptable.

The only difference between experiment 1 from [27] and experiment P1 from [45]
is that value iteration was run for 100 iterations for experiment P1: more
than were required for convergence. The best parameters for the modified base-
stock policy and the KPIs from evaluating the policies are the same for these
experiments, as we would expect.

| Service level (%) | Wastage (%) | Holding (units)
---|---|---|---
$m$ | Exp | Product | $\texttt{S}_{\text{best}}$ | VI | SO | VI | SO | VI | SO
2 | 1 | A | 13 |  95.5 $\pm$ 0.8 |  95.2 $\pm$ 0.8 |  6.0 $\pm$ 1.0 |  6.3 $\pm$ 1.0 |  2.7 $\pm$ 0.1 |  2.7 $\pm$ 0.1
|  | B | 12 |  94.9 $\pm$ 0.8 |  95.5 $\pm$ 0.7 |  4.2 $\pm$ 0.8 |  5.3 $\pm$ 0.9 |  2.1 $\pm$ 0.1 |  2.3 $\pm$ 0.1
| 2 | A | 18 |  96.9 $\pm$ 0.6 |  96.6 $\pm$ 0.6 |  4.2 $\pm$ 0.8 |  4.4 $\pm$ 0.7 |  3.7 $\pm$ 0.2 |  3.6 $\pm$ 0.2
|  | B | 7 |  91.5 $\pm$ 1.1 |  92.5 $\pm$ 1.0 |  6.5 $\pm$ 1.2 |  8.3 $\pm$ 1.3 |  1.2 $\pm$ 0.1 |  1.3 $\pm$ 0.1
3 | 1 | A | 15 |  98.3 $\pm$ 0.5 |  98.3 $\pm$ 0.5 |  2.2 $\pm$ 0.6 |  2.4 $\pm$ 0.6 |  4.9 $\pm$ 0.2 |  4.8 $\pm$ 0.2
|  | B | 14 |  98.4 $\pm$ 0.4 |  98.4 $\pm$ 0.4 |  1.5 $\pm$ 0.5 |  1.7 $\pm$ 0.5 |  4.1 $\pm$ 0.2 |  4.1 $\pm$ 0.2
| 2 | A | 21 |  99.1 $\pm$ 0.3 |  99.2 $\pm$ 0.3 |  1.2 $\pm$ 0.4 |  1.3 $\pm$ 0.4 |  6.7 $\pm$ 0.3 |  6.7 $\pm$ 0.2
|  | B | 8 |  96.6 $\pm$ 0.7 |  96.1 $\pm$ 0.7 |  3.3 $\pm$ 0.9 |  3.2 $\pm$ 0.8 |  2.4 $\pm$ 0.1 |  2.3 $\pm$ 0.1
| 3 | A | 15 |  98.3 $\pm$ 0.5 |  98.3 $\pm$ 0.5 |  2.2 $\pm$ 0.6 |  2.4 $\pm$ 0.6 |  4.9 $\pm$ 0.2 |  4.8 $\pm$ 0.2
|  | B | 14 |  98.4 $\pm$ 0.4 |  98.4 $\pm$ 0.4 |  1.5 $\pm$ 0.5 |  1.7 $\pm$ 0.5 |  4.1 $\pm$ 0.2 |  4.1 $\pm$ 0.2
| 4 | A | 21 |  99.1 $\pm$ 0.3 |  99.2 $\pm$ 0.3 |  1.2 $\pm$ 0.4 |  1.3 $\pm$ 0.4 |  6.6 $\pm$ 0.3 |  6.7 $\pm$ 0.2
|  | B | 8 |  96.6 $\pm$ 0.7 |  96.1 $\pm$ 0.7 |  3.3 $\pm$ 0.9 |  3.2 $\pm$ 0.8 |  2.4 $\pm$ 0.1 |  2.3 $\pm$ 0.1
Table 10: The best combination of order-up-to levels $\left(\texttt{S}^{a},\texttt{S}^{b}\right)_{\text{best}}$, fit using simulation optimization, and KPIs for policies fit using value iteration (VI) and simulation optimization (SO) for all of the experimental settings for Scenario B from [27] | Service level (%) | Wastage (%) | Holding (units)
---|---|---|---
$m$ | Exp | Product | $\texttt{S}_{\text{best}}$ | VI | SO | VI | SO | VI | SO
2 | P1 | A | 13 |  95.5 $\pm$ 0.8 |  95.2 $\pm$ 0.8 |  6.0 $\pm$ 1.0 |  6.3 $\pm$ 1.0 |  2.7 $\pm$ 0.1 |  2.7 $\pm$ 0.1
|  | B | 12 |  94.9 $\pm$ 0.8 |  95.5 $\pm$ 0.7 |  4.2 $\pm$ 0.8 |  5.3 $\pm$ 0.9 |  2.1 $\pm$ 0.1 |  2.3 $\pm$ 0.1
| P2 | A | 13 |  95.5 $\pm$ 0.8 |  95.1 $\pm$ 0.8 |  6.0 $\pm$ 1.0 |  6.2 $\pm$ 0.9 |  2.7 $\pm$ 0.1 |  2.6 $\pm$ 0.1
|  | B | 14 |  95.9 $\pm$ 0.6 |  95.5 $\pm$ 0.7 |  3.5 $\pm$ 0.7 |  3.7 $\pm$ 0.7 |  2.5 $\pm$ 0.1 |  2.5 $\pm$ 0.1
| P3 | A | 16 |  96.2 $\pm$ 0.7 |  96.8 $\pm$ 0.6 |  4.9 $\pm$ 0.9 |  6.0 $\pm$ 0.9 |  3.2 $\pm$ 0.2 |  3.4 $\pm$ 0.1
|  | B | 14 |  95.9 $\pm$ 0.6 |  95.8 $\pm$ 0.6 |  3.5 $\pm$ 0.7 |  3.7 $\pm$ 0.7 |  2.5 $\pm$ 0.1 |  2.5 $\pm$ 0.1
| P4 | A | 18 |  96.8 $\pm$ 0.6 |  96.7 $\pm$ 0.6 |  4.2 $\pm$ 0.8 |  4.5 $\pm$ 0.8 |  3.7 $\pm$ 0.2 |  3.7 $\pm$ 0.2
|  | B | 17 |  96.5 $\pm$ 0.6 |  97.1 $\pm$ 0.5 |  2.9 $\pm$ 0.6 |  3.8 $\pm$ 0.7 |  2.9 $\pm$ 0.2 |  3.2 $\pm$ 0.1
Table 11: The best combination of order-up-to levels
$\left(\texttt{S}^{a},\texttt{S}^{b}\right)_{\text{best}}$, fit using
simulation optimization, and KPIs for policies fit using value iteration (VI)
and simulation optimization (SO) for all of the experimental settings for
Scenario B from [45]

##  Appendix C Additional information for Scenario C

###  C.1 Scenario description

In this appendix we recast the problem formulated by [40] into a consistent
notation used for all three of the scenarios.

The state of the environment, $S_{t}$, comprises two components:
$\tau\in\\{0,1,...,6\\}$, representing the day of the week, and the units in
stock at the start of the day
$\underline{\text{X}_{t}}=\left[X_{m-1,t},X_{m-2,t},...,X_{1,t}\right]$. The
lead time, $L$, is always zero which means that the units ordered on day $t$
are received before any demand arises on day $t$. There are therefore only
$m-1$ elements in $\underline{\text{X}_{t}}$ and a total of $m$ elements,
including $\tau$, in $S_{t}$.

In the previous problems, the maximum value of each element of
$\underline{\text{X}_{t}}$ was $A_{\max}$, because all units arrived with the
same remaining useful life. All units received in the same period would be in
the same element of $\underline{\text{X}_{t}}$. In this scenario, the
remaining useful life on arrival is stochastic and therefore, depending on the
policy, a series of orders could be received such that an element of
$\underline{\text{X}_{t}}$ would exceed $A_{\max}$. We assume there is a
maximum capacity of $A_{\max}$ for stock of each possible value of remaining
useful life. Units received in excess of this limit are not accepted at the
point of delivery. The total number of possible states is therefore
$7\times\left(A_{\max}+1\right)^{m-1}$. The entries in
$\underline{\text{X}_{t}}$ are ordered by ascending age: the first element
represents stock with $m-1$ days before expiry, and the last element is the
stock that will expire at the end of day $t$. In Table 12 we present the
parameter values that are the same for all of the experiments for Scenario C.

|  $D_{\max}$ |  $A_{\max}$ |  $C_{f}$ |  $C_{h}$ |  $C_{s}$ |  $C_{w}$ |  $\gamma$ |  $\epsilon$
---|---|---|---|---|---|---|---|---
Value |  20 |  20 |  10 |  1 |  20 |  5 |  0.95 |  $1\times 10^{-4}$
Table 12: Parameter values that are consistent for all of the experiments for
Scenario C.

The stochastic elements in the transition are the daily demand $D$, and the
age profile of the units received:
$\underline{\text{Y}}=\left[Y_{m},Y_{m-1},...,Y_{1}\right]$.

The probability of a given random outcome $\omega$ is the product of the
probability of the demand given the state, and the probability of receiving
units with a specific age profile given the action:

$\displaystyle\text{Prob}\left(\Omega=\omega|S=s,A=a\right)$
$\displaystyle=P\left(\Omega=\left(d,\underline{\text{y}}\right)|S=s,A=a\right)$
(29)
$\displaystyle=P\left(D=d,\underline{\text{Y}}=\underline{\text{y}}|S=s,A=a\right)$
$\displaystyle=P\left(D=d|S=s\right)P\left(\underline{\text{Y}}=\underline{\text{y}}|A=a\right)$

Demand is modelled by truncated negative binomial distributions, one for each
day of the week. The demand distribution therefore only depends on the weekday
element of the state. The negative binomial distribution models the number of
failures, $x$, in a series of repeated Bernoulli trials before achieving a
specified number of successes. The probability that daily demand is equal to
$d$ on weekday $\tau$ is:

$\displaystyle P(D=d|S=s)=\begin{cases}P(d;n^{\tau},\delta^{\tau}),&\text{if
}d\in\\{0,1,...,D_{\max}-1\\}\\\
1-F(D_{\max}-1;n^{\tau},\delta^{\tau}),&\text{if }d=D_{\max}\\\ \end{cases}$
(30)

where $P(x;n^{\tau},\delta^{\tau})$ is the probability mass function of a
negative binomial distribution parameterised by a target number of successes
$n^{\tau}$ and a mean $\delta^{\tau}$ for weekday $\tau$ and
$F(x;n^{\tau},\delta^{\tau})$ is the corresponding cumulative distribution
function. The probability of success in an individual Bernoulli trial is
$p^{\tau}=\frac{n^{\tau}}{n^{\tau}+\delta^{\tau}}$. The parameters for each
day of the week are set out in Table 13.

$\tau$ |  0 |  1 |  2 |  3 |  4 |  5 |  6
---|---|---|---|---|---|---|---
$n^{\tau}$ |  3.5 |  11.0 |  7.2 |  11.1 |  5.9 |  5.5 |  2.2
$\delta^{\tau}$ |  5.7 |  6.9 |  6.5 |  6.2 |  5.8 |  3.3 |  3.4
Table 13: Parameters of the demand distribution for each weekday from Monday
($\tau=0$) to Sunday ($\tau=6$)

The remaining useful life of units on arrival is modelled by a multinomial
distribution with a number of trials equal to the order quantity $a$ and a
number of events equal to the maximum useful life $m$. The probability mass
function for the distribution is:

$\displaystyle P(\underline{\text{Y}}=\underline{\text{y}}|A=a)$
$\displaystyle=P(Y_{m}=y_{m},...,Y_{1}=y_{1}|A=a)$ (31)
$\displaystyle=\begin{cases}\frac{a!}{y_{m}!y_{m-1}!...y_{1}!}p_{m}(a)^{y_{m}}p_{m-1}(a)^{y_{m-1}}...p_{1}(a)^{y_{1}},&\text{if
$a=\sum_{i=1}^{m}y_{i}$}\\\ 0,&\text{if $a\neq\sum_{i=1}^{m}y_{i}$}\\\
\end{cases}$

The parameters of the multinomial distribution are modelled by an affine
function of the order quantity $a$:

$\log\left(\frac{p_{k}(a)}{p_{1}(a)}\right)=c_{0}^{k}+c^{k}_{1}a,\quad\forall
k\in\\{2,3,..,m\\}$ (32)

If the distribution of remaining useful life on arrival does not depend on
order quantity, and therefore the uncertainty is exogenous,
$c_{1}^{k}=0\quad\forall k\in\\{2,3,..,m\\}$. The values of $c_{0}^{k}$ and
$c_{1}^{k}$ for our experiments are set out in Table 14. These represent a
subset of the experiments run by [40]. The parameters for the two experiments
where $m=5$ were determined by [40] by fitting multinomial logistic regression
models to observed data from a hospital system in Ontario, Canada.

$m$ | Exp | $c_{0}^{2}$ | $c_{0}^{3}$ | $c_{0}^{4}$ | $c_{0}^{5}$ | $c_{0}^{6}$ | $c_{0}^{6}$ | $c_{0}^{8}$ | $c_{1}^{2}$ | $c_{1}^{3}$ | $c_{1}^{4}$ | $c_{1}^{5}$ | $c_{1}^{6}$ | $c_{1}^{7}$ | $c_{1}^{8}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
3 | 1 | 1.0 | 0.5 |  |  |  |  |  |  |  |  |  |  |  |
| 2 | 1.0 | 0.5 |  |  |  |  |  | 0.40 | 0.80 |  |  |  |  |
5 | 1 | 1.6 | 2.6 | 2.8 | 1.6 |  |  |  |  |  |  |  |  |  |
| 2 | 1.9 | 3.1 | 3.1 | 2.5 |  |  |  | -0.03 | -0.06 | -0.03 | -0.09 |  |  |
8 | 1 | 0.8 | 1.4 | 1.9 | 2.3 | 1.7 | 1.2 | 0.8 |  |  |  |  |  |  |
| 2 | 0.8 | 1.4 | 1.9 | 2.3 | 1.7 | 1.2 | 0.8 | -0.03 | -0.04 | -0.05 | -0.06 | -0.07 | -0.08 | -0.09

Table 14: Parameters for the affine function used to model the parameters of
the multinomial distribution of remaining useful life on arrival for each
experiment. These are a subset of the experiments described by [40]

.

The reward function comprises four components: a holding cost per unit in
stock at the end of the period ($C_{h}$), a shortage cost per unit of unmet
demand ($C_{s}$), a wastage cost per unit that perishes at the end of the
period ($C_{w}$) and a fixed ordering cost which is incurred when $A_{t}>0$
($C_{f}$). The single-step reward function after taking action $A_{t}$ in
state $S_{t}=\left(\tau_{t},\underline{\text{X}_{t}}\right)$, and observing
$\Omega_{t}=(D_{t},\underline{\text{Y}_{t}})$ is

$R_{t+1}=-C_{f}\mathbbm{1}_{A_{t}>0}-C_{h}\left[Y_{m,t}+\sum^{m-1}_{i=1}\min\left(X_{i,t}+Y_{i,t},A_{\max}\right)-D_{t}\right]^{+}\\\
-C_{s}\left[D_{t}-Y_{m,t}-\sum^{m-1}_{i=1}\min\left(X_{i,t}+Y_{i,t},A_{\max}\right)\right]^{+}-C_{w}\left[\min\left(X_{1,t}+Y_{1,t},A_{\max}\right)-D_{t}\right]^{+}$
(33)

Equation 34 shows how the elements of the state are updated, following a OUFO
policy:

$\displaystyle\tau_{t+1}$ $\displaystyle=(\tau_{t}+1)\mod 7$ (34)
$\displaystyle X_{j,t+1}$
$\displaystyle=\left[\min\left(X_{j+1,t}+Y_{j+1,t},A_{\max}\right)-\left[D_{t}-\sum^{j}_{k=1}\min\left(X_{k,t}+Y_{k,t},A_{\max}\right)\right]^{+}\right]^{+}\quad\forall
j\in\\{1,...,m-2\\}$ $\displaystyle X_{m-1,t+1}$
$\displaystyle=\left[Y_{m,t}-\left[D_{t}-\sum^{m-1}_{k=1}\min\left(X_{k,t}+Y_{k,t},A_{\max}\right)\right]^{+}\right]^{+}$

The scenario is an infinite horizon MDP with a discount factor and periodicity
because the demand depends on the day of the week. We take advantage of the
periodicity of the problem, and use a convergence test based on those
described by [54]. Performing this convergence test requires retaining at
least the last seven (as this is the periodicity) estimates of the value
function, and we can only test for convergence after we have run at least
seven iterations. We tested the following inequality at the end of each
iteration once $i\geq 7$:

$\displaystyle\Delta_{\text{max},i}$
$\displaystyle=\max_{s\in\mathbb{S}}\left[\sum_{j=0}^{6}\frac{1}{\gamma^{i-j-1}}\left(V_{i-j}(s)-V_{i-j-1}(s)\right)\right]$
(35) $\displaystyle\Delta_{\text{min},i}$
$\displaystyle=\min_{s\in\mathbb{S}}\left[\sum_{j=0}^{6}\frac{1}{\gamma^{i-j-1}}\left(V_{i-j}(s)-V_{i-j-1}(s)\right)\right]$
(36) $\displaystyle\Delta_{\text{max},i}-\Delta_{\text{min},i}$
$\displaystyle\leq
2\epsilon\min\left[|\Delta_{\text{max},i}|,|\Delta_{\text{min},i}|\right]$
(37)

When the inequality is met, the additional undiscounted reward being added to
each state in one whole cycle (one week) is approximately the same, subject to
our confidence level. In turn this means that for each weekday, every state is
being increased by the same amount and therefore the best action for each
state will not change. We therefore terminated value iteration when the
inequality was met. If there were no discounting, and so $\gamma=1$, the term
in the square brackets would be equal to $V_{i}(s)-V_{i-7}(s)$: the total
change in value from one cycle (in this case, one week). This convergence test
relies on checkpoints from previous iterations, and we therefore saved a
checkpoint every iteration.

###  C.2 Additional results

In Table 15 we present the best combination of parameters for the heuristic
policy fit using simulation optimization. In Table 16 we present the mean and
standard deviation of three KPIs calculated over 10,000 evaluation rollouts
for each policy. For consistency with the calculation of the reward function
in [40] the holding KPI includes the units that will expire at the end of the
current day. These units are excluded from the calculation of the stock
holding at the end of the day in the other scenarios

|  |  |  Weekday $\tau$
---|---|---|---
$m$ | Exp | Parameter | 0 | 1 | 2 | 3 | 4 | 5 | 6
3 | 1 | $\texttt{S}_{\text{best}}$ | 13 | 12 | 14 | 11 | 11 | 8 | 7
|  | $\texttt{s}_{\text{best}}$ | 6 | 7 | 7 | 6 | 6 | 3 | 3
| 2 | $\texttt{S}_{\text{best}}$ | 14 | 14 | 15 | 13 | 12 | 9 | 9
|  | $\texttt{s}_{\text{best}}$ | 7 | 7 | 7 | 7 | 6 | 3 | 4
5 | 1 | $\texttt{S}_{\text{best}}$ | 16 | 17 | 16 | 13 | 13 | 10 | 14
|  | $\texttt{s}_{\text{best}}$ | 7 | 8 | 8 | 7 | 7 | 3 | 3
| 2 | $\texttt{S}_{\text{best}}$ | 17 | 16 | 16 | 13 | 13 | 11 | 14
|  | $\texttt{s}_{\text{best}}$ | 7 | 7 | 9 | 8 | 8 | 3 | 4
8 | 1 | $\texttt{S}_{\text{best}}$ | 19 | 15 | 18 | 18 | 14 | 13 | 16
|  | $\texttt{s}_{\text{best}}$ | 8 | 8 | 8 | 7 | 8 | 3 | 4
| 2 | $\texttt{S}_{\text{best}}$ | 18 | 18 | 16 | 15 | 14 | 11 | 15
|  | $\texttt{s}_{\text{best}}$ | 9 | 8 | 9 | 7 | 7 | 3 | 4
Table 15: The best combination of parameters for the heuristic policy $\left(\left(\texttt{s}^{0},\texttt{S}^{0}\right),...,\left(\texttt{s}^{6},\texttt{S}^{6}\right)\right)_{\text{best}}$ fit using simulation optimization for each of our experiments for Scenario C, a subset of the experiments run by [40]. |  | Service level (%) | Wastage (%) | Holding (units)
---|---|---|---|---
$m$ | Exp | VI | SO | VI | SO | VI | SO
3 | 1 |  95.3 $\pm$ 0.9 |  95.3 $\pm$ 0.9 |  12.6 $\pm$ 1.3 |  12.6 $\pm$ 1.4 |  4.9 $\pm$ 0.1 |  4.9 $\pm$ 0.1
| 2 |  96.6 $\pm$ 0.8 |  96.2 $\pm$ 0.8 |  7.0 $\pm$ 1.1 |  7.2 $\pm$ 1.1 |  5.8 $\pm$ 0.1 |  5.7 $\pm$ 0.1
5 | 1 |  97.4 $\pm$ 0.7 |  97.0 $\pm$ 0.7 |  3.2 $\pm$ 0.8 |  3.0 $\pm$ 0.7 |  6.8 $\pm$ 0.1 |  6.7 $\pm$ 0.1
| 2 |  97.4 $\pm$ 0.7 |  97.5 $\pm$ 0.7 |  3.1 $\pm$ 0.7 |  3.4 $\pm$ 0.8 |  6.8 $\pm$ 0.1 |  7.0 $\pm$ 0.2
8 | 1 | — |  97.9 $\pm$ 0.6 | — |  0.7 $\pm$ 0.3 | — |  8.0 $\pm$ 0.2
| 2 | — |  97.7 $\pm$ 0.6 | — |  1.0 $\pm$ 0.4 | — |  7.5 $\pm$ 0.2
Table 16: KPIs for policies fit using value iteration (VI) and simulation
optimization (SO) for each of our experiments for Scenario C, a subset of the
experiments run by [40]. Value iteration was not feasible when $m=8$.

##  Appendix D Notation

In Table 17 we summarise the notation we have used to recast the problems
described in Scenarios A, B and C. In Scenario B we use a superscript $a$ for
product A and $b$ for product B if a variable is product-specific. In Scenario
C we use a superscript $\tau$ (or weekday index from 0 to 6 representing
Monday to Sunday, respectively) if a variable is weekday-specific. We drop the
subscript $t$ in some contexts where elements from different days do not
feature.

Markov decision process | $\mathbb{S}$ | Set of possible states
---|---|---
| $S_{t}$ |  State observed at the start of day $t$
| $s$ |  A specific element of $\mathbb{S}$
| $\mathbb{A}$ | Set of possible actions
| $A_{t}$ |  Action taken at the start of day $t$ after observing $S_{t}$
| $a$ |  A specific element of $\mathbb{A}$
| $\mathbb{\Psi}$ | Set of possible rewards
| $R_{t}$ |  Reward received when state $S_{t}$ is observed
| $r$ |  A specific element of $\mathbb{\Psi}$
| $\mathbb{\Omega}$ | Set of possible realisations of stochastic elements in a transition
| $\Omega_{t}$ |  Realisation of the stochastic elements in the transition between $S_{t}$ and $S_{t+1}$
| $\omega$ |  A specific element of $\mathbb{\Omega}$
| $\gamma$ | Discount factor
| $G_{t}$ |  Return, the discounted sum of rewards received after taking action $A_{t}$
| $\pi(s)$ | Policy, a function mapping a state to an action
| $\pi^{*}(s)$ | Optimal policy, policy with the maximum expected return from every state
| $V^{\pi}(s)$ |  Value function, expected return starting in state $s$ and following policy $\pi$
| $Q^{\pi}(s,a)$ |  State-action value function, expected return taking action $a$ in state $s$ and
|  |  following policy $\pi$ thereafter
| $T(s,a,\omega)$ | Deterministic transition function
Reward function components | $C_{v}$ | Variable ordering cost per unit
| $C_{f}$ | Fixed ordering cost
| $C_{w}$ | Wastage cost per unit
| $C_{s}$ | Shortage cost per unit
| $C_{h}$ | Holding cost per unit
| $C_{r}$ | Revenue per unit
Value iteration | $\epsilon$ | Tolerance for convergence test
Heuristic policy parameters | S | Order-up-to level
| s | Reorder point
Inventory control | $L$ | Lead time
| $m$ | Maximum useful life
| $D_{t}$ |  Demand on day $t$
| $D_{\max}$ | Maximum daily demand
| $A_{\max}$ | Maximum daily order quantity
| Xt |  Vector of stock on hand at the start of day $t$, ordered by ascending age
| $X_{i,t}$ |  Element of Xt with $i$ days of remaining useful life at the start of day $t$
| $X_{t}$ |  Total stock on hand at the start of day $t$
| $I_{t}$ |  Total stock on hand and in transit at the start of day $t$
Scenario A | Ot |  Vector of stock in transit at the start of day $t$, ordered by ascending age
| $O_{i,t}$ |  Element of Ot that will arrive in $i$ periods at the start of day $t$
| $O_{t}$ |  Total stock in transit at the start of day $t$
| $W_{t}$ |  Number of units that expire at the end of day $t$
| $\mu$ | Mean of the gamma distribution for daily demand
| $\frac{\mu}{\sigma}$ | Coefficient of variation of the gamma distribution for daily demand
Scenario B | $H_{t}$ |  Number of units of a product issued to fill demand arising on day $t$
| $D_{t}^{u}$ |  Excess demand for product B willing to accept product A on day $t$
| $D_{t}^{z}$ |  Total demand for product A, including any substitution, on day $t$
| $\mu$ | Mean of the Poisson distribution for daily demand