The Base-stock Model

All code in one file

This code illustrates the python implementation of the base-stock model as discussed in [hopp08:_factor_physic]. For the formulae I refer to Addendum to Chapter 2 of Factory Physics: Formulae.

import numpy as np
from scipy import stats

The basestock class

class BaseStock(object):
    def __init__(self, l, h, b,X):
        self.l = l # leadtime
        self.h = h # holding cost
        self.b = b # backorder cost
        self.X = X # demand distribution during lead time

        self.theta = self.X.mean()
        self.sigma = self.X.std()
        self.G = self.X.cdf

    def S(self,r): # service level
        return self.G(r)

    def B(self, r):  # backorder level
        return self.X.expect(lambda x: x - r - 1, r + 1, np.inf)

    def I(self,r):
        return r + 1. - self.theta + self.B(r)

    def Y(self,r):
        return self.h*self.I(r) + self.b*self.B(r)

    def rstar(self):
        return self.X.ppf(self.b/(self.b+self.h))

    def meetServiceLevel(self, level):
        # compute reorderlevel for a given service level
        return self.X.ppf(level)

Poisson demand

Let us consider the example of Superior Appliances of Section 2.4.2. The demand is Poisson distributed with mean \(\theta = 10\). Therefore,

theta = 10
l = 1 # leadtime
h = 15.
b = 25.
X = stats.poisson(theta)
basestock = BaseStock(l=l, h=h, b=b, X=X)

r = 14

The service level for this value of \(r\) is

print("Service level: ", basestock.S(r))
Service level:  0.916541527065

We can check this from the definition too.

for i in range(4, 18):
    print(r"g(%d) = %f G(%d) = %f"%(i, X.pmf(i), i, X.cdf(i)))
g(4) = 0.018917 G(4) = 0.029253
g(5) = 0.037833 G(5) = 0.067086
g(6) = 0.063055 G(6) = 0.130141
g(7) = 0.090079 G(7) = 0.220221
g(8) = 0.112599 G(8) = 0.332820
g(9) = 0.125110 G(9) = 0.457930
g(10) = 0.125110 G(10) = 0.583040
g(11) = 0.113736 G(11) = 0.696776
g(12) = 0.094780 G(12) = 0.791556
g(13) = 0.072908 G(13) = 0.864464
g(14) = 0.052077 G(14) = 0.916542
g(15) = 0.034718 G(15) = 0.951260
g(16) = 0.021699 G(16) = 0.972958
g(17) = 0.012764 G(17) = 0.985722

If we want a service level of at least \(0.9\), then we see that this corresponds to \(r^* = 14\) as, from the numbers above, we read that \(G(13) < 0.9 <G(14)\).

We can also apply the critical fractile formula to get the same result.

print("Minimal reorder level for a fillrate of 0.9: ",
basestock.meetServiceLevel(0.9))
Minimal reorder level for a fillrate of 0.9:  14.0

Now we turn to the average backlog level.

print("Expected backlog: ", basestock.B(r))
Expected backlog:  0.103478679787

Let us also compare this value with the other formulas we found in Addendum to Chapter 2 of Factory Physics: Formulae. Observe that the summation in the first of these equations runs from level \(r+1\) to \(\infty\). Clearly, this summation involves an infinite number of terms, which is somewhat problematic to handle numerically. There are several ways to get around this problem, but the easiest way is to simply truncate the summation when the terms becomes very small. Here we simply truncate at the length 100.

print(sum(X.sf(i) for i in range(r+1, 100)))
print(theta - sum(X.sf(i) for i in range(0, r+1)))
0.103478679787
0.103478679787

We get the same results, as expected.

The average on-hand inventory level is this:

print("Expected on-hand inventory: ", basestock.I(r))
Expected on-hand inventory:  5.10347867979

The total cost is given by

h = 15.
b = 25.
Y = lambda r: h*basestock.I(r) + b*basestock.B(r)

for r in range(5, 20):
   print(r, Y(r), X.cdf(r), b/(b+h))
