On Optimal Policies for Production-Inventory Systems with Compound Poisson Demand and Setup Costs

This page illustrates how to find, by means of the concept of g-revised cost, optimal policies for production/inventory systems of single items with setups. The algorithm is based on methods described in [B1]; the details are expounded there.

The page consists of two parts. The first part concentrated on the python code. The second part shows how to apply the method to a number of cases.

Code and Experiments

All code in one file

Standard imports

import sys
import numpy as np
import pylab as pl

from scipy.stats import  uniform, expon
from scipy.optimize import brentq

The Algorithm

The class StochasticEPQ contains the main algorithms. Most important are the methods computeGamma and computeU.

Recall that \(\gamma(x) = h(x) + \lambda \int_0^\infty \gamma(x-y) G(y)\). To obtain good starting values, approximate \(\gamma(u)\) by \(\gamma(z)\) for all \(u\leq z\) and for \(z\) sufficiently small. Then the integral equation reduces at \(z\) to \(\gamma(z) = h(z) + \rho \gamma(z)\) since \(\lambda \int_0^\infty G(y) dy = \lambda E D = \rho\). Thus, as starting values for \(\gamma\) we take \(\gamma(z) = h(z)/(1-\rho)\).

class StochasticEPQ(object):
    def __init__(self, **args): #r, labda, rv, K, h, tol = 1e-3):
        self.__dict__.update(args)
#         self.r = r
#         self.labda = labda
#         self.rv = rv
#         self.h = h
#         self.K = K
#         self.tol = tol

        self.load = self.labda*self.rv.stats('m')
        self.checkLoad()
        self.determineDistributions()
        self.computeGamma()

    def __repr__(self):
        return "s = %4.3f; S = %4.3f; g =  %4.3f"%(self.s, self.S, self.g)

    def checkLoad(self):
        if self.r - self.load < 0.01:
            print >> sys.stderr, "load is too high."
            print >> sys.stderr, "r = ", self.r
            print >> sys.stderr, "labda ES) = ", self.load
            quit()

    def determineDistributions(self):
        leftIndex = np.searchsorted(grid, 0)
        end = self.rv.ppf(1. - delta/10.)
        rightIndex =  np.searchsorted(grid, end)
        support = grid[leftIndex:rightIndex]
        self.G = self.rv.sf( support )
        self.f = self.rv.pdf( support ) * delta

    def computeGamma(self):
        f1 = self.r-self.labda*self.G[0]*delta
        self.gamma0 = self.h(grid)/(self.r - self.load)
        l = len(self.G) # find the support of G
        for j in range(l+1,len(grid)):
            res = self.h(grid[j])
            res += self.labda*delta*np.dot(self.gamma0[j-1:j-l:-1],self.G[1:])
            res /= f1
            self.gamma0[j] = res

    def computeGammag(self):
        # subtract the revision cost g from gamma0
        self.gammag = self.gamma0 - self.g/(self.r - self.load)

    def finds(self):  # find small s, not S
        if self.gammag[0] < -np.finfo(float).eps:
            # g is too large, hence s is not in the grid
            return -1, None
        if self.gammag.min()>-np.finfo(float).eps: # g is too small
            return None, None

        for i, x in enumerate(self.gammag):
            if x < 0: #find s
                return i, grid[i]
    @property
    def s(self):
        si, s = self.finds()
        return s

    def findS(self):  # find large S, not s
        i = self.u.argmin()
        return i, grid[i]

    def isSOk(self):
        # it may occur that at the end of the grid, gamma is still
        # negative. Then the computation of S is not necessarily
        # reliable.
        return self.gammag[-1] > np.finfo(float).eps # to be on the safe side

    @property
    def S(self):
        Si, S = self.findS()
        return S

    def computeU(self, g):
        self.g = g
        self.computeGammag()

        si, s = self.finds()
        if si == None: # g is too small
            return self.K
        if si < 0:
            # g is too large. Hence, we can safely return -1
            return -1.
        if not self.isSOk():
            print >> sys.stderr, "Right end point of grid is too small."
            print >> sys.stderr, "gamma[-1] < 0-"
            quit()
        l = len(self.f) # l is the support of the distribution
        if si<l:
            print >> sys.stderr, "Left starting point of grid is too large."
            print >> sys.stderr, "The support of the distribution is larger than
