Simulating Single-Item Inventory Systems

1. Simulating Single-Item Inventory Systems

In this section, we develop a simulation environment to analyze a single-item inventory system under periodic review controlled by an \((s,S)\) policy.

The Lighthouse Company sells ceiling lamps and, since the replenishment cost \(K=50\) euro is relatively high, it considers using an \((s,S)\) policy to control the inventory. The company wants to use simulation to find sound values for \(s\) and \(S\).

The details of the case are as follows. Ceiling lamps sell for 100 euro, while the purchase price is 40 euro. The monthly holding cost is estimated at 50% of the purchase price. To prevent customers from buying from a competitor (e.g., online), the Lighthouse Company offers a discount of 20% per day on the regular selling price when demand cannot be met from on-hand stock. The pmf \(f_{i} = \P{D=i}\) of daily demand is given by

\begin{align*} f_{0} &= 1/6, & f_{1} &= 1/5, & f_{2}&= 1/4, & f_{3} &= 1/8, & f_{4} &= 11/120, & f_{5} &= 1/6. \end{align*}

The lead time \(\leadtime=2\) days.

Here are the model parameters. The daily demand is a RV. The cost parameters must have consistent units; here we convert costs to a per-day basis.

"""Lighthouse case: model parameters."""
import random_variable as rv

D = rv.RV({0: 1 / 6, 1: 1 / 5, 2: 1 / 4, 3: 1 / 8, 4: 11 / 120, 5: 1 / 6})

l = 2  # lead time
h = 40 * 0.5 / 30  # daily holding cost
b = 100 * 0.2  # daily backlog cost
K = 50  # Order cost, used for sS and Qr
N = 100  # simulation duration

if N <= l:
    print(f"The simulation duration {N=} is shorter than the lead time {l=}.")
    quit()

With the code below, we set up a simulation to estimate the performance measures for one specific \((s,S)\) policy. For the KPIs we introduce two general helper functions, and write these to the file functions.py. By importing this file, we can use these functions at other places, just like the RV class.

import numpy as np


def Plus(x):
    if isinstance(x, (int, float)):
        return max(x, 0)
    elif isinstance(x, (list, tuple)):
        return [Plus(v) for v in x]
    else:
        return np.maximum(x, 0)


def Min(x):
    if isinstance(x, (int, float)):
        return Plus(-x)
    elif isinstance(x, (list, tuple)):
        return [Min(v) for v in x]
    else:
        return Plus(-x)

We now proceed with the inventory model. Most of the following modules should be familiar. The last module loads the RV class, which was discussed in \cref{sec:prob_arithmetic}.

import numpy as np
from numpy.random import default_rng

import random_variable as rv
from  functions import Plus, Min

from lighthouse_case import D, K, N, l, h, b

This function computes, for the demand D, the inventory positions and inventory levels under specified s and S. For the inventory levels, we combine \cref{eq:i5} into one recursion. The first for loop for the inventory levels computes the levels for \(t < \leadtime\), since no replenishments occur in those periods.

def simulate_sS(demands, s, S):
    P = np.zeros(N)  # inventory position
    Q = np.zeros(N)  # order quantities
    IL = np.zeros(N)  # inventory level
    P[0] = IL[0] = S  # initial values

    for t in range(1, N):
        Pprime = P[t - 1] - demands[t - 1]
        Q[t] = (S - Pprime) * (Pprime <= s)
        P[t] = Pprime + Q[t]

    for t in range(1, l):
        IL[t] = IL[t - 1] - demands[t - 1]
    for t in range(min(l, N), N):
        IL[t] = IL[t - 1] - demands[t - 1] + Q[t - l]

    return P, IL, Q

For the service levels, the cycle service level \(\alpha_{c}\) is of particular interest; we compute it using slicing, a common Python technique.¹¹ Ask ChatGPT for explanation if you find this line hard. I found it hard to get it right.

def cost(IL, Q):
    return K * (Q > 0).mean() + h * Plus(IL).mean() + b * Min(IL).mean()


def alpha(IL):
    return (IL >= 0).mean()


def beta(demands, IL):
    return np.minimum(demands, Plus(IL)).sum() / demands.sum()


def alpha_c(IL, Q):
    return ((Q[:-l] > 0) * (IL[l - 1 : -1] >= 0)).sum() / (Q[:-l] > 0).sum()

We can now run a simulation to compute the KPIs.

rng = default_rng(3)  # set the seed
demands = D.rvs(size=N, random_state=rng)

P, IL, Q = simulate_sS(demands, s=3, S=20)

print(f"{demands.mean()=:.2f}, {P.mean()=:.2f}, {IL.mean()=:.2f}, {Q.mean()=:.2f}")
print(f"{Plus(IL).mean()=:.2f}, {Min(IL).mean()=:.2f}")
print(f"{cost(IL, Q)=:.2f}")
print(f"{alpha(IL)=:.2f}, {beta(demands, IL)=:.2f}, {alpha_c(IL, Q)=:.2f}")

demands.mean()=2.32, P.mean()=13.11, IL.mean()=8.64, Q.mean()=2.32
Plus(IL).mean()=8.71, Min(IL).mean()=0.07
cost(IL, Q)=13.71
alpha(IL)=0.95, beta(demands, IL)=0.88, alpha_c(IL, Q)=0.67

It remains to find good policy parameters. In principle, this is simple: vary s and S over a range of values to find the cost-minimizing combination.²² Finding the optimal \(s\) and \(S\) efficiently is not an easy problem. If you are interested, consult me for details. For this we can use the next general grid search algorithm, which we store in function.py; below we apply it to our case

def grid_search(grid, f):
    "Minimize a f on the grid."
    best_x, best_value = None, float('inf')
    for x in grid:
        if not isinstance(x, (tuple, list)):
            # allow unpacking *x for non-iterable x
            x = (x,)
        value = f(*x)
        if value < best_value:
            best_x, best_value = x, value
    return best_x, best_value

The results show that the best \(s\) and \(S\) lie in the search area. This suggests they are optimal.

def cost_of(demands, s, S):
    P, IL, Q = simulate_sS(demands, s, S)
    return cost(IL, Q)

from functions import grid_search

pairs = ((s, S) for s in range(-5, 10) for S in range(s + 1, 30))
(s_opt, S_opt), V = grid_search(pairs, lambda s, S: cost_of(demands, s, S))

print(f"{s_opt=}, {S_opt=}, {V=:.2f}")

s_opt=5, S_opt=25, V=13.00

Exercise 1:

What inventory policy does the following code try to implement; what is wrong; and what are its policy parameters?

import numpy as np

gen = np.random.default_rng(seed=42)

l = 2
n = 100
Q_size = 10
s = 20
labda = 3
D = gen.poisson(labda, size=n)

P = np.zeros_like(D)
Q = np.zeros_like(D)
for t in range(1, len(D)):
    Q[t] = Q_size * (P[t - 1] <= s)
    P[t] = P[t - 1] + Q[t] - D[t]

IL = np.zeros_like(D)
for t in range(l, len(D)):
    IL[t] = IL[t - 1] + Q[t - l] - D[t]

Solution

Did you actually try? Maybe see the ‘hints’ above!:

Solution, for real

The \((Q,r)\) policy with \(Q=10\) and \(r=s=20\). However, this version is correct; the loop in the problem statement is wrong.

for t in range(1, len(D)):
    Pprime = P[t - 1] - D[t - 1]
    Q[t] = Q_size * (Pprime <= s)
    P[t] = P[t - 1] + Q[t] - D[t]

This work is licensed by the University of Groningen under a Creative Commons Attribution-ShareAlike 4.0 International License.