An Inventory with exponential leadtimes under the Base-Stock Policy

All code in one file

We analyze an inventory process subjected to a base-stock policy. Orders arrive according to a Poisson process with rate \(\lambda\) at an inventory. If orders cannot be met from on hand inventory, they are backlogged up to some level \(R\). Orders that arrive when the backlog level is equal to \(R\) are rejected. When the inventory position hits the reorder level \(S-1\) from above, a replenishment order is issued. The replenishment lead time is assumed to be exponentially distributed with expectation \(L\).

It is well known that the inventory level is equal to the number of customers \(p_n\) in the \(M/G/\infty\) queue. Hence,

\[p_n = e^{-\rho} \frac{\rho^n}{n!}\]

where \(\rho = \lambda L\).

We implement this model as a subclass of the from the class MarkovChain, see Markov Chains.

Standard imports

from collections import namedtuple
import numpy as np
from scipy.stats import poisson
from scipy.misc import comb
from markovChain import MarkovChain

A state is only characterized by its inventory level \(i\).

State = namedtuple('State', "i")

Observe that when \(k\) replenishments are outstanding the return rate must be \(k/L\). Also, when the inventory level is \(i\), the number of replenishments must be \(k = S-i\).

class BaseStock(MarkovChain):
    def __init__(self, labda, L, S, R):
        super(BaseStock, self).__init__()
        self.labda = labda
        self.L= L
        self.S = S
        self.R = R # (minus) the max number of backorders
        self.initstate = State(S)

    def mu(self, i):
        return 1.*(self.S-i)/self.L

    def transitionStates(self, state):
        i = state.i
        tostates = {}

        if i < self.S:
            tostates[State(i+1)] =

        # down
        if i > self.R:
            tostates[State(i-1)] = self.labda

        return  tostates

    def doAll(self):
def base():
    kwds = {'labda': 1., 'L': 3., 'S': 4, 'R': -15}
    base = BaseStock(**kwds)
    print "Poisson: ", poisson(kwds['labda'] * kwds['L']).pmf(0)
    rho = kwds['labda'] * kwds['L']

The last few lines are intended as a check.

Pweave does not run this code if we call base() within the if-clause. Hence we quote it out.

#if __name__ == "__main__":
State(i=-15) 4.75690876824e-10
State(i=-14) 3.01270923359e-09
State(i=-13) 1.80762578309e-08
State(i=-12) 1.02432141821e-07
State(i=-11) 5.46304831998e-07
State(i=-10) 2.73152453842e-06
State(i=-9) 1.2747116279e-05
State(i=-8) 5.52375115303e-05
State(i=-7) 0.00022095007674
State(i=-6) 0.00081015039365
State(i=-5) 0.0027005016864
State(i=-4) 0.00810150618188
State(i=-3) 0.0216040194789
State(i=-2) 0.0504093857697
State(i=-1) 0.100818785511
State(i=0) 0.16803133247
State(i=1) 0.224041807674
State(i=2) 0.224041838721
State(i=3) 0.149361246512
State(i=4) 0.0497870890702
Poisson:  0.0497870683679