"
            print >> sys.stderr, "s - grid starting point."
            quit()

        u = np.zeros_like(grid)
        u[:si] = self.K # u(x) = K on x <= s

        # integrate on x > s
        gamma = self.gammag
        l = len(self.f)
        for j in range(si,len(grid)):
            res = self.r*self.gammag[j]/self.labda
            res += np.dot(u[j-1:j-l:-1],self.f[1:])
            res /= 1.-self.f[0]
            u[j] = np.minimum(self.K, res)
        self.u = u
        return u.min()

    def computeOptimalg(self):
        # first locate a  proper lower bound gl and upper bound gu  for g
        gl = 0.; gu = 1.
        while self.computeU(gu) >= 0.+ self.tol:
            gl = gu; gu *= 2.
        # now search for g between the lower bound and the upper bound
        self.g = brentq(self.computeU, gl, gu, xtol = self.tol)

    def saveDataForGnuplot(self, filename):
        res = np.transpose((grid, self.h(grid), self.gammag, self.u))
        np.savetxt(filename, res)

    def computePolicyCost(self, s, S):
        # compute cost for a given s, S policy
        # not yet checked.
        f1 = self.r-self.labda*self.G[0]*delta
        time = ones_like(grid)/(self.r - self.load)
        l = len(self.G) # find the support of G
        for j in range(l+1,len(grid)):
            res = 1.
            res += self.labda*delta*np.dot(time[j-1:j-l:-1],self.G[1:])/f1
            time[j] = res

        si = searchsorted(grid, s)
        Si = searchsorted(grid, S)

        u = zeros_like(grid)
        # integrate on x > s
        l = len(self.f)
        for j in range(si,len(grid)):
            res = self.r*time[j]/self.labda
            res += dot(u[j-1:j-l:-1],self.f[1:])
            res /= 1.-self.f[0]
            u[j] =  res
        print u[Si], u[si]
        time = u[Si]


        u = zeros_like(grid)
        l = len(self.f)
        for j in range(si,len(grid)):
            res = self.r*self.gamma0[j]/self.labda
            res += dot(u[j-1:j-l:-1],self.f[1:])
            res /= 1.-self.f[0]
            u[j] = res
        print u[Si], u[si], u[si+1]
        cost= u[Si] + self.K
        print cost/time

Cases

Bisecting the revision cost

Here we show how the bisection on \(g\) works. When \(g=2.1\) it is too small, so that the minimum of the value function \(u\) is still positive. When \(g=2.5\), it is too large as the minimum of \(u\) is negative. The bisected value \(g=2.3\) is just about right.

The demand is uniform on \([0, 1]\); holding costs are linear.

def bisect(filename, g):
    global delta
    global grid

    delta = 0.005
    grid = np.arange(-20.,10.,delta)
    def h(x): return (-3.*np.minimum(0,x) + 1.*np.maximum(0,x))

    demand = uniform(loc = 0., scale = 1.)

    c = StochasticEPQ(r=1, K=5., labda=3./2., rv=demand, h=h, tol = 1e-3)

    # g too small
    c.computeU(g )
    c.saveDataForGnuplot(filename)
    pl.title("$s$ = %3.2f, $S$ = %3.2f, $g$ = %3.2f"%(c.s, c.S, c.g))
    pl.plot(grid, c.gammag, label = 'gamma')
    pl.plot(grid, c.u, label = 'u')
    pl.plot(grid, h(grid)-c.g, label = 'h')
    pl.grid(True)
    pl.legend()
    pl.xlim((-5., 10));
    pl.ylim((-2, 5.5))
    pl.show()

# g too small
bisect("bisectTooSmall.dat", g = 2.1 )
../../_images/stochasticEPQ_figure3_1.png
# g just about right
bisect("bisectRight.dat", g = 2.3)
../../_images/stochasticEPQ_figure4_1.png
# g too large
bisect("bisectTooLarge.dat", g=2.5)
../../_images/stochasticEPQ_figure5_1.png

The case of Graves and Keilson with with exponential demand

Graves and Keilson, in MOR 6(2), 1981, 246–262, consider a production/inventory problem with exponentially distributed demand.

