Gavish-Graves1980: A Production-Inventory Problem under a continuous review policy

All code in one file

The code below implements the production/inventory problem described by Gavish and Graves, A one-product production/inventory problem under continuous review policy, Operations Research, 28(5), 1980, 1228–1236, and in Gavish and Graves, Production/inventory systems with a stochastic production rate under a continuous review policy, Comput. & Op. Res. Vol 8(3), 169–183, 1981.

Single-unit demands arrive in accordance with a Poisson process with rate \(\lambda\) at a production/inventory system. Demand is met from on-hand inventory; otherwise demand is backlogged. Costs in inventory and backlog level are linear. Production can be switched on at the expense of a setup cost \(K\). When production is on, items are produced at a constant rate of \(\mu\). When production is off, inventory is not replenished. The problem is to determine the optimal two-critical number policy, i.e., a level \(I^{**} = Q_2\) at which to switch off produciton, and a level \(I^* = Q_2 - Q1<I^{**}\) at which to switch on production.

Some preliminaries

This is just boilerplate code.

from __future__ import division

import numpy as np
from  scipy.integrate import quad
import scipy.stats
from scipy.stats import poisson, uniform
from nickysLibrary import memoize

Computing the average costs

Here is the class that computes the expected average cycle costs.

To find the optimal two-critical number policy (the best \(s = I^*\) at which to switch on the server, and \(S = I^{**}\) at which to switch it off) we simply search over all \(Q1\) and \(Q2\) that seem relevant, and print the smallest value found. As the memoization makes the numerical work blazingly fast, there is no real need to use the search procedure advocated in Gavish and Graves. (We simply take the lazy path, and see what happens. Computers now are not what they were in 1980s. Some things do become better over time :-))

class GavishGraves1980(object):
    def __init__(self, md1, h, b, K):
        self.md1 = md1
        self.h= h
        self.b = b
        self.K = K

    def CF(self, Q1, Q2):
        if Q2 >= Q1:
            res = self.h*(2*Q2-Q1+1.)/2.
        else:
            res = self.h*(Q2*(Q2+1))
            res += self.b*(Q1-Q2)*(Q1-Q2-1)
            res /= 2.*Q1
        return res

    def CN(self, Q1, Q2):  # eq 16
        res = (self.h+self.b)*self.md1.theta(Q1, Q2)
        res += self.h*(Q2-self.md1.qB - (Q1-1.)/2.)
        return res

    def C(self, Q1, Q2):   # eq 2
        res = self.CF(Q1,Q2)*self.md1.TF(Q1)
        res += self.CN(Q1,Q2)*self.md1.TN(Q1)
        res += self.K
        res /= self.md1.TN(Q1) + self.md1.TF(Q1)
        return res

    def optimize(self):
        self.Q1opt = self.Q2opt = 20
        self.Min = self.C(Q1=self.Q1opt, Q2=self.Q2opt)
        # a better algorithm should improve the speed considerably
        for Q2 in range(0, 40):
            for Q1 in range(1, 40):
                res = self.C(Q1, Q2)
                if res < self.Min:
                    self.Q1opt, self.Q2opt, self.Min = Q1, Q2, res

    def optimizeSmarter(self):
        raise NotImplementedError

    @property
    def s(self):
        return self.Q2opt - self.Q1opt

    @property
    def S(self):
        return self.Q2opt

Experiments

Comparison to the N-Policy for the M/D/1 queue

When we set Q2=0 the inventory model discussed here should reduce to the M/D/1 queue under the N-policy, see [tijms94:_stoch_model_algor_approac] Section 1.3.5. Let us compare the numerical results to check whether we did not make any errors in the programming or analysis.

We have to implement the deterministic random variable ourselves. With this we can build the M/D/1 queue.

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")

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

B = deterministic(loc = 3)
npol = NPolicy(labda = 0.3, demand = B, K = 500, h = 2.)
md1 = MG1(labda = npol.labda, B = npol.demand)

gg = GavishGraves1980(md1, h = 1., b= npol.h, K = npol.K)

for N in range(1, 10):
    print(N, gg.C(Q1=N,Q2=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

Computing the busy-period steady state probabilities

Gavish and Graves (1981) describe in their Section 5 a procedure to compute the busy period probabilities \(\pi(n)\). I use a method based on a recursion developped in [tijms03:_first_cours_stoch_model]. I show in Computation of the M/G/1 State Probabilities that these methods give the same results. The method of Tijms appears to be faster, and more robust.

Large Holding Costs

When we take the holding cost h to be very large, it must be optimal to set Q2 to zero, since maintaining a positive amount of stock costs too much. Thus, the system must also reduce the M/D/1 queue under an N-policy.

gg = GavishGraves1980(md1, h = 1000., b= 1., K =100.)

for Q2 in range(0, 5):
    print(N, gg.C(Q1 = 1, Q2=Q2))
9 7.95
9 253.156271427
9 636.14033833
9 1134.28954378
9 1726.33658846

Setting the arrival rate to zero

I suspect that \(s\to -1\) and \(S\to 0\) as \(\lambda \to 0\), the intuition being that when hardly any customers arrive, it is best to keep an inventory level of zero, and accept the customers as backorders. In essence, this is a make-to-order environment (keeping stock is simply too expensive when customers arrive very seldom). In the limit \(\lambda =0\), there is simply no reason to keep positive inventory, since customers never arrive for the product.

for labda in np.linspace(0.3, 0.001, 10):
    md1 = MG1(labda = labda, B = deterministic(loc = 3))
    gg = GavishGraves1980(md1, h = 1, b= 10, K =50)
    gg.optimize()
    print(labda, gg.s, gg.S)
0.3 8 13
0.266777777778 2 7
0.233555555556 1 5
0.200333333333 0 4
0.167111111111 -1 3
0.133888888889 -1 2
0.100666666667 -1 2
0.0674444444444 -1 1
0.0342222222222 -1 1
0.001 -1 0

We see indeed that \(s\) becomes \(-1\) and \(S\) decreases to \(0\).

Comparison to the model of Lee and Srinivasan 1987

#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)

B = uniform(loc = 2., scale = 2.)
mg1 = MG1(labda = 0.1, B = B)
gg = GavishGraves1980(mg1, h = 2., b= 20., K =3000.)
print(gg.C(Q1=10,Q2=9))
36.0388253938

Comparing this to Lee-Srinivasan 1987: Continuous review (s,S) policy for production-inventory systems we see a difference. Why?