5 104.39955586 0.067085962879 0.625
6 84.6052126955 0.130141420882 0.625
7 68.4140385596 0.220220646602 0.625
8 56.7268257096 0.332819678751 0.625
9 50.0440142885 0.457929714472 0.625
10 48.3656042962 0.583039750193 0.625
11 51.2366501483 0.696776146303 0.625
12 57.8989092041 0.791556476395 0.625
13 67.4774861089 0.864464422619 0.625
14 79.1391471915 0.916541527065 0.625
15 92.1895310593 0.951259596696 0.625
16 106.107866668 0.972958390215 0.625
17 120.536762124 0.985722386403 0.625
18 135.24930194 0.992813495396 0.625
19 150.111128261 0.996545658024 0.625

As demand is discrete, there may not be an \(r\) such that the critical fractile formula can be solved. We see in the table above that

\[G(10) < 0.625 < G(11),\]

and that \(Y(10)\) is minimal. Hence, it is necessary to compute both \(Y(10)\) and \(Y(11)\), i.e., at \(r\) and \(r+1\) such that \(G(r) \leq b/(b+h) \leq G(r+1)\) to find out the minimal cost. If we apply this formula right away, we get:

print("Optimal reorder level r*: ", basestock.rstar())
Optimal reorder level r*:  11.0

We get the same results, as expected.

Normal demand.

When the expected replenishmenttime demand \(\theta\) is not too small, the normal distribution can serve as a good approximation of the Poisson distribution. Lets see what we get if we apply this approximation.

theta = 10
sigma = np.sqrt(theta)
X = stats.norm(loc = theta, scale = sigma)
basestock = BaseStock(l=1, h = 15., b = 25.,  X = X)

r = 14

print("Service level: ", basestock.S(r))
print("Expected backlog: ", basestock.B(r))
print("Expected on-hand inventory: ", basestock.I(r))
print("Expected cost: ", basestock.Y(r))
print("Optimal reorder level r*: ", basestock.rstar())
print("Minimal reorder level for a fillrate of 0.9: ",
basestock.meetServiceLevel(0.9))
Service level:  0.897048394634
Expected backlog:  0.07682904031971743
Expected on-hand inventory:  5.07682904032
Expected cost:  78.0731616128
Optimal reorder level r*:  11.0076261423
Minimal reorder level for a fillrate of 0.9:  14.0526218861

An Efficient Frontier

Factory Physics for Managers discusses the concept of efficient frontier for inventory systems. To idea is to compare the performance of a company’s inventory position to the best possible performance. Suppose here that a certain item is controlled by a base stock policy. We can vary the reorder point \(r\) to control the fillrate and the average inventory level. Assuming that we want to achieve a certain fillrate, we like to know the (minimal) average on-hand inventory that is required for this fillrate. Once we know this, we can compare the company’s actual fillrate and the inventory in stock to see how things can be improved. For further discussing see Factory Physisc, or Factory Physics for managers.

To make an efficient frontier for the Superior Appliances example of Section 2.4.2, we vary \(r\) from 0 to some number \(N\), and compute the service level \(S(r)\) and average on-hand inventory \(I(r)\). Then we plot \(I(r)\) as a function of \(S(r)\), so that by choosing a certain service level it is easy to see what on-hand inventory is associated with that service level.

import matplotlib.pylab as plt
import seaborn; seaborn.set()

theta = 10
l = 1 # leadtime
h = 15.
b = 25.
X = stats.poisson(theta)
basestock = BaseStock(l=l, h=h, b=b, X=X)

N = 30
S = [0]*N # Store the service levels
I = [0]*N # Store the inventory levels
for r in range(N):
    S[r] = basestock.S(r)
    I[r] = basestock.I(r)

plt.plot(S, I)
plt.xlim((-0.01, 1.01))
plt.ylim(ymin=-1)
plt.xlabel("Service level")
plt.ylabel("Average on-hand inventory")
plt.show()
../../_images/baseStock_figure13_1.png

Here we can see that to achieve a service level of 90% we need an average on-hand inventory of about 5 items. If we are not prepared to invest in this amount of inventory, we have to accept that service will decrease. On the other hand, assume that the current service level is 80% and that the (realized) average on-hand inventory is 10 items. Then it must be possible to realize a much higher service level, for instance, by improving the inventory control policy used.