A Simulation of the Base-Stock Policy

Here we simulate an inventory process subjected to a base-stock policy with a numpy device calles structured arrays.

Orders arrive according to a Poisson process with rate \(\lambda\) at an inventory. When the inventory level is positive the order is satisfied, otherwise the order is backlogged. For every order a replenishment order is issued to replenish the inventory at a lead time \(L\) later. The order-up-to level is \(S\).

All code in one file

The Simulator

import numpy as np
import pylab as pl
from scipy.stats import expon, poisson


N = 1e6
labda = 3.2
leadtime = 1.2
S = 3

N is the simulation length, that is, the number of demands.

dtype = [('time', float), ('amount', int)]

We implement the demands, replenishments and the inventory level as structured arrays. The first field is the time, the second the ordered amount. Therefore we use the above dtype.

demands = np.zeros(N, dtype = dtype)
interarrivalTimes = expon.rvs(size = N, scale = 1./labda)
arrivaltimes = np.cumsum(interarrivalTimes)
demands['time'] = arrivaltimes
demands['amount'] = -1

We first generate the inter-arrival times. By adding these, by means of cumsum, we get the arrival times of the orders themselves. Since each demand reduces the inventory level by 1, we set the amount to -1.

replenishments = np.zeros(N, dtype = dtype)
replenishments['time'] = demands['time'] + leadtime
replenishments['amount'] = 1

Each replenishment arrives a lead time after the demand that triggered the replenishment. Since replenishments increase the inventory by 1, the amount is set to 1.

epochs = np.sort(np.concatenate((demands, replenishments)), order = 'time')

An epoch is a point in time in which the inventory level changes by 1 or -1. We therefore sort the epochs at which the demands and the replenishments arrive in increasing order of the time.

init = np.zeros(1,dtype=dtype)
init[0] = (0,S)
inventory = np.concatenate((init, epochs))
inventory['amount'] = np.cumsum(inventory['amount'])

The inventory level starts at level S at time 0. Some thought will reveal that the cumulative sum over the amounts results in the inventory level.

A sample path:

jj = 100
pl.plot(inventory['time'][:jj], inventory['amount'][:jj])


We now compute some probabilities.

m =  min(inventory['amount'])

For the purpose of binning the results we need the lowest level m of the inventory.

We first compute the fraction of time spent the inventory level spends at a certain level.

bins = np.zeros(S+1-m)

prev = inventory[0]
for epoch in inventory[1:]:
    bins[epoch['amount']] += epoch['time'] - prev['time']
    prev = epoch

bins = np.roll(bins, -m)
bins = bins/np.linalg.norm(bins, 1) # normalize

We use the inventory level as an index to the bin, and sum the amount of time spent at each inventory level. We need to roll the elements of bins to shift the elements such that the bin corresponding to inventory level m is the left-most array element, and the bin corresponding to level S is the right-most element.

Next, we compute use histogram to compute the number of epochs that see the inventory at a certain level. Note that these epochs include demand arrivals and replenishments.

hist = np.histogram(inventory['amount'], bins = np.arange(m-0.5, S+1.5), \
                    density = True)

As np.histogram requires the boundaries of the bins, we use boundaries in between the integers.

Finally, theory tells us that \(P(I = S-n) = P(D= n)\), i.e., the steady state probability that the inventory level \(I = S-n\) must be equal to the probability that the demand \(D\) during the lead time is equal to \(n\). Since the interarrival times are exponential, the number of arrivals during the lead time is Poisson. We therefore compute

\[P(D=n ) = \frac{(\lambda L)^n}{n!}\exp(-\lambda L)\]
probs = poisson.pmf(range(S+1-m), labda*leadtime)
pl.plot(range(m,S+1), bins, label = 'B')
pl.plot(range(m,S+1), hist[0], label = 'S')
pl.plot(range(m,S+1), probs[::-1], label = 'P')

It is interesting to think about why the plot with label S appears to be shifted to the left. Is there something wrong with the PASTA property?

Epoch type given the inventory level

What can we say about the type of the arrival epoch given the inventory level? Clearly, when \(I_t=S\) at time \(t\), only demands can arrive, so that \(P(T=D|I=S) = 1\), where \(T\) is the type of the epoch and \(D\) denotes a demand. Also, when \(I_t= m\), where \(m\) is the minimal inventory level as observed during the simulation, it must be that \(P(T=R|I=m) = 1\), where \(R\) denotes a replenishment epoch.

We now estimate \(P(T=D|I)\).

bins = np.zeros((S+1-m,2))

prev = inventory[0]
for epoch in inventory[1:]:
    bins[prev['amount'], (epoch['amount'] - prev['amount']+1)/2] += 1
    prev = epoch
bins = np.roll(bins, -m, axis = 0)

We map a demand, which occurs when epoch[‘amount’] - prev[‘amount’] = -1, to an index 0, and a replenishment, which occurs when epoch[‘amount’] - prev[‘amount’] = 1, to an index 1.

bins = [b[0]/sum(b) for b in bins] # normalize
pl.plot(range(m,S+1), bins, label = 'Type is D')