# 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 = {}

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

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

return  tostates

def doAll(self):
self.generateStateSpace(self.initstate)
self.generateTransitionMatrix(True)
self.computePi()

def base():
kwds = {'labda': 1., 'L': 3., 'S': 4, 'R': -15}
base = BaseStock(**kwds)
base.doAll()
base.printPi()
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__":
base()

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