def gravesKeilson():
    global delta
    global grid

    delta = 0.005
    grid = np.arange(-20.,30.,delta)
    def h(x): return (-2.*np.minimum(0,x) + np.maximum(0,x))

    demand = expon(scale=1.)

    c = StochasticEPQ(r=1, K=25., labda=0.9, rv=demand, h=h, tol = 1e-3)

    # Let's take g = 4.68 as in GK
    c.computeU(11.44)
    print c,
    print "theta = S - s = %4.3f"%(c.S - c.s)

    # the optimal result:
    c.computeOptimalg()
    print c

gravesKeilson()
s = 1.320; S = 15.345; g =  11.440 theta = S - s = 14.025
s = 4.925; S = 10.375; g =  9.158

We do not get exactly the same results as Graves and Keilson, which we suspect is due to the fact that Graves and Keilson restrict \(s, S\) to the integers, while we leave them free.

Comparison to the D-Policy of Tijms

[tijms03:_first_cours_stoch_model], Section 8.2, computes the optimal switching level D for the D-policy. Lets check that we get the same results.

def dPolicy():
    # see Tijms
    global delta
    global grid

    delta = 0.0005
    grid = np.arange(-20.,1.,delta)
    def h(x): return (-1.*np.minimum(0,x) + 10000.*np.maximum(0,x))

    def fitErlang(mu,scv):
        from scipy.stats import erlang
        alpha = 1./scv
        labda = alpha/mu
        rv = erlang(alpha, scale = 1./labda)
        return rv

    demand = fitErlang(1, 0.5)

    c = StochasticEPQ(r=1, K=25., labda=0.5, rv=demand, h=h, tol = 1e-3)

    c.computeOptimalg()
    print c

dPolicy()
s = -1.288; S = -0.000; g =  1.204

Tijms gets D=2.847, which is nearly the same as our result.

A dent in the value function

def dent():
    global delta
    global grid

    delta = 0.005
    grid = np.arange(-20.,10.,delta)
    def h(x): return (-4.*np.minimum(0,x) + np.maximum(0,x))

    demand = uniform(loc=1.02, scale = 0.1)

    c = StochasticEPQ(r=1, K=5., labda=0.2, rv=demand, h=h, tol = 1e-3)
    c.computeU(g = 1.7)
    c.saveDataForGnuplot("dent.dat")
    print c
    pl.title("$s$ = %3.2f, $S$ = %3.2f, $g$ = %3.2f"%(c.s, c.S, c.g))
    pl.plot(grid, c.gammag, label = 'gamma')
    pl.plot(grid, c.u, label = 'u')
    pl.plot(grid, h(grid)-c.g, label = 'h')
    pl.grid(True)
    pl.legend()
    pl.xlim((-10,10));
    pl.ylim((-10, 20))
    pl.show()

dent()
s = -0.280; S = 1.090; g =  1.700
../../_images/stochasticEPQ_figure8_1.png

General holding costs

Now the inventory cost function is much more general function, and the demand is exponentially distributed. See E.L. Porteus, Stochastic Inventory Theory, Stanford University Press, 2002, Exercise 9.5, for the explanation behind the cost function.

def notsS():
    global delta
    global grid

    delta = 0.005
    grid = np.arange(-40.,15.,delta)
    #def h(x): return (-4.*(x<=0)+ (x>0))*np.ceil(x)
    def h(x): return (-4.*np.minimum(x,0))+ 4.*(x<=0) + \
        (x>0)*(np.ceil(x)-0.25*(x-np.floor(x)))

    demand = expon(scale = 2.)

    c = StochasticEPQ(r=1, K=5., labda=0.2, rv=demand, h=h, tol = 1e-3)
    c.computeOptimalg()
    print c
    c.saveDataForGnuplot("notss.dat")
    pl.title("$s$ = %3.2f, $S$ = %3.2f, $g$ = %3.2f"%(c.s, c.S, c.g))
    pl.plot(grid, c.gammag, label = 'gamma')
    pl.plot(grid, c.u, label = 'u')
    pl.plot(grid, h(grid), label = 'h')
    pl.grid(True)
    pl.legend()
    pl.xlim((-10,10));
    pl.ylim((-10, 20))
    pl.show()

notsS()
s = 1.710; S = 2.995; g =  5.447
../../_images/stochasticEPQ_figure9_1.png