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.

We follow their notation.

All code in one file

Code

We use memoization for the recursions.

import numpy as np
import scipy
from  scipy.integrate import quad
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.