Stochastic Optimization of Inventory at Large-scale Supply Chains

Zhaoyang Larry Jin, Mehdi Maasoumy, Yimin Liu, Zeshi Zheng, Zizhuo Ren
C3 AI
{[Link], [Link], [Link], [Link], [Link]}@[Link]
arXiv:2502.11213v1 [[Link]] 16 Feb 2025

Abstract ply chains exposed the inherent volatility of such networks.

With this volatility to be continued in the post-pandemic
Today’s global supply chains face growing challenges era, a solution is required to handle uncertainties in the
due to rapidly changing market conditions, increased net- global supply chains. While battling the uncertainties by in-
work complexity and inter-dependency, and dynamic uncer- creasing capital spent on raising the inventory level, the en-
tainties in supply, demand, and other factors. To combat terprises begin to pay more attention at the concept of lean
these challenges, organizations employ Material Require- manufacturing which includes maximizing the resilience of
ments Planning (MRP) software solutions to set inventory the inventory availability while minimizing the amount of
stock buffers – for raw materials, work-in-process goods, capital expenditure on inventory held.
and finished products – to help them meet customer service Current practice in today’s supply chains suffer from
levels. However, holding excess inventory further compli- numerous issues covering aspects including suppliers,
cates operations and can lock up millions of dollars of cap- transportation, manufacturing, distribution and customers
ital that could be otherwise deployed. Furthermore, most (Fig. 1). Perturbations in each of the aspects are hardly spe-
commercially available MRP solutions fall short in consid- cific only to itself. The network nature of the supply chains
ering uncertainties and do not result in optimal solution for makes any issue in each of these aspects easily propagate
modern enterprises. and amplify upstream or downstream. On one direction,
At C3 AI, we fundamentally reformulate the inventory dynamic demand and order changes would impact produc-
management problem as a constrained stochastic optimiza- tion schedules and inventory management at manufacturing
tion. We then propose a simulation-optimization framework sites or distribution centers. In addition, available inven-
that minimizes inventory and related costs while maintain- tory and production requirements would inform supplier or-
ing desired service levels. The framework’s goal is to find ders. The equipment reliability, product quality issues, and
the optimal reorder parameters that minimize costs sub- shipment uncertainties would impact part orders. On the
ject to a pre-defined service-level constraint and all other other direction, supplier uncertainties and delays would af-
real-world operational constraints. These optimal reorder fect production schedules, while order delays and quality
parameters can be fed back into an MRP system to drive issue would risk customer satisfaction. Furthermore, mate-
optimal order placement, or used to place optimal orders rial movement delays would strain the entire supply chain.
directly. This approach has proven successful in reducing The Lack of end-to-end operational visibility would lead to
inventory levels by 10-35 percent, resulting in hundreds of delays across supply chains, excess inventory levels, quality
millions of dollars of economic benefit for major enterprises issues, poor OTIF, or low customer satisfaction.
at a global scale.
Material Requirements Planning (MRP) is at the heart of
the process in today’s supply chain systems. MRP uses a
1. Introduction deterministic logic that calculates the time and quantity of
the orders that need to be placed, based on net material re-
Global supply chain faced enormous challenges due to quirement, re-order parameters, current inventory, demand
uncertainties introduced by the shocks in both supply side forecast and other inputs (Fig. 2). MRP computes material
and demand side during the COVID-19 pandemic. A sup- requirements based on demand forecast and schedules the
ply shock hit as the pandemic first emerged from main- arrivals as close as possible to the time they will be con-
land China and disrupted global supply. Over the ensu- sumed. Re-order parameters such as Safety Stock Value
ing months, global shutdowns and stay-at-home orders trig- (SSV) and Safety Time (ST) for the Safety Stock MRP sys-
gered a second shock, this time to the global demand. These tem, or Reorder Point (RoP) for the Reorder Point MRP
two waves of unprecedented disruptions to the global sup- system are used as inputs to MRP to address the uncertain-

1
 heuristics, the solution of the MRP in existing software are
not optimal.
Commercially available MRP solutions suffer from three
major drawbacks:
• The dynamic programming solution to the optimal
control formulation of MRP algorithms were perceived
to be too complex. Hence, commercial MRP systems
Figure 1: Common issues in manufacturing ranging from rely on grossly simplified approximate heuristics solu-
suppliers to customers tions to the problem which do not guarantee optimality
of the MRP solution.
• MRP relies on input parameters to account for uncer-
ties in both supply and demand by providing buffers in the
tainties associated with supply and demand functions
system.
that are set by end-users in an ad-hoc fashion, without
MRP receives as input, 1) the bill of materials (BOM)
rigorous analysis of uncertainty associated with histor-
which is a list of the materials, components and sub-
ical supply and demand.
assemblies required to make each product, 2) End Products
which is also called independent demand or level ”0” on • Approximate solutions to MRP do not learn from
the BOM, 3) Quantity that is required at a given time, 4) nor adapt to the uncertainties and the time-varying
Time at which the quantities are required to meet demand, and customer- or supplier-dependent operational con-
5) Inventory status records, which is a record of net mate- straints surrounding the manufacturing process.
rials available for use already in stock and materials on or-
der from suppliers, and 6) Planning data which includes all Such limitations along with disparate, siloed data across
the restraints and directions to produce such items as rout- multiple enterprise resource planning (ERP) systems result-
ing, labor and machine standards, quality and testing stan- ing in stale and limited insights for the planners have led to
dards, lot sizing techniques (i.e. fixed lot size, lot-for-lot, overly conservative inventory estimates due to hard to quan-
economic order quantity), scrap percentages, and other in- tify uncertainties.
puts. Furthermore, a lack of end-to-end visibility and a
Demand forecast reflects the independent demand which system-level optimal decision making makes today’s supply
is the demand for the end product such as a computer or a chains suffers form several issues. Delayed arrivals of pur-
bicycle. Dependent demand, on the other hand, is demand chase parts by supplier, causes production delay for down-
for component parts or sub-assemblies. For example, this stream components and delay in production of components
would be the microchips in the computer, the wheels on the results in production delay for top-level products. There-
bicycle. Quantities for dependent demand are derived from fore, purchase part delays lead to overall trapped inventory
independent demand using the Bill of Material (BOM). The and increased inventory holding cost. Most manufacturing
planning system needs to consider not only the quantities of processes suffer from poor forecasts. Demand forecasts end
each of the component parts needed, but also the lead times up being either under-forecasts which lead to stock-outs and
needed to produce and receive the dependent demand items. poor OTIF, or over-forecasts which lead to excess inventory
and hundreds of millions of dollars of capital locked up in
MRP can be expressed as an optimal control problem
unused inventory.
and the resulting optimal control problem can be solved
by dynamic programming (Zipkin, 2000). These dynamic
programming methods (e.g., dynamic lot-sizing) were per-
ceived to be too complex, A few approximate heuristics
(e.g., the Silver-Meal heuristic, 1973) were developed for
the problem. Such heuristics which are used in MRP sys-
tems to date, rely on gross simplifications. For example,
Silver–Meal heuristic is a forward-looking method that re- Figure 2: MRP follows a deterministic logic to place or-
quires determining the average cost per period as a function ders based on demand forecast, current day inventory, plan-
of the number of periods the current order is to span and ning calendar, order parameters and previously scheduled
stopping the computation when this function first increases. arrivals (ELT: expedited lead time, LT: lead time, PTF: plan-
Well-known methods to find order quantities are Dynamic ning time fence)
lot-sizing (Wagner and Whitin, 1958) introduced in 1958,
Silver–Meal heuristic (Silver, 1973) composed in 1973, and These challenges are presented as daunting problems in
Least-unit-cost heuristic. As a result of such simplifying supply chain management and we propose to solve them

2
with a stochastic inventory optimization (SIO) approach. In able additional relaxing assumptions have been made which
this article we will discuss a new formulation of inventory cannot be made on real examples in real world operations.
optimization problem as a constrained stochastic optimiza- Simulation-based method was brought into the picture be-
tion problem in Section 3. We discuss the data integra- cause of its capability of taking into account specific busi-
tion and model deployment platform that enabled solving ness restrictions, control variables and uncertainties of the
this problem at large scale in Section 4. We present results problem (Köchel and Nieländer, 2005). However, a pure
showing inventory reductions across multiple systems at a simulation-based approach won’t lead to identification of
global scale in Section 5. best parameters to operate the system with, which could
lead to optimal performance. As a result, a simulation-
2. Related Work optimization approach has has been adopted recently to ad-
equately model the complex system-level interactions and
Solutions the MRP algorithm provides a heuristic ap- constraints and solve for optimal re-order parameters which
proach for computing a sub-optimal but feasible solution leads to optimized inventory level and cost (Fu, 2002).
for optimization problem, given SSV, ST, demand forecast,
Olafsson and Jumi Kim (2002) provided an overview of
and other ordering parameters. To account for such limita-
simulation optimization methods which were later adopted
tions of MRP, new variations of MRP has been introduced.
to solve the supply chain optimization problems. Given the
For example, Ptak and Smith (2011) introduced a new type
involvement of multiple decision makers in the whole sup-
of MRP called demand driven MRP (DDMRP). DDMRP
ply chain network. Swaminathan et al. (1998) proposed
leverages knowledge from theory of constraints, traditional
a multi-agent approach to represent critical supply chain
MRP & distribution resource planning (DRP), Six Sigma
components (e.g., retails, manufacturers) and simulate the
and lean manufacturing. It leverages MRP for planning, and
supply chain based on their decisions. Mele et al. (2006)
kan-ban techniques for execution across multi-echelon sup-
further extended agent-based simulation-optimization ap-
ply chains, which means that it has the strengths of both but
proach by incorporating genetic algorithm. Joines et al.
also the weaknesses of both and it remains a niche solution
(2002) leveraged a multi-objective GA methodology with
(Velasco Acosta et al., 2020; Smith and Smith, 2014).
an existing supply chain simulator to optimize the system
Since the multi-commodity distribution system design
parameters. Besides directly optimize the inventory which
was first introduced by the work of Geoffrion and Graves
may be limited by operation constraints, a simulation-
(2010), multiple optimization-based approaches have been
optimization framework was proposed by Jung et al. (2008)
proposed for the design of supply chain networks (Aikens,
to optimize the safety stock levels with consideration of pro-
1985; Geoffrion and Powers, 1995; Vidal and Goetschal-
duction capacity.
ckx, 1997). Most of these research, however, made the
The majority of methods mentioned above did not prior-
assumption that the supply chain are deterministic without
itize treatment of the stochastic nature of inventory above
considering the uncertainties from different sources such as
the optimization. For the ones that formulate the issue as
demands and resource capacity. Supply chain disruption
stochastic programming problem, they lack the details in
due to the ignorance of uncertain operating conditions will
the modeling to properly represent the dynamics of inven-
cause huge economic impact. A stochastic programming
tory in supply chain. Those drawbacks make it difficult to
model was proposed by Santoso et al. (2005) to solve the
productionize those solutions at a large scale.
supply chain network design problem.
Since then, many attempts have been made to solve
stochastic nature of supply chain network by leveraging 3. Methodology
mathematical methods. formulated A two-stage stochas-
3.1. Problem formulation
tic program targeting at an optimal operation plan was for-
mulated by Ierapetritou and Pistikopoulos (1994) and a We cast the problem at hand as a stochastic constrained
decomposition-based optimization approach was proposed optimization problem and solve for the optimal set of re-
to solve the program. Further improvement on the approach order parameters that minimize inventory levels subject to
was made by incorporating a MILP planning model (Gupta a given service level constraint. The objective function
and Maranas, 2000) . A non-linear mixed integer formula- in the optimization problem includes costs associated with
tion was proposed by Gupta and Maranas (2000) with non- holding inventory, ordering and transportation. This cost
convex objective and constraint functions. is optimized with respect to operational constraints. These
The main issue with the methods above is that they have operational constraints include planning calendar, planning
oversimplified the supply chain in order to model it as a time fence, (expedited) lead time, rounding values, mini-
mathematical optimization, which has resulted in the omis- mum order quantity, customer-dependent operational con-
sion of critical real-world operational complexities (Agar- straints, complex material pricing (e.g., volume-based pric-
wal, 2019). Furthermore, it order to make the problem solv- ing). Another important constraint is often the service level

