# 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.determineDistributions()
self.computeGamma()

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

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

# g just about right
bisect("bisectRight.dat", g = 2.3)

# g too large
bisect("bisectTooLarge.dat", g=2.5)


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


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