# Lee-Srinivasan 1987: Continuous review (s,S) policy for production-inventory systems¶

## Introduction¶

Here we implement the model of H.S. Lee and M.M. Srinivasan, The Continuous review (s,S) policy for production-inventory systems with Poisson demands and arbitrary processing times, Tech. Rep. 87-33, Ann Arbor, 1987.

All code in one file

## Code¶

We use memoization for the recursions.

import numpy as np
import scipy
from scipy.stats import uniform, poisson

def memoize(f, cache={}):
def g(*args, **kwargs):
key = ( f, tuple(args) )
if key not in cache:
cache[key] = f(*args)
return cache[key]
return g


The following functions are direct implementations of the formulae in Lee and Srinivasan.

In the recursive definition of the function h in Eq (2.8) we use the function Q which is defined as

$Q(i) = 1 - \sum_{j=0}^i q_j.$

We implement this recursively by setting $$Q(i) = Q(i-1) - q_i$$ if $$i\geq 0$$, and $$Q(-1) := 1$$.

It turns out that the recursion in Lemma 2.3 as proposed by the authors can be considerably simplied, at least for numerical purposes. Let $$\delta(k) = f_{k,k+1}$$, where $$f_{k,k+1}$$ is the expected cost to move from inventory level $$k$$ to $$k+1$$, and more generally $$f_{i,j}$$ is the expected cost to move from level $$i$$ to $$j$$ while production is on. Then from their Equation (2.5)

$\begin{split}\begin{split} \delta_k &= E_k + \sum_{j=1}^\infty q_j f_{k+1-j, k+1}\\ &= E_k + \sum_{j=1}^\infty q_j \sum_{i=0}^{j-1} \delta(k+1-j-i)\\ &= E_k + \sum_{i=0}^\infty \sum_{j=i+1}^{\infty} q_j\delta(k+1-j-i)\\ &= E_k + \sum_{i=0}^\infty \sum_{u=0}^{\infty} q_{u+i+1}\delta(k-u)\\ &= E_k + \sum_{u=0}^{\infty} \delta(k-u) \sum_{i=u+1}^\infty q_{i}\\ &= E_k + \sum_{u=0}^{\infty} \delta(k-u) (1-\sum_{i=0}^u q_{i})\\ &= E_k + \sum_{u=0}^{\infty} \delta(k-u) Q_u \end{split}\end{split}$

where we use the above definition for $$Q_u = Q(u)$$. In the recursion it is necessary to take the term with $$u=0$$ out of the summation leading to

$\delta_k = [1-Q(0)]^{-1}(E_k + \sum_{u=1}^{\infty} \delta(k-u) Q_u).$

The last step is to deal with the summutation to infinity. It is safe to cut it off at 100 as $$Q(100)$$ is very small. Since $$Q$$ is a monotone decreasing function I simply cut off the summation at this point.

class Lee(object):
def __init__(self, labda, demand, ch, cb, K):
self.labda = labda
self.demand = demand
m = demand.stats("m")
rho = labda*m
v = demand.stats("v")
m2 = v+m*m    # the second moment of the demand distribution
self.WT = 1./(1.-rho) * ( labda*m2/(2*(1.-rho)) + m )  # 2.9
self.m = m
self.rho = rho
self.m2 = m2

self.ch = ch # holding cost
self.cb = cb # backlog cost
self.K = K   # switching cost

def TC(self, r,S): # 2.1
return (self.C(r,S) + self.K)/self.L(r,S)

def C(self, r,S):
return self.CN(r,S) + self.CP(r,S)

def L(self, r,S):  # 2.14
return r/( (1.-self.rho) * self.labda)

def CN(self, r,S):  # 2.3
def g(k):
if k>= 0: # 2.2a
return self.ch*k/self.labda
else:     # 2.2b
return -self.cb*k/self.labda
return np.sum(g(k) for k in range(S-r+1, S+1))

def CP(self, r,S): # 2.4
return np.sum(self.delta(k) for k in range(S-r, S))

@memoize
def q(self, i):   # 2.6
res = quad(lambda t: poisson.pmf(i, self.labda*t)*self.demand.pdf(t), 0, np.Inf)
return res[0]