3
constraint which is defined as expected probability of not γp (·) offers more flexibility than aggregating with the ex-
hitting a stock-out during the next replenishment cycle. pectation (i.e., E[·]), although the latter is more often used
The stochasticities, or uncertainties, in the constrained as the as the aggregation function for stochastic optimiza-
optimization problem come partly from the fact that there tion problems.
are often unexpected delays and shortages on the arriving As mentioned in Section 1, in all large scale operational
orders. In addition, the miscellaneous and blocked material settings where we intend to deploy our C3 SIO solution, the
movements that cannot be anticipated ahead of time are also reordering of materials is done via legacy MRP systems.
categorized as uncertainties. Finally, the optimization result In order to respect that constraint and to ease the customer
is required to be frequently adjusted based on the latest es- adaptation, we optimize for the reorder parameters, which
timation of the future demand, which is subject to change are the inputs of the MRP system. Therefore, the decision
from one day to the next and is often different from the ac- variables in this formulation are denoted by R:H , which rep-
tual demand, leading to another source of uncertainty. resents the reorder parameters of the MRP system up to step
To address the constraints and the uncertainties men- t + H. To be specific, R:H is expressed as Eq. 2 below
tioned above, we formulate the problem in a model pre-
dictive control (MPC) setting, which is a well-established R:H = [rt|t , rt+1|t , . . . , rt+H−1|t ] ∈ RH×2 (2)
framework in the process control discipline Garcia et al. where the variable r at a given timestamp is a single pair
(1989); Mayne et al. (2000). At each timestamp t (also re- of reorder parameters for the MRP system. Note that our
ferred to as the sampling time), our problem can be formu- formulation would work for MRP systems with a variety of
lated as the following optimization problem replenishment strategies. However, throughout the paper,
"H−1 # we limit the scope of our discussion within the Safety Stock
X
min γp t+i t+i
xt+i|t · πh + at+i|t · πo,t , MRP system. Therefore, the reorder parameters here in-
R:H clude the safety stock value (SSV) and the safety time (ST),
i=0
s.t. :H
ψ(X ) − SLmin ≥ 0, and the variable rt+i|t can be represented by Eq. 3 below

xt+i+1|t = xt+i|t + at+i|t − ct,t+i

DF +Θ rt+i|t = [SSVt+i|t , STt+i|t ] ∈ R2 (3)
∀i ∈ [0, . . . , H − 1], (1)
The variable X:H in the first constraint of Eq. 1 repre-
i+1:H t+i|t t+i t+i|t
A = MRP(x , CDF ,r , Ai:H ) sents the vector of stochastic inventory level in the predic-
∀i ∈ [0, . . . , H − 1], tion horizon H, and is clearly defined in Eq. 4 below:

r t+i+1|t
= rt+i|t ∀i ∈ [0, . . . , H − 1], X:H = [xt|t , xt+1|t , . . . , xt+H−1|t ] ∈ RH (4)
:H
R ≥ 0, where the bold lower-case notation x ∈ R1 with super-
scripts represents the stochastic inventory level for a spe-
where the variable H denotes the length of the prediction
cific date, predicted at time t. The variable SLmin represents
horizon, which limits the length of the inventory simulation
the minimum target service level. The function ψ(·) calcu-
we perform at sampling time t. The superscript represents
lates the service level within the prediction horizon. It can
the timestamp, and the t + i|t notation in the superscript
be expressed as Eq. 5 below
depicts the value on timestamp t + i predicted at the times-
tamp t. The variable x ∈ R1 denotes the stochastic in-
"P #
H t+i|t
:H i=1 δ(x )
ventory level, and a ∈ R1 depicts the stochastic expected ψ(X ) = γpSL ∈ R1 , (5)
H
order arrival quantity. Note that the expected order arrival
quantity can be further split into the expected standard ar- where the notation δ(·) is a level set function, which is de-
rival (SA) and the expected expedited arrival (EA), which fined as
will be discussed in Section 3.2.1. The variable πh is the
(
1, x ≥ 0
inventory holding cost, and πo,t represents both ordering δ(x) = (6)
0, otherwise
and transportation costs. Note that objective function could
take different forms depending on the business setting our Note here the service level is computed as the percentage
customer has. of days on which the inventory level is not negative. There
The notation γp (·) in the objective function of Eq. 1 is are other definitions of the service level, such as fill rate
the aggregation function that takes the specific percentile of (details in Chu et al. (2015)), which is also supported by
a given uncertainty distribution, where the variable p repre- this formulation. Also note that the aggregation functions
sents the percentile of the uncertainty distribution, (e.g., if γp (in Eq. 1) and γpSL (in Eq. 5) are identical, but p and pSL
p = 50%, γp returns the median of the uncertainty distribu- could take different values, depending on the risk profile of
tion). We want to point out that aggregating with function the customer. We use different notation here for clarity.

4
 The second constraint in Eq. 1 is a simple inventory up- results when we optimize the decision variables is the ex-
t,t+i
dating logic, where cDF ∈ R1 denotes the demand fore- istence of the uncertainties (represented by Θ). Should the
cast for timestamp t + i made at timestamp t, the variable uncertainties in the formulation be completely removed, the
Θ denotes a collection of stochastic variables representing decision variables will be all 0 since legacy MRP systems
various sources of uncertainties. It includes the movement are capable of handling the deterministic inventory planning
uncertainty UM M , the supplier quantity uncertainty USQ , without any buffer specified.
the supplier time uncertainty UST , and the demand forecast
uncertainty UDF , all of which will be discussed in detail in 3.2. Algorithm overview
Section 3.2.2.
Once the objective function and the constraints are con-
The function notation MRP(·) in the third constraint
figured (e.g., the hyper-parameters including γp in Eq. 5 de-
in Eq. 1 represents the constraints imposed by the legacy
termined), Eq. 1 can be solved for at each sampling time.
MRP system. It takes the inventory level up to the current
However, in practice, it is difficult for the customer to pro-
timestamp xt+i|t , the demand forecasts at the current times-
t+i vide those hyper-parameters beforehand. Therefore, as part
tamp CDF , the reorder parameters at the current timestamp
t+i|t of the solution, we need to help the customer to configure
r , and the expected order arrivals starting at the current
the hyper-parameters that fit best for their business inter-
timestamp Ai:H , then it updates (overwrites) the expected
est. In addition, MRP systems can be highly customized
order arrivals starting at the next timestamp to the end of the
t+i depending on the customer. There are pieces of logic such
predictive horizon (i.e., Ai+1:H ). The variable CDF ∈ RH
as order cancellation and order merging (i.e., flush window)
is the demand forecast made at timestamp t + i, and it can
that are difficult if possible at all to formulate as simple
be expressed as
linear constraints for the optimization problem. This fact
t+i
CDF = [ct+i,t+i t+i,t+i+H−1
, . . . , cDF ] ∈ RH . (7) prevents us from directly feeding the optimization problem
DF
formulated in Eq. 1 in a traditional mixed-integer linear pro-
The notation Aτ :H represents the stochastic expected or- gramming (MILP) solver.
der arrivals (include both standard arrivals and expedited To handle the challenges mentioned above, the C3 AI In-
arrivals) starting at timestamp t + τ (∀τ ∈ [0, . . . , H − 1]) ventory Optimization application introduces an AI Stochas-
to the end of the predictive horizon (i.e., timestamp t + H). tic Inventory Optimization Algorithm (“the algorithm”).
More specifically, The algorithm learns a material plant’s historical uncertain-
ties in material movements - supply, demand, and other -
Aτ :H = [at+τ |t , . . . , at+H−1|t ] ∈ RH−τ . (8) over a historical time period (“training period”) to select
the best hyper-parameters that are feasible to the customer.
Note that the initial value for Ai:H in Eq. 1 while i = 0
Then, the algorithm perturbs these learned material move-
(noted as Ã:H ) comes from the the solution of Eq. 1 on the
ment uncertainties into the material planning simulations
previous sampling time t̃, which can be explicitly written as
in the future (i.e., “validation period”, and “live produc-
Ã:H = [at|t̃ , . . . , at+H−1|t̃ ] ∈ RH . (9) tion”). Based on those simulations and uncertainty realiza-
tions in the future, and with the help of the selected hyper-
More details on this MRP constraint can be found in Algo- parameters, the algorithm recommends optimized reorder
rithm 1. parameters (i.e., SSV and ST for the Safety Stock MRP sys-
The fourth constraint in Eq. 1 represents the limitation tem) that minimize inventory costs while meeting the target
of the current MRP system, where only a single value can service level constraint.
be specified for each of SSV and ST through out the entire Fig. 3 demonstrates the general workflow of the algo-
prediction horizon. Note that the reorder parameters are not rithm. The algorithm has two phases: a training (offline)
required to be identical when the sampling time t evolves. phase and a validation or live production (online) phase.
For example, rt+i|t is not necessarily equal to rt+i|t+j for The training phase happens on a period of time further in the
j > 0 and i > j. The fifth constraint in Eq. 1 indicates that history while the validation phase is on a more recent time
the reorder parameters are always non-negative. period. The term “simulation period” in this paper refers to
The formulation presented in Eq. 1 will be solved itera- the timestamps (i.e., days) either in the training or in the val-
tively for different sampling time t. The reorder parameters idation phase depending on the context. Note that the live
recommended at the sampling time t will be preserved till production is a special case of the validation phase where
it is updated at the next sampling time t + j. Note that the the optimization is performed at “today” (the last date we
sampling time increment j (also referred to as the optimiza- have the data available). Therefore, the reorder parameters
tion frequency) is not required to be 1. More details on the are only optimized for one timestamp in the live production
optimization frequency will be discussed in Section 3.2.6. phase.
Also note that the reason we would expect to see non-zero We start the discussion of Fig. 3 from the validation

5
 Figure 3: Workflow of the C3 AI Stochastic Inventory Optimization algorithm

phase, where the optimization problem formulated in Eq. 1 proach configurable by the user) of the realizations during
is solved repeatedly at each sampling time t to generate the the validation period.
final simulated inventory for the validation period. The in- The live production (or, operation) is essentially identical
ventory simulation of the validation phase has three key to the validation phase. The only distinction is that there is
components (steps): (1) uncertainty distribution step, (2) only one timestamp (day) for reorder parameter recommen-
reorder parameters estimation step, and (3) order simula- dation in the live production. Also, since we will not have
tion step. During the uncertainty distribution step, all the the actual data (e.g., actual inventory, actual arrivals, actual
uncertainties of the model inputs are collected within the consumption, etc.) for the live production, there is no need
uncertainty sampling window (USW). Those uncertainties to compare the simulated results with the actual inventory.
are saved for the sampling procedure that happens in the The purpose of the training phase is to determine the
reorder parameters estimation step. In the reorder param- hyper-parameters, which helps the inferring of the reorder
eters estimation step, the MRP simulations are performed parameters at each sampling time t in the validation phase.
repeatedly starting at each sampling time t to infer the best The training phase (shown in blue in Fig. 3) is essentially a
reorder parameters. The inference of the re-order parame- grid search process repeating the validation procedure mul-
ters is accomplished by the “k-iteration” process (with the tiple times in the training period, each time with a differ-
hyper-parameters determined by the training phase), which ent combination of hyper-parameters. Then the best hyper-
will be discussed in Section 3.2.5. The reorder parameters parameter is determined based on the inventory savings and
will be retained till they are overwritten by the inference service level achieved by each combination. Note that the
made at the next sampling time. The algorithm is designed grid search approach is chosen among all optimization op-
with this overwriting process to effectively incorporate the tions to ensure the stability of the algorithm, which will be
new information that becomes available at each sampling discussed in Section 3.2.3.
time t. To summarize, the deployment of the algorithm includes
three key steps:
Subsequently, in the order simulation step, the future or-
der placements will be generated based on a series of re- 1. Train the model: the model learns from material move-
order parameters. Finally, the algorithm incorporates the ments and associated uncertainties over the historical
uncertainty distributions again to generate multiple realiza- training period and optimizes the hyper-parameters.
tions of the simulated inventory. The inventory cost of each
realizations can then be computed accordingly. Inventory 2. Validate the model: the trained model (with optimized
savings can be computed by comparing the actual cost of hyper-parameters) is validated in a production-like op-
the inventory with the median (or any other aggregation ap- eration over a period of recent history to evaluate real

6
 performance (i.e., inventory cost and service level). simulates inventory for the horizon x̂t+i , i = 0, ..., H. No-
tice that hat accent is used to denote inventory position
3. Operate the model (live production): The model is run within MRP. When ELT <= i < LT , where ELT and
in production, simulating uncertainties into the future LT are expedited lead time and lead time, respectively, the
to provide recommendations on the optimal reorder pa- algorithm places new expedited orders if projected inven-
rameters. The model’s performance is monitored; it tory x̃t is about to go below 0. When LT <= i < H, the
may be automatically or manually retrained as needed. algorithm places new standard orders if projected inventory
is about to go below safety stock value SSV . The quantity
In the following parts of this section, we will cover top- of new standard orders are modified to comply with mini-
ics including the MRP simulation, the modeling of uncer- mum order quantity M O and rounding value RV .
tainties, optimization of the hyper-parameters, the forward-
looking simulation, the “k-iteration” algorithm, and other
Algorithm 1 Simplified Safety Stock MRP
algorithm treatments.
Input: starting inventory xt , demand forecast CDF t
,
t+i
expected standard order arrivals SA , expected expedited
3.2.1 Inventory simulation with MRP order arrivals EAt+i , where i = 0...H, and H is the
planning time horizon, lead time LT , expedited lead time
As introduced previously, the MRP system sits at the core
ELT , safety stock value SSV , minimum order quantity
of C3 AI SIO Algorithm. For any given sampling time t,
M O, rounding value RV
the MRP takes the reorder parameters along with other in-
t+i Output: updated expected arrivals of standard orders
puts including demand forecast (CDF ), lead time (LT), ex-
SAt+i and expedited orders EAt+i , i = 0...H
pedited lead time (ELT), planning calendar (PC), item vol-
ume pricing (IVP), etc. within a horizon (corresponds to
the predictive horizon H in Eq. 1, then outputs the order 1: Apply order cancellation on expected standard order ar-
placement in that horizon. rivals (SA) based on Algorithm 2
2: x̂t = xt
As is explained in Section 3.1, the optimization frame-
3:
work introduced by C3 AI SIO could support MRP systems
with a variety of replenishment strategies. However, we 4: for i = 0, . . . , ELT − 1 do
limit the scope of our discussion in this paper to the Safety 5: x̂t+i+1 = x̂t+i + SAt+i + EAt+i − ct,t+i
DF
6: end for
Stock MRP where the reorder parameters are safety stock
7:
value (SSV) and safety time (ST).
8: for i = ELT, . . . , LT − 1 do
9: x̃t+i+1 = x̂t+i + SAt+i + EAt+i − ct,t+i
DF
Safety Stock MRP 10: if x̃t+i+1 < 0 then
11: EAt+i += −x̃t+i+1
The Safety Stock MRP system places orders for a win- 12: end if
dow in the future based on the current inventory position, 13: x̂t+i+1 = x̂t+i + SAt+i + EAt+i − ct,t+i
DF
existing orders and demand forecast. The MRP projects the 14: end for
inventory position for every day in this future window. Or- 15:
ders are placed to arrive when the projected inventory drops 16: for i = LT, . . . , H − 1 do
below a predefined safety stock value, subject to reorder 17: x̃t+i+1 = x̂t+i + SAt+i + EAt+i − ct,t+i
DF
constraints such as lead time, factory calendar, planning cal- 18: if x̃t+i+1 < SSV then
endar, minimum order quantity and rounding values, etc. 19: ∆x = SSV − x̃t+i+1
Algorithm 1 provides a simplified view of the Safety 20: SAt+i += max(⌈ ∆x−M RV
O
⌉, 0) × RV + M O
Stock MRP system, where the key pieces of logic are de- 21: end if
scribed in details. However, some auxiliary pieces logic in- 22: x̂t+i+1 = x̂t+i + SAt+i + EAt+i − ct,t+i
DF
cluding item volume pricing (IVP), lot size in days (LS), 23: end for
planning calendar (PC), order flushing (or, merging) are in- 24:
tentionally left out in this algorithm to improve the readabil- 25: return SAt+i , EAt+i i = 0...H
ity of the algorithm. Order cancellation logic, on the other
hand, is explained in the next sub-section (Algorithm 2).
Input data to Algorithm 1 include starting inventory xt ,
t
demand forecast CDF , the expected arrivals of existing stan- Order cancellation
t+i
dard orders SA and existing expedited orders EAt+i , Orders sent to vendors can be cancelled based on a can-
where i = 0, ..., H, and H is the horizon. The algorithm cellation window (e.g., 30 days prior to arrival). Usually

7
the orders can be cancelled when the orders’ expected ar- Depending on the accuracy of the demand forecast
rival date is later than the planning-time-fence (P T F ) day, model, the predicted demand can be quite different from the
counting from the current day. Within the MRP module, we actual demand/consumption. Demand forecast uncertainty
are enabled to cancel orders when the P T F is greater than models the distribution of the difference between demand
the lead time (LT ). However, in the case that P T F be- forecast and actual consumption.
ing smaller than LT , MRP does not cancel the orders that There are two types of material movements, miscella-
will be arriving in between [P T F, LT ] if more quantity is neous movements, and blocked movements. These are typi-
ordered than needed. cally not considered in traditional MRPs and are not part of
We developed a pre-processing algorithm (Algorithm 2) the amount represented in the demand forecast. Miscella-
for MRP to handle the cancellation of the excessive orders. neous movements happen due to movement between facili-
The algorithm modifies the orders that were placed pre- ties or customer returns. Blocked movements happen when
viously but are expected to arrive in between [P T F, LT ]. parts arrived but cannot be consumed (e.g., due to a quality
Rather than iterating incrementally from P T F to LT , it it- issue). As these material movements are generally unpre-
erates in the backward order from LT to P T F to make sure dictable, they are treated as sources of uncertainty. Material
the cancelled order does not impact future inventory levels. movement uncertainty models the distribution of the sum
Within each iterated date, it computes the surplus of inven- of miscellaneous movements and blocked movements for a
tory, which is denoted by x̃t+k in the algorithm. If the sur- given SKU.
plus is greater than zero, the expected arrivals for that date Supplier quantity uncertainty models the distribution of
would be reduced to offset the surplus. the difference between the planned quantity in the purchase
orders and the actual quantity received for a given SKU.
Algorithm 2 The order cancellation logic Supplier time uncertainty models the distribution of the
Input: inventory xt of the day t, expected standard time difference of the planned delivery date and the actual
order arrivals SAt+i , expedited order arrivals EAt+i and delivery date of orders for a given SKU.
demand forecast ct,t+i
DF where i = 0...H, and H is the plan-
In reality, different suppliers can have different behav-
ning time horizon for a single round of MRP simulation, iors. One supplier may be more likely to delay and reduce
lead time LT , planning time fence P T F , and safety stock quantity than another supplier. However, we usually don’t
value SSV have enough data points to accurately model the behavior of
Output: modified expected order arrivals SAt+i where each supplier. Therefore, we average out different behaviors
i = 0...H of all suppliers for a given SKU into one distribution.
1: Set k = LT We estimate the distribution of each source of uncer-
2: while k > P T F do tainty using an empirical distribution. Each empirical distri-
3:
Pi=k Pi=k
x̃t+k = xt + i=0 SAt+i + i=0 EAt+i − bution is represented by a list of historical data points. We
Pi=k t,t+i assume that all historical data points are equally likely, then
i=0 cDF − SSV we can sample the empirical distribution from the list of
4: if x̃t+k ≤ 0 then
historical data points. The list of historical data points is re-
5: k =k−1
trieved from relevant quantities within a certain uncertainty
6: continue to next iteration
sampling window (USW).
7: end if
The USW for training is the entire training period. The
8: if x̃t+k > SAt+k then
USW during validation is a growing window with a fixed
9: Set SAt+k := 0
start date, where the initial length of the USW L0USW is given
10: else
by
11: Set SAt+k := SAt+k − x̃t+k
12: end if LUSW
0 = max(LUSW min , LT + b
USW
), (10)
13: k =k−1 where LT represents lead time, bUSW represents the buffer
14: end while (or, the offset) of the USW length and LUSW min denotes the
15: return SA minimum initial length of USW.
The empirical distribution for material movement uncer-
tainty is represented by UMM , which is a list of data points
3.2.2 Modeling of the uncertainties in the uncertainty sampling window, where each data point
is the sum of miscellaneous movement and blocked move-
In our algorithm, we consider four major sources of un- ment at that time step. To be more specific,
certainties, including demand forecast uncertainty, material
movement uncertainty, supplier quantity uncertainty and USW
supplier time uncertainty. UMM = [..., cτMIM + cτBM , ...] ∈ RL , ∀τ ∈ USW, (11)

8
where LUSW is the length of the uncertainty sampling win- representation of the demand forecast, where LT and ST
dow, cτMIM ∈ R1 and cτBM ∈ R1 denotes the value of stand for lead time and safety time, respectively. The 1-
consumption due to miscellaneous movement and blocked dimensional representation can be expressed as
movement at time τ , respectively. The variables cτMIM and
cτBM are positive when materials move in and negative when USW

materials move out. CDF, 1D = [..., cτDF−LT−ST,τ , ...] ∈ RL , ∀τ ∈ USW. (15)

We use USQ to denote the list of historical data points that
Here we use demand forecast made on lead time plus safety
represents the empirical distribution of the supplier quan-
time days before as we believe this most accurately reflects
tity uncertainty. Each data point in USQ is the difference
the actual demand forecast uncertainty. Note that CDF, 1D
between planned quantity and actual received quantity of
effectively comes from the off-diagonal terms of the matrix
an order that is planned to be delivered in the uncertainty
CDF, 2D .
sampling window. It can be expressed as
In practice, demand forecast oftentimes comes in on a
weekly basis or monthly basis. To convert the demand fore-
USQ =[..., min(soA − soP , 0), ...] ∈ RL
USW
, cast into a daily one, so that it is comparable to the ac-
(12) tual consumption, we apply a convolution step on CDF, 1D
∀o ∈ Orders planned to arrive in USW, that acts as a moving average smoother. Specifically, when
demand forecast comes in once everyday d, we convolve
where soA ∈ R1 and soP ∈ R1 denote the actual and planned 1 1
supplied quantity from order o, respectively. To make sure CDF, 1D with a d-dimensional constant vector a = [ , ..., ]
d d
our algorithm is conservative in nature, when the actual sup- to obtain the smoothed 1-dimensional demand forecast
plied quantity SAo is larger than the planned supplied quan- ∗
CDF, 1D .
tity SPo , we do not take this favorable condition into con- i+d/2
i,∗ 1 X j
sideration when quantifying the uncertainty. Therefore, we CDF, 1D = CDF, 1D , (16)
will replace the positive number with 0 in the empirical dis- d
j=i−d/2
tribution. j
We use UST to denote the list of historical data points that where CDF, 1D represents the j-th element in CDF, 1D de-
i,∗
represents the empirical distribution of the supplier time un- fined in Eq. 15, and CDF, 1D denotes the i-th element in the
∗
certainty. Each data point in UST is the difference (in num- smoothed 1-dimensional demand forecast CDF, 1D .
ber of days) between planned delivery date and actual de- Let
livery date of an order that is planned to be delivered in the USW

uncertainty sampling window. CA = [..., cτA , ...] ∈ RL , ∀τ ∈ USW (17)

be the list of actual consumption in the uncertainty sampling

USW window. We compute the demand forecast uncertainty with
UST =[..., max(toA − toP , 0), ...] ∈ RL ,
(13) USW
∀o ∈ Orders planned to delivered in USW, ∗
UDF = CDF, 1D − CA ∈ R
L
. (18)

where toA ∈ R1 and toP ∈ R1 denote the actual and planned Note that in practice, we would encounter outliers in
delivery time (in days) from order o, respectively. When the CDF , so pre-processing steps (discussed in Section 3.2.6)
actual delivery date toA is earlier than the planned delivery are necessary before computing UDF or consuming CDF in
date toP , again, this is a favorable condition that we do not the algorithm.
take into consideration. Therefore, we will instead put a
value of 0 into the empirical distribution. 3.2.3 Optimization of the hyper-parameters
To construct the demand forecast uncertainty UDF , we
need to compress the 2-dimensional demand forecast ex- As mentioned earlier in the algorithm overview, we will
pressed in Eq. 14 help the customers to decide the risk profile that fits best
to their business interest. The risk profile can be quanti-
USW
τ ×H fied by hyper-parameters including pSL for γpSL in Eq. 5.
CDF, 2D = [..., CDF , ...] ∈ RL , ∀τ ∈ USW. (14)
This hyper-parameter is also referred to as service level per-
into a 1-dimensional array before subtracting it by the ac- centile, or SLP. As will be discussed later, in the context
tual consumption. To be specific, following Eq. 7, let of Safety Stock MRP system, the SLP directly affects the
to ,tf
cDF ∈ R1 be the demand forecast (planned consumption) choice of SSV in the C3 SIO algorithm, but it does not
for time tf made at time to . For each time τ in the uncer- have any impact on the choice of ST. To help with quan-
tainty sampling window, extract the demand forecast value tifying customer’s risk profile on the time delay and recom-
made for time τ on time τ − LT − ST to get a 1-dimensional mending ST, we introduce another hyper-parameter, safety

9
time percentile (or, STP). Although those two variables are STP, on the other hand, describes the percentile in the
called “hyper-parameters”, they are analogous to the “pa- supplier time uncertainty distribution. The strategy used in
rameters” in a regular machine learning model, which is op- the Safety Stock MRP system to de-risk the supplier de-
timized during the training process and directly used in the lay is to place the purchase order earlier. The safety time
model for the test period. They are referred to as “hyper- is the variable that determines how early the order should
parameters” in C3 SIO algorithm to distinguish from other be placed compared to the default order placement date in
model parameters that are not optimized during the training the MRP output. Therefore, the choice of the most suitable
process. safety time depends on the supplier time uncertainty profile,
To distinguish SLP from the service level (SL), we note which could vary over the simulation period. A fixed choice
that SL describes the performance of a single simulated in- of STP could map to different safety time while the supplier
ventory curve. However, SLP describes the performance of time uncertainty profile changes over time.
multiple simulated inventory curves, and it is not a replace- Fig. 5 demonstrates what distribution are we looking at
ment but a complement of SL. SLP describes the percent- while deciding the safety time percentile. The plot on the
age of inventory realizations that meets the minimum target left shows the delay of each arriving order in the USW (i.e.,
service level (i.e., SLmin in Eq. 5). For example, if the min- the delays are represented by the differences between the
imum target service level is 96% and the total number of green peaks and the orange peaks). Those delays are trans-
uncertainty realization of simulated inventory is 100, which lated into an empirical distribution of supplier time uncer-
means there will be 100 simulated inventory curves. A ser- tainty as is shown in the right plot. Safety time will then
vice level percentile of 50% means that 50 out of the total be determined based on the STP and the empirical distri-
100 realizations will have a service level above 96%. bution. For example, a STP of 50% would translate into
a safety time of 3 days with the empirical distribution pro-
SLP will help the algorithm to recommend SSV at each
vided in Fig. 5, since the 50th percentile of the empirical
sampling time during the “k-iteration” process, which will
distribution is 3 days, and a STP of 100% would translate
be described in detail in Section 3.2.5. SLP decides how
into safety time of 6 days. Similar to the choice of SLP, a
conservative or aggressive the SSV choice will be. Fig. 4
high STP represents a conservative strategy since it would
provides a visual illustration on the discrepancies between
make MRP to place all the orders so early that even the
a high SLP and a low SLP with the grey areas represent the
longest delay would be prevented. This STP would result
uncertainty range of the simulated inventory. Note that we
in a high service level, however, the inventory level will be
are using a minimum target service level of 100% while dis-
piling up due to the early placement of orders and the in-
cussing about SLP in Fig. 4. A high SLP (e.g., 100% as is
ventory cost will be high. On the other side, a low STP is
shown on the left of the figure) puts the entire uncertainty
a aggressive strategy since it will result in a smaller safety
range of the simulated inventory above 0 (i.e., for 100 out
time. Therefore, the inventory is prone to stock-out, but the
of 100 uncertainty realizations, each of the realization has a
inventory cost will be lower.
service level of 100%), while a low SLP (e.g., 50% on the
right side of the figure) only puts the median curve above 0.
The low SLP is a more aggressive strategy, which means the
inventory will have a higher chance of going below the min-
imum target service level (high risk for stock-out), but the
inventory cost will be much lower (high return for inventory
savings). For now, we pick the SLP in the training based
on the inventory cost of the median realizations (curves in
black, which means we set p = 50% for γp in Eq. 1), but it
is configurable by the customer. Figure 5: Safety time percentile is the percentile on the sup-
plier time uncertainty

In the training period color-coded in blue in Fig. 3, the

best combination of SLP and STP is selected by the grid
search approach (illustrated in Fig. 6) so that the median of
the simulated inventory realizations would have the lowest
inventory cost while its service level is above the target ser-
vice level. Note that we apply the objective function and
the constraint on the median of the realizations as a design
Figure 4: Comparison of high versus low service level per- choice, which is configurable to fit the user’s preferences.
centiles Fig. 7 provides a holistic picture of how the hyper-

10
 Figure 6: Grid search of hyper-parameters
Figure 7: Selection and consumption of the hyper-
parameters
parameters are selected and consumed in the algorithm.
During the training period, the median realization of the
simulations for all the combination of SLP and STP are gen- 1. The demand forecast made on t for the forward-
erated, as is shown in the top part of the figure. Among all t
looking window CDF ∈ RH ,
the combinations, one combination with SLP equal to a and
STP equal to b gives the lowest inventory cost while satis- 2. The existing orders scheduled to arrive within a seed-
fies the constraint on the service level. Therefore, this com- ing window
bination of (a, b) is selected as the best hyper-parameters.
s −1
This hyper-parameter pair is then used in the validation pe- Atseed
s
= [atseed , . . . , at+t
seed ] ∈ Rts , (19)
riod. In the validation period, the SLP and STP are trans-
lated to SSV and ST recommendations at each timestamp where ts is the length of the seeding window, and
of the simulation. Then the inventory is simulated based aτseed , ∀τ ∈ [t, . . . , t + ts − 1] represents the seeded
on those recommendations to generate simulated inventory arrival on timestamp τ ,
curve shown as the light blue curve in the bottom figure.
3. The candidate SSV and ST,
3.2.4 Forward-looking simulations 4. The initial inventory position at sampling time xt0 (see
Section. 3.2.6 for more details).
After selecting the hyper-parameters, we enter the reorder
parameter estimation step. To find the optimal reorder pa- Note that the size of the seeding window ts varies depend-
rameters at a given sampling time t, we simulate inventory ing on the customers business logic, but oftentimes it is set
within a forward-looking window from t to H, where H to lead time (LT ).
is the horizon. This simulation is done with the MRP sys- Given the new orders scheduled by MRP, we generate a
tem and needs to be repeated multiple times, incorporating total of Nr realizations of inventory curves. The inventory
different samples from the uncertainty distributions. We re- at time τ + 1 for a realization r is given by
fer to this process as the forward-looking simulation, which
includes two steps:
xτ +1 (r) = xτ (r) − cτA (r) − cτMM (r) + sτs (r) ∈ R1
1. Run a MRP simulation within the forward-looking (20)
window; ∀τ ∈ [t, . . . , t + H − 1],

2. Generate multiple uncertainty realizations of the simu- where xτ (r) is the inventory at time τ for the r-th realiza-
lated inventory by perturbing the demand forecast, the tion, cτA (r) is the sampled consumption at time τ , cτMM (r) is
delay and quantity shortage of the arriving orders. the sampled material movement at time τ , sτs (r) is sampled
purchase order arrivals at τ . Note that the inventory level x
The main inputs to the MRP simulation include: here is no longer stochastic (not bold) with specific samples.

11
 Specifically, the sampled consumption at time τ is given Algorithm 3 The k-iteration algorithm (for the Safety Stock
by MRP)
cτA (r) = ct,τ 1
DF − uDF (r) ∈ R , (21) 1: Set ST to the value corresponding to the STP (which is
where ct,τ is the demand forecast made on time t for time invariant of the k-th iteration)
DF
τ , uDF (r) is the r-th sample from UDF (Eq. 18). 2: Set SSV = 0, h > 0
The sampled material movement at time τ is given by 3: while k < max iteration AND h > 0 do
4: Get Nr samples from each of the empirical uncer-
cτMM (r) = uMM (r) ∈ R1 , (22) tainty distribution (i.e., UST , USQ , UM M , UDF )
5: Apply the forward-looking simulation to get
where uMM (r) is the r-th sample from UMM (Eq. 11). X :H (r), ∀r ∈ [1, . . . , Nr ]
o
Let o1 , o2 , ..., on be the set of orders MRP placed, tej 6: Calculate the inventory cost and the SL associated
oj
and se , ∀j ∈ [1, ..., n] be the expected arrival time and ex- with each of the inventory realization r
pected quantity for order oj . The r-th sampled arrival time 7: if Eq. 27 does not hold then
o
for order oj , denoted by tsj (r) is given by 8: Compute the inventory deficit h (Fig. 8)
9: else
tosj (r) = toej + uST (r) ∈ R1 , (23) 10: Set h = 0
11: end if
where uST (r) is r-th sample from UST (Eq. 12). The r-th
o 12: Update SSVk+1 = SSVk + h and add incremental
sampled quantity for order oj , denoted by ssj (r) is given
counter k = k + 1
by
13: end while
sosj (r) = soej + uSQ (r) ∈ R1 , (24)
14: return (SSV, ST)
where uSQ (r) is the r-th sample from USQ (Eq. 13). Then
the sampled purchase order arrivals at τ , denoted by sτs (r),
is given by total of Nr realizations of the simulated inventory curves
X are generated (i.e., X :H (r), ∀r ∈ [1, . . . , Nr ]) with the
sτs (r) = sosj (r) ∈ R1 , ∀tosj (r) = τ. (25) forward-looking simulation procedure. The inventory cost
oj
and SL associated with each of the realization of the sim-
The inventory realizations generated by the forward- ulated inventory (i.e., X :H (r)) is computed. An inventory
looking simulation can be expressed by deficit h is calculated to quantify how much lift of the in-
ventory is needed to make sure ψ(X:H ) is greater or equal
X :H (r) = xt (r), . . . , xt+H−1 (r) ∈ RH , to SLmin . For example, in Fig. 8, if the minimum target SL
 
(26) and the SLP are both assumed to be 100%, the inventory
∀r ∈ [1, . . . , Nr ].
deficit h would be the height of the red-boxed area, since
Those realizations X :H (r), ∀r ∈ [1, . . . , Nr ] are used to that is how much the inventory level needs to be lifted so
represent the stochastic variable X:H and will be consumed that 100% (SLP) of the realization curves would satisfy a
by the k-iteration algorithm to produce recommendations minimum target SL of 100%. Note that for the MRP simula-
for the reorder parameters. tions in the k-iteration algorithm, the arrivals in [0, LT ) (LT
represents lead time) is seeded with the simulated arrivals
from the MRP orders placed during the previous sampling
3.2.5 The k-iteration algorithm
time, since the order placed at the simulation of the cur-
The k-iteration algorithm (also included in the reorder pa- rent sampling time will not arrive the till after LT . There-
rameter estimation step shown in Fig. 3) is used to map fore, here the service level is computed in the time period
the hyper-parameters to the reorder parameter recommen- of [LT, min(2LT + ST, H)] due to the fact that the service
dations at each sampling time t. This algorithm is triggered level in [0, LT ) is not in any way affected by the choice of
in both training and validation (live production) phases of SSV at sampling time t, and the service level for later than
the algorithm workflow. 2LT + ST is likely overwritten in the next sampling time.
The procedure of the k-iteration algorithm can be found Once the inventory deficit h is determined, the SSV for
in Algorithm 3. At a sampling time t, the ST is determined the next iteration would be calculated by adding h to the
by the STP and the empirical distribution of the supplier current SSVk (k denotes the iteration). The new SSVk+1
time uncertainty UST in the corresponding USW. For SSV, would be used to run the MRP simulation and the uncer-
it is initialized to 0. Then we enter the loop (i.e., the k- tainty realizations for the next iteration k + 1. As is shown
iteration) for increasing SSV incrementally till the SLmin in Fig. 9, the lifted SSV is usually identical to the lift of
constraint in Eq. 1 is met. In every iteration of this loop, a inventory level since it would result in a large order arrival

12
 in more details in the next section.

3.2.6 Other algorithm treatments

There are a few key treatments of the workflow that would
have significant impact on the performance and outputs of
the algorithm. They include the batched uncertainty real-
izations with the MRP simulation, the “steady-state” initial
inventory of the k-iteration algorithm, the optimization fre-
quency, and the demand forecast pre-processing. We will
explain those treatments in this section.

Batched uncertainty realizations with MRP

In the k-iteration algorithm described above, each uncer-

tainty realization r does not trigger another MRP simula-
tion, which means for Nr simulated inventory curves, the
Figure 8: The k-iteration algorithm: compute inventory MRP simulation is only performed once. The rest of the
deficit h curves are generated by simply perturbing the MRP result
with the uncertainties sampled from the empirical distribu-
right after the lead time (LT) to lift the over all inventory tions as a post-processing step. This process is referred to as
level. The SLP will increase as a result. The iteration would the batched uncertainty realizations with MRP simulation.
continue till the max number of iteration is arrived or there The benefit of such an approach is the reduction of com-
is no need to lift anymore (h = 0), which at each sampling putation time, since the uncertainty sampling is done much
time t, is equivalent to the following expression (Eq. 27) quicker in comparison to the MRP simulation. However,
utilizing the notations developed in Eq. 1 and 5. those realizations are referred to as the proxies since there
could be difference between the results of this approach and
"P
min(2LT+ST,H)
# those results if we sample uncertainty distributions during
i=LT δ(xt+i ) the MRP simulation processes (i.e., sequential uncertainty
γpSL − SLmin ≥ 0, (27)
min(LT + ST, H − LT) realizations with MRP).
The proxies of the simulated inventory is considered to
where pSL is the SLP.
be sufficiently accurate for the realizations of the simulated
inventories, since those realizations are only used to deter-
mine how much SSV should be lifted, rather than accurately
computing inventory cost and SL.

“Steady-state” initial inventory

Intuitively, the initial inventory used for the k-iteration

algorithm at sampling time t is taken from the simulated in-
ventory computed from the previous sampling time, shown
in Eq. 28. Note that at the beginning of the simulation pe-
riod (for either training or validation), the simulated inven-
tory is seeded with the actual inventory.

xt0 = xtsim , (28)

where xt0represents the initial inventory for the k-iteration
algorithm starting at sampling time t. The variable xtsim
represents the inventory simulated at the previous sampling
Figure 9: The k-iteration algorithm: lift SSV time.
The initial inventory determined by Eq. 28 is sensitive to
Note that in the k-iteration procedure, we have config- the simulation process. Specifically, if the simulated inven-
urable choice of the initial inventory, which we will discuss tory is high at t, the initial inventory of the k-iteration would

13
be high. Therefore, all the realizations in the forward- frequency adds another constraint on the optimization of the
looking simulation would be far above zero (satisfying the reorder parameters, we state here that the more frequent the
SLmin constraint) even when we set SSV to 0. As a result, reorder parameters are allowed to change, the more room
the algorithm will give a zero SSV recommendation as long there will be for the algorithm to improve the results. To ra-
as the current simulated inventory is high. tionalize this statement, we setup experimental simulations
Ideally, the customers would expect the SSV recommen- with various optimization frequencies for a total length of
dations be independent of the current inventory level. To 90 days. Fig. 10 shows the SSV and ST recommendations
achieve this goal, we propose an approach that always use with optimization frequencies of 90 days, 30 days, and 15
a “steady-state” inventory as the initial inventory fro the k- days. The blue lines are the SSV and ST recommendations
iteration algorithm. The initial inventory of the k-th itera- with a frequency of 90 days, which remains unchanged dur-
tion is denoted in Eq. 29. ing the entire simulation period. The yellow and red lines
represents the reorder parameters with frequency of 30 and
LT-1
X 15 days, respectively. The changes of lines in those colors
t,t+i t+i
xt0,k = SSVk + (cDF − cA ), (29) are clearly observed.
i=0

where xt0,k represents the initial inventory for the k-th iter-
ation starting at sampling time t. The variable SSVk is the
SSV candidate at the k-th iteration. The variable LT denotes
the lead time. The variable ct,t+i
DF ∈ R1 denotes the demand
forecast for timestamp t + i made at t, and ct+i
A ∈ R1 rep-
resents the simulated arrivals at timestep t + i.
There are two notes for computing the initial inventory
with Eq. 29. First, since the arrivals in [0, LT ) is seeded
with the simulated arrivals, the inventory level at LT in the Figure 10: The SSV and ST recommendations with differ-
k-iteration procedure will be exactly at the level of the can- ent optimization frequencies
didate SSV (i.e., SSVk ) regardless of the simulated inven-
tory level at the sampling time t. Second, Eq. 29 does not Fig. 11 demonstrates the median of the simulated inven-
enforce the initial inventory to be larger or equal to zero, tory levels that correspond to each optimization frequency.
which means negative inventories could occur in [0, LT ) In this example, the simulation inventory level corresponds
due to the choice of this initial inventory. However, since to 90 days is the highest amount all three curves, and the one
the inventory level will be brought back to positive at LT corresponds to 15 days is the lowest, while all three curves
and the service level is only computed in [LT, 2LT + ST ], are 100% above zero. The statement we made holds true in
the potential negative inventories in [0, LT ) will not affect general, since the more frequent the reorder parameters are
the SSV recommendations. allowed to change, the less constraint there will be for the
optimization, and therefore, a better result we would expect.
Optimization frequency To summarize, there is a clear trade-off between the algo-
Another key treatment is the frequency for which the re- rithm performance and the customer adoption for the choice
order parameters are recommended. This frequency also of optimization frequency. In reality, the decision made on
corresponds to the frequency for which Eq. 1 is solved. Ide- this parameter would vary case-by-case.
ally, solving Eq. 1 at every sampling time t would give an
optimal result since new information is usually provided on
a daily basis. However, it would be challenging (if possi-
ble at all) for the customers to adopt reorder parameter rec-
ommendations on a daily basis. Therefore, we designed a
feature in the algorithm which allows the optimization fre-
quency to be configurable by the user.
If the optimization frequency in days is set to be a num- Figure 11: Simulated inventory levels with different opti-
ber j larger than 1 day (e.g., 30 days). The reorder parame- mization frequencies
ter recommendations are updated every j days. For the days
that those reorder parameters are not updated, they will de-
fault to the value computed at the previous sampling time. Demand forecast pre-processing
However, the MRP simulations in the order simulation step One more treatment to be discussed in this section is the
are still performed on a daily basis. Since the optimization pre-processing of the demand forecast data. The demand

14
forecast provided by the customer usually comes with out- Specific to the C3 SIO Application, we leveraged the
liers. Those outliers are either due to input errors or inac- Platform to integrate various formats of datasets from dif-
curate forecasts. The outliers usually are with value that ferent data sources to the Suite such that the data can be
is orders of magnitude larger than the regular data points, easily consumed and orchestrate with the deployed opti-
which could have large impact on the optimization result. mization framework. The computational intensity from pro-
To mitigate the impact of those outliers, we apply a clipping posed framework can be easily identified based on the na-
step on the demand forecast data, which is demonstrated in ture of multiple simulation runs. Typically, customers from
Eq. 30. different industries have hundreds of thousands to millions
types of SKUs to optimize, leading to necessity of effective
 large scale optimization. The Model Deployment Frame-
ηCDF + nc · σCDF ,
 work in C3 AI Suite enable users to scale-out up to mil-
C̃DF, 1D = if CDF, 1D > ηCDF + nc · σCDF , (30) lions of ML models within a single application through pre-


CDF, 1D , otherwise packaged automation including automatic scoring, compute
resource auto-scaling, version management, alerting, hyper-
parameter optimization, and auto retraining. By leveraging
where C̃DF, 1D denotes the clipped demand forecast, and
the model deployment framework on C3 AI Suite, it allows
CDF, 1D is the raw demand forecast data computed in Eq. 15.
us to perform inventory optimization and simulation experi-
The variables ηCDF and σCDF represent the median and the
ments across millions of SKUs, at scale and in a distributed
standard deviation of the demand forecast data, respectively.
fashion.
The variable nc is the a multiplier that controls the strictness
of the clipping. It is a parameter that is configurable by the 4.1. Data integration
user. In our use cases shown in Section 5, we default it to
5.0. The C̃DF, 1D is then used to calculate the convolved de- With the available data integration techniques supported
mand forecast CDF,∗
1D mentioned in Eq. 16.
on the Platform, we have configured the C3 Inventory Op-
The outliers in the raw demand forecast data will also timization Application to directly load the data from the
result in the over-estimation of the demand forecast uncer- enterprise resource planning (ERP) system to facilitate the
tainty. Therefore, the demand forecast uncertainty derived data consumption of the SIO algorithm, as is shown in Fig-
from the clipped demand forecast will go through another ure 13. The C3 Source Collection Types are used to map
step of clipping to further reduce those impacts. The clip- the source data locations and it can be configured to retrieve
ping applied is shown in Eq. 31. the source files on regular basis automatically with only re-
trieving the increments. The C3 Canonical Types are used
( for mapping the table headers of different data formats, and
ηUDF + nu · σUDF , if UDF > ηUDF + nu · σUDF translating them into the data format that is native to the C3
ŨDF = AI Suite. Specifically to the C3 SIO Application, we have
UDF , otherwise
two layers of Canonicals, with the first layer mapping di-
(31)
rectly to the ERP source data, and the second layer mapping
where UDF is the demand forecast uncertainty defined in
∗ from a set that are transformed from the first layer, which
Eq. 18 with CDF derived from C̃DF (instead of from CDF
follow closely to the business logic and are in closer form
directly). The variables ηUDF and σUDF represent the median
to the entity types in the Application. This provides users
and the standard deviation of the demand forecast uncer-
flexibility to provide data in either of the Canonical forms,
tainty, respectively. Similar to that for C̃DF , the variable nu
depending on data availability. The C3 Transform Types
is the a multiplier that controls the strictness of the clipping.
are used for converting the data from Canonical Types and
In our use cases shown in Section 5, we default it to 1.0.
mapping the converted ones either to an additional layer of
Note that the movement uncertainty UMM can be treated in
Canonicals or directly to entity types. Within the C3 Inven-
the same manner to generate clipped movement uncertainty
tory Optimization Application, these entity types hold the
ŨMM .
information associated with certain business modules in re-
source planning, e.g. inventory level, material movement
4. Inventory optimization at scale transactions, purchase orders, production orders, demand
C3 AI Suite (the Platform), serves as the platform that forecasts, etc.
hosts both data and computational frameworks that back the
4.2. Model deployment
C3 Stochastic Inventory Optimization Application, is capa-
ble to provide enterprise users with data integration, ma- By leveraging C3 AI Suite, we first incorporated var-
chine learning development and operation tools, in a gener- ious features based on the feature engineering tools and
alized application-development software stack (Figure 12). developed C3 AI ML pipeline to incorporate the whole

15
 Figure 12: Flowchart for model deployment process in C3 AI Suite

Figure 13: Flowchart for data integration process of ERP data

simulation-optimization work flow. In order to deploy the and use 54 representative SKUs selected by the same client.
ML pipelines on large number of SKUs, we created ML We discuss the impact of some of the model parameters on
segments to support the model to SKUs relationships and the performance of the algorithm. In Section 5.3, we illus-
then run distributed jobs to deploy ML pipelines on a large trate the generality of the algorithm by showing its perfor-
number of SKUs or a group of SKUs that belong to one mance on the customer data across a variety of verticals.
segment. We can also configure champion and challenger
model to make sure the live model is always providing the 5.1. Optimization results for sampled SKUs
optimal results. Model operations and model governance In this part, we apply C3 SIO on two representative
are critical steps in managing the model recommendation SKUs (i.e., SKU A and SKU B) that come from one of our
and explaining the model predictions. clients in the energy and manufacturing industry. The min-
imum target service level for SKU A and SKU B are set to
5. Case study results 91% and 96%, respectively. We have the data available for
both SKUs from 2018-06-01 to 2020-10-01. The training
The C3 Stochastic Inventory Optimization algorithm has period for C3 SIO starts at 2018-06-01 and ends at 2020-
been tested on customer data from a variety of verticals. 03-01, while the validation period starts at 2020-04-01 and
Here, we show the performance of C3 SIO with different ends at 2020-10-01. Since 2020-10-01 is the last date with
levels of details and demonstrate the inventory savings that data available, the operation date (or, live production date)
can be achieved by applying our algorithm. is set to 2020-10-01. The optimization frequency is set to
In Section 5.1, we talk about two representative SKUs 30 days, which means our algorithm will be able to recom-
come from one of our clients in the energy and manufactur- mend new reorder parameters every 30 days.
ing industry. We demonstrate in detail how the algorithm According to the requirements of the client, we did not
manages to recommend a lower SSV while keeping a satis- change STP and always set ST = 0. Therefore, for the
fying service level. In the Section 5.2, we expand the dataset grid search in the training period, we will only select the

16
optimal SLP from the SLP candidate list. The SLP candi- creases of the blue curve in Plot 1. Comparing the blue
date list is set to a default value of [50, 70, 90, 92.5, 95] here. curves and the orange curves (planned arrivals) in Plot 3,
The default multiplier for clipping the demand forecast (nc we will have a sense of the supplier time uncertainty (how
in Eq. 30) is set to 5.0. The default multiplier for clipping much the actual arrivals are late than the planned arrivals)
the demand forecast uncertainty (nu in Eq. 31) is set to 1.0. and the supplier quantity uncertainty (how many quantities
The default minimum length of uncertainty sampling win- the actual arrivals are less than the planned arrivals) for the
dow LUSWmin (Eq. 10) is 30 days and the USW buffer b
USW
actual inventory. Similarly, the dashed red curve in plot 3
(Eq. 10) is set to 14 days. We set uncertainty realizations in corresponds to the simulated arrivals, which corresponds to
the forward-looking simulation procedure (Nr in Eq. 26) to the increases of the dashed orange curve in Plot 1. Com-
100. We will alter those default values later in Section 5.2 to paring the red dashed curve and the green curve, we will
evaluate their impact on the performance of the algorithm. know about the supplier time and quantity uncertainties for
Note that during the order simulation step in the C3 SIO the simulated inventory.
workflow, we also need to generate multiple realizations Plot 4 in Fig. 14 shows the simulated SSV recommenda-
(we note it as Nos here) of the inventory to compare with tion (orange) and the actual SSV (blue). The simulated SSV
the actual inventory the customer has. The variable Nos is is overall much lower than the actual SSV resulting in less
set to 10 in this case. Also, the client intends to focus on order placed in the validation period. This is the main rea-
the savings for the inventory holding cost. Therefore, the son why the overall simulated inventory level is lower than
variable πo,t in Eq. 1 is set to 0 in this case. the actual in Plot 1. Note that the orange curve is stair-like
Figure 14 shows the results for the inventory optimiza- because we set the optimization frequency to 30 days.
tion on SKU A. Plot 1 on the top of Fig. 14 shows the com-
parison between the median of the simulated inventory re-
alizations (dashed orange curve) and the actual inventory
(blue curve) over the validation period. The gray area sur-
rounding the orange curve is the uncertainty range gener-
ated from the 10 uncertainty realizations. Plot 2 shows
the demand forecast (black curve) and the actual consump-
tion (red dashed curve). Plot 3 shows comparisons for the
planned arrivals and the actual arrivals in the simulation and
in actual, respectively. Plot 4 of Fig. 14 compares the simu-
lated SSV recommendations (orange curve) with the actual
SSV (blue curve).
Looking at Plot 1 in Fig. 14, it is obvious that the sim-
ulated inventory in the dashed orange curve is overall at a
lower level in comparison with the actual inventory (blue Figure 14: Simulation result for SKU A with Plot 1-4 from
curve). The gray area shows the uncertainty range for the top to bottom. Plot 1: The simulated inventory and the ac-
simulated inventory with 10 uncertainty realizations. As we tual inventory. Plot 2: The demand forecast and the actual
can see, a small portion of the gray area drops below zero consumption. Plot 3: Planned arrivals and the actual ar-
around June 1st, 2020, indicating a slight out-of-stock risk rivals for both simulated and the actual inventories. Plot 4:
(potentially caused by some outliers in the uncertainties). The simulated SSV recommendations and the actual SSV.
However, the dashed orange curve (the median of all the
curves in the gray area) is far away from 0, which gives us To find the reason behind the low SSV recommenda-
more confidence in the overall healthiness of the inventory tion shown in Plot 4 of Fig. 14, we need to look into the
level. uncertainty sampling process and the k-iteration algorithm.
The black curve and the dashed red curve in the Plot 2 Fig. 15 demonstrates the uncertainties sampled at the begin-
of Fig. 14 indicates that the demand forecast is significantly ning of the validation period. Note that the lead time for this
more than the actual consumption. We refer to this scenario SKU is 6 days, with the USW buffer bUSW set to 14 days,
as the over-forecast. Note that Plot 2 only shows the de- the uncertainty sampling window spams from 2020-03-12
mand forecast made on April 1st, 2020. If we look into the to 2020-04-01. The plot on the left in Fig. 14 shows the his-
demand forecasts on other days in the validation period, we togram of the demand forecast uncertainty. There are more
will find that they are consistently over forecasting. positive than the negative samples in the histogram, which
Plot 3 in Fig. 14 represents the arrival information for is consistent with the over-forecast we observed in Plot 2 of
both simulated inventory and the actual inventory. The ac- Fig. 14.
tual arrival is shown in blue, which corresponds to the in- In the k-iteration algorithm, the demand forecast is used

17
 istics. Fig. 17 shows the result for the inventory optimiza-
tion on SKU B, which is similar to Fig. 14 for SKU A with
all the plots revealing the same information. We want to
point out that there is one big difference for SKU B com-
pared with SKU A regarding the demand forecast. As we
can see in Plot 2 of Fig. 17, the black curve is significantly
less than the dashed red curve, which indicates a significant
under-forecast for the demand forecast. Note that the black
curve in Plot 2 only shows the demand forecast made on
April 1st, 2020. In fact, the demand forecast made on most
of the days in the validation period are all zero. This is not
a rare case in the customer’s data. We observe significant
Figure 15: Sampled uncertainties within the USW ending under-forecast or complete missing of demand forecasts for
at 2020-04-01 for SKU A. From upper left to lower right: the vast majority of the in scope SKUs provided by the cus-
demand forecast uncertainty, movement uncertainty (which tomer, due to the fact that the customer lacks the capability
includes both the miscellaneous and blocked movements), to predict the consumption for the SKUs it holds. Never-
supplier time uncertainty and the supplier quantity uncer- theless, as we will demonstrate with SKU B here, C3 SIO
tainty. algorithm is capable to maintaining the service level with
SSV recommendations even when the demand forecast is
very sparse or completely missing.
as the input for the MRP run. Then, the output of the MRP Plot 4 in Fig. 17 indicates that the algorithm overall rec-
run is “corrected” to the actual consumption by sampling ommends a much higher SSV compared to the actual SSV
the demand forecast uncertainty. If we have more positive consumed by the customer. Yet, Plot 1 shows that that
samples in the demand forecast uncertainty, it means the higher SSV recommendation ends up with a overall lower
actual consumption is more likely to be less than the de- simulated inventory level compared to the actual inventory,
mand forecast used as the input of the MRP. Therefore, the while still satisfies the service level constraint. It is not sur-
simulated inventory level is more likely to be pushed up prising that the actual SSV is much lower than the simulated
in the uncertainty realizations. Fig. 16 shows the simula- SSV, but the actual inventory ends up in a much higher po-
tion with uncertainty realizations in the first iteration of the sition compared to the simulated inventory. Since the actual
k-iteration procedure. The gray curves represent the uncer- inventory is managed by the human operator, and they may
tainty realizations of the simulated inventory. As we can not place order exactly following the MRP output. In the
see, most of the uncertainty realizations are above 0 due to case of SKU B, it is likely that the human operator placed
the reason we discussed above. Therefore, the inventory more orders than what MRP recommends. Should the or-
deficit is small and we do not need to lift too much on SSV ders be placed exactly following the actual SSV, the actual
to satisfy the constraint on the service level. inventory will have significant service level issue (i.e., seri-
ous stock-out).
It is quite intuitive that when the demand forecast is
missing, we maintain a high SSV level to account for the
unexpected consumption. Fig. 18 shows the uncertainties
sampled at the beginning of the validation period. The left
plot with the demand forecast uncertainty indicates an over-
Figure 16: The forward-looking simulation result in the first all under-forecast, since most of the non-zeros samples are
iteration of the k-iteration procedure at 2020-04-01 for SKU negative. The negative samples in the demand forecast un-
A certainty would pull the inventory level lower in the un-
certainty realizations of the k-iteration procedure, forcing
As we can see in this example, C3 SIO algorithm of- a higher SSV to maintain the service level.
fered an objective view for potential risks involved every Fig. 19 shows the uncertainty realizations in the last iter-
time when we recommend the SSV. This view is obtained ation of the k-iteration procedure on 2020-04-01. The gray
from MRP run with multiple uncertainty realizations in the curves in the plot are the uncertainty realizations. Noticing
k-iteration procedure. With SKU A, such a view helped us that the overall downward trend of the inventory is due to
reduced the SSV significantly compared to the actual SSV the negative samples of the demand forecast uncertainty. It
used by the customer. indicates that the demand forecast fed to the MRP run is an
Now let’s look at another SKU with different character- under-forecast, so the “corrections” are applied to pull the

18
 5.2. Sensitivity study on the model parameters
Different from the hyper-parameters, which specifically
refer to SLP and STP. There are a few model parameters
which are also critical to the performance of C3 SIO. Those
model parameters include the multiplier for clipping the de-
mand forecast (i.e., the variable nc in Eq. 30), the multiplier
for clipping the demand forecast uncertainty (i.e., the vari-
able nu in Eq. 31), the candidate list for SLP, the length
of the offset (buffer) for USW (i.e., the variable bUSW in
Eq. 10), and the number of uncertainty realizations (i.e., the
variable Nr in Eq. 26).
To evaluate the impact of the model parameters on the
model performance, we expand the experiment subjects to
Figure 17: Simulation result for SKU B with Plot 1-4 from an ensemble of 54 SKUs, which comes from same client
top to bottom. Plot 1: The simulated inventory and the ac- as the SKUs shown in Section 5.1. Those 54 representative
tual inventory. Plot 2: The demand forecast and the actual SKUs are hand-picked by the client’s subject matter experts
consumption. Plot 3: Planned arrivals and the actual ar- with an intention to cover various scenarios they encounter
rivals for both simulated and the actual inventories. Plot 4: in the business operation.
The simulated SSV recommendations and the actual SSV. For each of the study, we compare a few key perfor-
mance metrics including the percentage of SKUs meeting
the target service level (i.e., service level adherence rate,
rAD ), inventory holding cost saving percentage averaged
over the entire validation period (s̄inv ), safety stock holding
cost saving percentage averaged over the entire validation
period (s̄ss ), safety stock holding cost saving percentage at
the operation date (sss, op ).
The service level adherence rate is defined as follows:
PNSKU
i=1 δ(SLi − SLi,target )
rAD = , (32)
NSKU
where variable NSKU represents the total number of SKUs
in the scope. In our experiment setting, NSKU = 54. The
notation δ(·) is a level set function defined in Eq. 6. The
Figure 18: Sampled uncertainties within the USW ending
variables SLi and SLi,target denotes the service level from
at 2020-04-01 for SKU B. From upper left to lower right:
C3 SIO simulation and the target service level for SKU i,
demand forecast uncertainty, movement uncertainty (which
respectively.
includes both the miscellaneous and blocked movements),
The inventory holding cost saving percentage averaged
supplier time uncertainty and the supplier quantity uncer-
over validation period is defined as follows:
tainty.
PNT t
(xa − xtsim )πht
s̄inv = t=1 PNT t t , (33)
t=1 xa πh
inventory level down. Therefore, to meet the service level
constraint, a higher SSV is recommended. where the variable t denotes the timestep in the validation
period, while NT represents the total number of days in the
validation period. The variables xta and xtsim depict the ac-
tual and the simulated inventory level at timestamp t, re-
spectively, and πh represents the inventory holding cost.
Similarly, the safety stock holding cost saving percent-
age averaged over validation period, and the safety stock
Figure 19: The forward-looking simulation result in the fi- holding cost saving percentage at operation date are defined
nal iteration of the k-iteration procedure at 2020-04-01 for in Eqs. 34 and 35 below.
SKU B PNT
(SSVt − SSVtsim )πht
s̄ss = t=1 PNaT , (34)
t t
t=1 xa πh

19
 op op
Fig. 20. We compared cases with nu set to 0.25, 0.5 and
NT (SSVop
a − SSVsim )πh 2.0. As we increase the nu , the rAD will improve at the
sss, op = PNT t t , (35)
t=1 xa πh
cost of reducing savings. As nu gets larger than 1.0, the im-
provement for rAD stopped. The results here indicates that
where the notation SSV represents the safety stock value,
nu = 1.0 is a reasonable choice, if improving rAD is more
and the superscript op represents the timestep at the opera-
important increasing savings.
tion date, which usually is the date right after the validation
period finishes. Note that in Eqs. 34 and 35, the denomina-
tor is the total inventory holding cost (rather than the total
safety stock holding cost) over the validation period. This
is worth to note that s̄ss and sss, op are more important to this
client compared to s̄inv in this case.
The results of the sensitivity study for nc (the multiplier
for clipping the demand forecast) is shown in Fig. 20. The
we set nc = 5.0 as the baseline case, and experimented with
both nc larger and smaller than the baseline. The x-axis
shows nc of 1.0, 2.0, 5.0, and 10.0, successively. The blue
bars represent rAD , the orange bars represent s̄inv , the green
bars represent sss, op , and the red bars represent s̄ss . With nc
larger than the baseline, we see a clear dropping for all the Figure 21: Sensitivity study on the clipping multiplier for
saving metrics (i.e., s̄inv , s̄ss , sss, op ) with no improvement demand forecast uncertainty and movement uncertainty
at rAD . On the other hand, if we make the clipping criteria
more strict and have nc smaller than the baseline, we can The hyper-parameter search results depends on the ini-
observe a clear improvement of the saving metrics, espe- tial candidate list we provided to the algorithm. For exam-
cially for the inventory holding cost savings (orange bars). ple, if the grid search of SLP starts with different candidate
However, the service level adherence rate will be lower. SLP lists, the optimized SLP will be different and the re-
Obviously, there is a trade-off for between the service sulting inventory or SSV savings would be affected. Fig. 22
level and the saving. Starting from nc = 1.0, if we aim demonstrates the results for the sensitivity study on the SLP
to improve rAD by increasing nc , the improvement would candidate lists, since in this study we only optimizing on
come at the cost of reducing the savings. However, after SLP as hyper-parameter. Each item on the x-axis notations
nc = 5.0, increasing nc will no longer provide improve- a SLP candidate list. The mapping between the x-axis no-
ment on rAD . Therefore, with the presented data, we would tation in Fig. 22 and the SLP candidate list can be found
stop at nc = 5.0 to avoid additional reduction on the sav- in Table 1. For example, “Baseline” corresponds to the ini-
ings. tial SLP candidates of 50%, 70%, 90%, 92.5%, 95%. For
experiments with “Baseline”, the best SLP will be chosen
from the above 5 candidates.

Notation SLP candidates

Baseline [50, 70, 90, 92.5, 95]
SLP5 [50, 55, 60, 65, 70, 75, 80,
85, 90, 92.5, 95]
SLPMax975 [50, 70, 90, 92.5, 95, 97.5]
SLPMax100 [50, 70, 90, 92.5, 95, 97.5, 100]
SLPMin10 [10, 30, 50, 70, 90, 92.5, 95]
SLPMin10Max975 [10, 30, 50, 70, 90, 92.5, 95, 97.5]
Figure 20: Sensitivity study on the clipping multiplier for
raw demand forecast data Table 1: SLP candidates shown in Fig. 22

Similar study can be applied to nu (the multiplier for Figure 22 clearly indicates that the service level adher-
clipping the demand forecast uncertainty). Fig. 21 com- ence rate is not affected as we change the SLP candidate
pares the service level adherence rate and the saving met- list, for the 54 SKUs tested. In addition, as we increase the
rics for inventory optimization results with various nu . The upper bound of the candidate list from 95% to 97.5% and
baseline for nu is 1.0. The bars map to the same metrics as 100%, the saving metrics drop significantly. It is expected

20
since a high SLP value means the recommended SSV will
be more likely to be affected by the outliers of the uncer-
tainties. For example, if we have an outlier in the supplier
quantity uncertainty (a large shortage), setting SLP equal
to 100% means the recommended SSV will ensure the in-
ventory meets target service level even when this outlier is
present, which will end up with a very high SSV recom-
mendation. On the other hand, if the lower bound of the
candidate list is extended, we would see a slight increase
of the savings. Since selecting a very low SLP candidate
involves larger risks for service level hits, the decision of
Figure 23: Sensitivity study on the buffer for the uncertainty
extending the lower bound or not will be made by the cus-
sampling window (USW) size
tomers based on their risk profiles.

Figure 24: Convergence study on the number of uncertainty

Figure 22: Sensitivity study on the candidate list for service
realizations in the forward-looking simulation procedure
level percentile (SLP)

Figure 23 shows the results for cases with various buffers results of the algorithm are summarized in Tab. 2. The
for the USW (bUSW ). The baseline buffer is set to 14 days. columns in Tab. 2 refers to the total number of SKUs in-
Then the results with buffer of 30, 90, 180 and 365 days scope, the actual inventory cost for the customer for the val-
are plotted successively in the figure. Despite some minor idation period, the simulated inventory cost generate by C3
inconsistencies between the inventory saving and the SSV SIO in the validation period, the actual SSV cost for the cus-
saving, we observe a general downward trend for both sav- tomer, the cost of SSV recommended by C3, the percentage
ing metrics and the service level adherence rate rAD . This of savings (either with inventory savings or SSV savings),
result is expected since as the USW gets longer, outliers the actual service level maintained by the customer, and the
are more likely to appear in the sampling window, causing service level for the simulated inventory (with SSV recom-
problems for both savings and service level. Therefore, it mended by C3).
is a wise choice to keep a shorter USW as long as enough Note that for some of the customers (i.e., the one in en-
samples could be collected within the window. ergy & manufacturing vertical), they measure the savings
Figure 24 shows a convergence study on the number of by the percentage reduction on the SSV cost (defined in
uncertainty realizations in the forward-looking simulation Eq. 34). We want to distinguish it from other inventory sav-
(Nr ). The baseline in the figure represents 100 realizations. ings while interpreting the results.
As is observed from the figure, results experience large os- We demonstrated in Tab. 2 that C3 SIO algorithm is ef-
cillation with smaller number of realizations. As the num- fective for large-scale inventory optimization across a vari-
ber of realizations approaches 500, both saving metrics and ety of verticals. It is capable of reducing the inventory and
the service level adherence rate rAD stabilize. the safety stock costs by 10-20% while maintaining the ser-
5.3. Summary of results for customers in different vice level similar to its current position.
verticals
6. Conclusion and discussions
In addition to the customer in the energy and manu-
facturing industry, we have applied C3 SIO algorithm to In this paper, we introduced the challenges for the large-
other customer data in a variety of verticals including build- scale supply chains, and proposed to solve them with C3
ing, manufacturing, healthcare, and energy in general. The Stochastic Inventory Optimization (SIO) algorithm. The C3

21
 Customers #SKUs Inva InvC3 SSVa SSVC3 Saving % SLa SLC3
Customer 1 501 $10M $7.5M - - 25% 99% 98.0%
Customer 2 20.6k $64M - $11M $1.9M 14.2%* 99.9%** 98.7%
Customer 3 1.4k $24M $21M - - 14% 88.6% 90.3%
Customer 4 1.3K $7.7M $6.5M - - 15% 97.6% 91.0%
Customer 5 1k $11.7M $9.2M - - 21.2% 98.6% 94.2%

Table 2: Inventory optimization results for different customers (k: thousand, M: million) with EUR converted to USD at a
rate of 1.1. *SSV saving percentage is calculated following Eq. 34. **The actual SL is unreasonably high due to human
intervention with extra information

SIO algorithm works well with a variety of MRP systems, lize the classical solutions and further improve C3 SIO.
and we demonstrated in detail how it works with a specific
type of MRP system (i.e., the safety stock MRP). The C3 Acknowledgements
SIO algorithm has a training phase that helps the customer
to decide on the hyper-parameters that best describe their We would like to thank [Link] who provided all the re-
risk profile. The algorithm also has a validation (live pro- sources that made this research possible. We would like to
duction) phase that recommends the reorder parameters in thank our customers to worked with us over the last few
real-time. years, providing valuable feedback to help improving our
C3 Inventory Optimization application. We would also like
Note that the algorithm allows the customer to adjust
to thank our internal experts Alex Amato, Nikhil Krish-
the objective function based on their business situation.
nan, Burak Gundogdu for their valuable inputs. Our special
Also, the customers have the option to bypass the training
thanks go to Stefano Zavagli who provided insightful com-
phase and enter the live production with predefined hyper-
ments and contributed to some of the figures in the paper.
parameters if they already have a good understanding of
their risk profile. We pointed out that some idealized as-
sumptions for the uncertainties in the literature oftentimes Nomenclature
does not hold. The C3 SIO algorithm illustrated a practical Acronyms
workflow to handle the uncertainties in real-time produc-
tion. BOM Bill of materials
The performance of C3 SIO algorithm is demonstrated in
EA Expected expedited arrival
detail with data provided by one of the clients in the energy
and manufacturing industry. We also summarize the perfor- ELT Expedited lead time
mance of the algorithm on customers’ data across a variety
of verticals. We showed that C3 SIO algorithm is effective ERP Enterprise resource planning
in reducing the inventory and the safety stock costs while
maintaining a satisfying service level. FLS Forward-looking simulation
A sensitivity study is conducted showing that some of the KIA The k-iteration algorithm
model parameters will impact the optimization results. The
sensitivity study provided reasonable choices for model pa- LT Lead time
rameters as default values while deploying the algorithm in
scale. Nevertheless, we still encourage the customers to test MO Minimum order quantity
on a sub-sample of SKUs to find the best model parameters
MRP Material requirements planning
for their dataset before the live production.
Finally, we want to note here that C3 SIO algorithm is PTF Planning time fence
distinct from the solutions built with MILP in the litera-
ture due to its compatibility with the MRP systems, which RV Rounding value
carries critical value in terms of the customer adaptation.
SA Expected standard arrival
Should the MRP limitation be removed and the customers
are open to adapt an end-to-end solution (i.e., instead of ac- SIO Stochastic inventory optimization
cepting reorder parameters, but rather the order placement
recommendations), we would have more flexibility to uti- SKU Stock keeping unit

22
SL Service level H Predictive horizon

SLP Service level percentile LUSW

min The minimum length of the uncertainty sampling
window
SSV Safety stock value
Nr Number of uncertainty realizations in the forward-
ST Safety time
looking simulation procedure
STP Safety time percentile
Nos Number of uncertainty realizations in the order simu-
USW Uncertainty sampling window lation step

Variables pSL Service level percentile (SLP)

Θ Set of stochastic variables r The reorder parameters for one day
δ Level set function R:τ The 1d array of reorder parameters for a period of time,
γp Aggregation function that takes the p-th percentile of the which starts at t (inclusive) and ends at t + τ (exclusive)
uncertainty distribution
UDF The 1d array to represent the demand forecast uncer-
a The stochastic expected order arrival for one day (scalar) tainty within the uncertainty sampling window

Aτ :H The 1d array of the stochastic expected order arrivals UMM The 1d array to represent the material movement un-
starting at timestamp t + τ (inclusive) and ends at t + H certainty within the uncertainty sampling window
(exclusive)
USQ The 1d array to represent the supplier quantity uncer-
x The stochastic inventory level for one day (scalars) tainty within the uncertainty sampling window

X:τ The 1d array of stochastic inventory level for a period UST The 1d array to represent the supplier time uncertainty
of time, which starts at t (inclusive) and ends at t + τ (ex- within the uncertainty sampling window
clusive)

πh Inventory holding cost

References
Agarwal, A. (2019). Multi-echelon supply chain inventory plan-
πo,t Inventory ordering and transportation cost ning using simulation-optimization with data resampling. arXiv
preprint arXiv:1901.00090.
ψ Function that calculates the service level
Aikens, C. (1985). Facility location models for distribution plan-
Ã:H The 1d array of the stochastic expected order arrival
ning. European Journal of Operational Research, 22(3):263–
for the current predictive horizon predicted at a prior sam- 279.
pling time t̃. The array starts at t (inclusive) and ends at
t + H (exclusive) Chu, Y., You, F., Wassick, J. M., and Agarwal, A. (2015).
Simulation-based optimization framework for multi-echelon in-
bUSW The buffer for the length of the uncertainty sampling ventory systems under uncertainty. Computers & Chemical En-
window gineering, 73:1–16.

cτA Actual consumption at timestamp τ Fu, M. C. (2002). Optimization for simulation: Theory vs. prac-
tice. INFORMS J. on Computing, 14(3):192–215.
CDF, 1D Demand forecast aggregated to 1-dimension (from
the off-diagonal elements in the 2d matrix) within the un- Garcia, C. E., Prett, D. M., and Morari, M. (1989). Model pre-
certainty sampling window dictive control: Theory and practice—a survey. Automatica,
25(3):335–348.
CDF, 2D The 2-dimensional demand forecast matrix within
the uncertainty sampling window Geoffrion, A. and Graves, G. (2010). Multicommodity Distribution
System Design By Benders Decomposition, volume 20, pages
τ
CDF Demand forecast (1d array) made at the date τ for a 35–61.
length of predictive horizon H
Geoffrion, A. and Powers, R. F. (1995). Twenty years of strategic
t ,t
o f
cDF Demand forecast made at the date to for the date tf distribution system design: An evolutionary perspective. Inter-
(scalar) faces, 25:105–127.

23
Gupta, A. and Maranas, C. (2000). A two-stage modeling and Vidal, C. J. and Goetschalckx, M. (1997). Strategic production-
solution framework for multisite midterm planning under de- distribution models: A critical review with emphasis on global
mand uncertainty. Industrial & Engineering Chemistry Re- supply chain models. European Journal of Operational Re-
search, 39:3799–3813. search, 98(1):1–18.

Ierapetritou, M. and Pistikopoulos, E. (1994). Novel optimization Wagner, H. M. and Whitin, T. M. (1958). Dynamic version of the
approach of stochastic planning models. Industrial & Engineer- economic lot size model. Management science, 5(1):89–96.
ing Chemistry Research, 33:1930–1942.
Zipkin, P. H. (2000). Foundations of inventory management.
Joines, J. A., Gupta, D., Gokce, M. A., King, R. E., and Kay, M. G.
(2002). Supply chain multi-objective simulation optimization.
In Proceedings of the Winter Simulation Conference, volume 2,
pages 1306–1314 vol.2.

Jung, J. Y., Blau, G., Pekny, J. F., Reklaitis, G. V., and Eversdyk,
D. (2008). Integrated safety stock management for multi-stage
supply chains under production capacity constraints. Comput-
ers & Chemical Engineering, 32(11):2570–2581. Enterprise-
Wide Optimization.

Köchel, P. and Nieländer, U. (2005). Simulation-based optimisa-

tion of multi-echelon inventory systems. International Journal
of Production Economics, 93-94:505–513. Proceedings of the
Twelfth International Symposium on Inventories.

Mayne, D. Q., Rawlings, J. B., Rao, C. V., and Scokaert, P. O.

(2000). Constrained model predictive control: Stability and op-
timality. Automatica, 36(6):789–814.

Mele, F. D., Guillen, G., Espuna, A., and Puigjaner, L. (2006).

A simulation-based optimization framework for parameter op-
timization of supply-chain networks. Industrial & Engineering
Chemistry Research, 45(9):3133–3148.

Olafsson, S. and Jumi Kim (2002). Simulation optimization. In

Proceedings of the Winter Simulation Conference, volume 1,
pages 79–84 vol.1.

Ptak, C. and Smith, C. (2011). Orlicky’s material requirements

planning. McGraw-Hill Education.

Santoso, T., Ahmed, S., Goetschalckx, M., and Shapiro, A. (2005).

A stochastic programming approach for supply chain network
design under uncertainty. European Journal of Operational Re-
search, 167(1):96–115.

Silver, E. A. (1973). A heuristic for selecting lot size quantities

for the case of a deterministic time-varying demand rate and
discrete opportunities for replenishment. Prod. Inventory Man-
age., 2:64–74.

Smith, D. A. and Smith, C. (2014). Demand driven performance:

using smart metrics. McGraw-Hill Education.

Swaminathan, J. M., Smith, S. F., and Sadeh, N. M. (1998). Mod-

eling supply chain dynamics: A multiagent approach. Decision
sciences, 29(3):607–632.

Velasco Acosta, A. P., Mascle, C., and Baptiste, P. (2020). Ap-

plicability of demand-driven mrp in a complex manufacturing
environment. International Journal of Production Research,
58(14):4233–4245.

24