@memoize
def E(self, k):  # 2.7
if k > 0:
return self.E(k-1) + self.h(k)*(self.ch+self.cb) - self.m*self.cb
elif k == 0:
return self.labda*self.m2/2. * self.cb
else:
return self.E(k+1) + self.m*self.cb

@memoize
def h(self, k):  # 2.8
if k <= 0:
return 0.
else:
res = self.h(k-1) + self.Q(k-1)/self.labda
return res

@memoize
def Q(self, i):
return (i==-1) and 1 or self.Q(i-1) - self.q(i)

@memoize
def delta(self, k):
if k <= -2:  # 2.11.b
return self.delta(k+1) + self.m/(1.-self.rho) * self.cb
elif k == -1: # 2.9
return self.WT*self.cb
elif k == 0: # 2.11.b
return  self.delta(-1) - self.m/(1.-self.rho) * self.cb
else:        # 2.5
res = self.E(k) + np.sum(self.delta(k-u)*self.Q(u) for u in range(1, 100))
res /= 1.-self.Q(0)
return res


## Examples¶

### Example of Lee and Srinivasan¶

The long-run average cost of a production cycle are given by TC(r,S). For the values chosen above we obtain the following result.

The parameters correspond to their Example 2

B = uniform(loc = 2., scale = 2.)
lee = Lee(labda = 0.1, demand = B, K = 3000, ch = 2., cb = 20.)

print lee.TC(r=10,S=9)
print lee.TC(r=11,S=10)
print lee.TC(r=20,S=18)

30.2455485443
29.2474742768
28.7594680316


This is precisely equal to their result in Table 2.

### Comparion to an M/G/1 queue controlled by an N-Policy¶

When we take the holding cost ch to be very large, it must be optimal to set S to zero, since maintaining a positive amount of stock costs too much. In this case, the model discussed here should reduce to the M/G/1 queue under the N-policy, see [tijms03:_first_cours_stoch_model] Section 2.6. Let us compare the numerical results to check whether we did not make any errors in the programming or analysis.

In Section 2.6 of [tijms03:_first_cours_stoch_model] the average cost is computed where holding costs are only due for items in queue. [tijms94:_stoch_model_algor_approac] derives in Section 1.3.4 an expression for the average cost when holding costs are due for all items in the system, not just for items in queue.

class NPolicy():
def __init__(self, labda, demand, h, K):
self.labda = labda
self.demand = demand
m = demand.stats("m")
rho = labda*m
v = demand.stats("v")
m2 = v+m*m    # the second moment of the demand distribution
self.WT = 1./(1.-rho) * ( labda*m2/(2*(1.-rho)) + m )  # 2.9
self.m = m
self.rho = rho
self.m2 = m2

self.h = h # holding cost
self.K = K   # switching cost

def longRunCostInQueue(self, N):    # See Tijms 2003, 2.6.4
res = self.labda*(1.-self.rho)*self.K/N
res += self.h* self.labda*self.labda*self.m2/(2.*(1.-self.rho))
res += self.h*(N-1)/2.
return res

def longRunCostInSystem(self, N):  # See Tijms 1994, stochastic models 1.3.5
res = self.labda*(1.-self.rho)*self.K/N
res += self.h*(self.rho + self.labda*self.labda*self.m2/(2.*(1.-self.rho)))
res += self.h*(N-1)/2.
return res

class deterministic_gen(scipy.stats.rv_continuous):
def _cdf(self, x ): return np.where(x<0, 0., 1.)
def _stats(self): return 0., 0., 0., 0.

deterministic = deterministic_gen(name="deterministic")

B = deterministic(loc = 3)

npol = NPolicy(labda = 0.3, demand = B, K = 500, h = 2.)
lee = Lee(labda = npol.labda, demand = B, K = npol.K, ch = 2., cb = npol.h)

for N in range(1, 10):
print N, lee.TC(r=N,S=0), npol.longRunCostInQueue(N), npol.longRunCostInSystem(N)

1 24.9 23.1 24.9
2 18.4 16.6 18.4
3 16.9 15.1 16.9
4 16.65 14.85 16.65
5 16.9 15.1 16.9
6 17.4 15.6 17.4
7 18.0428571429 16.2428571429 18.0428571429
8 18.775 16.975 18.775
9 19.5666666667 17.7666666667 19.5666666667


Clearly, we should compare to the holding cost of the customers in the system, not just in queue.