1. INTRODUCTION
The simple economic order quantity (EOQ) model is the most fundamental
of all inventory models (see, e.g., [15]). It lucidly describes the trade-off between
the constant setup cost and the variable holding cost. Its formula still
serves as an effective approximation for much more complicated EOQ-type
models. Although the simple EOE (SEOQ) is a deterministic model with
constant demand rate, most of its ramifications include added factors of
randomness. For these stochastic EOQ models, three types of approximation
are generally used:
- In jump models (see, e.g., [8,13,14]) discrete amounts of input and output enter
and leave the buffer one by one or in small batches; the content level
processes generated in this case are step functions.
- In heavy-traffic models [3,4], both the gross
input and the gross output are assumed to be very large (approaching
infinity) over any time interval, with the ratio of input rate and output
rate being close to unity and their difference being almost constant.
Generally, the buffer content processes include Brownian components.
- Fluid models [5,6] are characterized by large amounts of inflow
and outflow over any time interval, with the ratio of input rate and
output rate not equal to unity; the stochastic fluctuations of the input
and the output streams are replaced by their local drifts, which might
also depend on external factors changing randomly over time.
In this article we introduce a fluid EOQ (FEOQ) model with a two-state
random environment that is modeled as a continuous-time Markov chain
alternating between a good and a bad state. The net
change rate of the system is assumed to be constantly equal to a
in the good state and to b in the bad state. It can be assumed
that a > b, which is intuitive since the sales during
good periods are higher than those during bad ones, but this is not
necessary for the analysis. Good and bad periods follow each other
according to an alternating renewal process. In this article we restrict
attention to the Markovian case of independent and exponentially
distributed random variables with parameters λ and μ,
respectively.
The two-state model reflects situations in which the demand rate for a
certain commodity undergoes periodically recurring changes. For example,
the demand for beef went down dramatically after the occurrence of BSE
(mad cow disease) cases and then after some time of alertness, it swung
back to its normal rate. Similar phenomena connected to diseases or health
consciousness (e.g., popular diet plans) can be observed for poultry,
pork, and other foodstuffs. The sales of furs are periodically influenced
by campaigns of animal rights groups, and those of mineral water and
canned food are influenced by terror alerts. Indeed, demand rates for many
goods go up and down between different levels due to fashion or other
recurring external effects. Models including a multistate Markovian
environment can be suitable for such situations; the two-state case
presented in this article could serve as a first approximation. We denote
by X = {X(t) : t ≥ 0} the content
level process of the inventory operating under the alternating good and
bad states. X decreases linearly at rates a or
b depending on the state of the environment. Once the system
becomes empty, a replenishment order, with a negligible lead time, of a
certain controllable size is placed. Note that as demand is assumed to
arrive in infinitesimal portions at rate a or b and the
lead times are zero, there is no backlogging. Let I =
{I(t) : t ≥ 0} be the status
process of the environment; that is, I(t) = 1 (0)
if at time t the state of the environment is good (bad). The
two-dimensional process (X,I) is Markovian and the
content level process X alone is regenerative. There are several
options for defining regeneration points for it. The cycle structure we
will use is given by the periods between replenishments occurring in good
periods. Let (X(0),I(0)) =
(qg,1) (i.e., X starts at level
qg in the good state). We define T
to be the time of the first replenishment taking place in a good period;
then [0,T) is the first cycle. One expects X to
proceed swiftly toward zero during good periods, whereas during bad
periods, it will probably decrease more sluggishly. If this is the case,
the marginal revenue is higher during good periods than during bad ones.
Whenever the content level process hits zero at good (bad) periods, the
controller places an order of size qg
(qb), where qg
and qb are controllable parameters. Note that
the successive order sizes form a random sequence of values in
{qg,qb},
determined in accordance with the random environment. However, in contrast
to random yield models [10,11] in which the random order sizes are
determined externally (“by nature”), in this model they are
selected optimally.
In the monograph [17] a large
variety of inventory models is presented in detail. The stochastic models
are based on point processes for the demand arrivals in random
environments (in the book called “world-driven”). The fluid
systems in [17] are deterministic EOQ
models with the classical extensions such as planned backorders, limited
capacity, quantity discounts, and imperfect quality. In the deterministic
setting, time-varying demands are considered also, however without
multiple-order quantities. The stochastic FEOQ model with more than one
order quantity expounded in this article seems to be new.
The controller's objective is to maximize the long-run average
revenue by selecting qg and
qb so as to properly balance between the
setup cost K, due each time an order is placed, and the
proportional holding cost h per unit time and per unit of stored
items. Let p and c be the sale price and the purchase
price of one unit, respectively. The long-run average holding cost is
hEX, where X is a random variable having the
steady-state distribution of X. (Note that X(t)
→ X in distribution, by the limit theorem for regenerative
processes [1, p.170], and that 0 ≤
X(t) ≤
max[qg,qb]
so that also EX(t) → EX.) Let N
be the number of replenishment orders in one cycle (say, the first one)
that are issued while the environment is in the bad state. Then the number
of orders in this cycle is N + 1 and it is easily seen that
N has a geometric distribution. The expected setup cost in a
cycle is KE(N + 1); the mean cycle length is, of course,
ET. The expected total sold output during a cycle is
qg + qb
EN. Combining all of these functionals, it follows from the
renewal reward theorem [16] that the
long-run average profit is given by
Of course, the expected values EN, EX, and
ET are all functions of qg and
qb.
Because the long-run average income generated by the sales depends
only on the proportions of times in good and bad periods, it is
independent of the order sizes, implying that the maximization of
R(qg,qb)
is equivalent to the minimization of the cost functional
It is easy to see that the long-run average profit and the long-run
average cost are related by
An intuitive conjecture is that if a > b, the
optimal (long-run average cost minimizing) replenishment levels
qg* and qb*
satisfy qg* >
qb*. The numerical analysis in Section 4
confirms this expectation.
In order to maximize (1.1) (or minimize (1.2)) we first need to
compute the functionals EN, EX, and ET. To this
end, we use an extension of the so-called level crossing technique (see,
e.g., [2,4,6]). In our context, this approach turns out to
be more intricate than in most other applications because the exact
calculation of the upcrossing rates requires some refined arguments.
The level crossing theory (LCT) was introduced in [7] for regenerative dam processes of the
GI/G/1 type and generalized in [9] to stationary dam processes. The general
approach makes it possible to construct Khintchine–Pollaczek
formulas for dam processes by equating the long-run average number of
upcrossings and downcrossings of an arbitrary level. It is based on the
idea that the long-run average of upcrossings (or downcrossings) of level
x gives the value of the steady-state density of the dam process
at x multiplied by the local release rule of the dam. In
[7,9] and
the applications following these two pioneer studies, attention is
restricted to the case that the release rule of the dam is a deterministic
function. In the model considered here, the demand rate is not state
dependent but stochastically changing according to the Markov environment,
which makes a rigorous analysis more complicated.
The article is organized as follows. In Section 2 the dynamics of the
FEOQ model are described in a formal manner. In Section 3 all relevant
functionals are derived in closed form from a steady-state analysis of the
model. Based on these explicit results, Section 4 provides numerical
examples in which we consider, in particular, the sensitivity of the
optimal ordering policies with respect to the model parameters.
2. MODEL DYNAMICS
Let
U1,V1,U2,V2,…
be the lengths of the alternating time periods in which the content
decreases at rates a and b, respectively. For the
analysis, we do not need the assumption a > b. We
assume that (U1,U2,…) and
(V1,V2,…) are two
independent sequences of independent and identically distributed (i.i.d.)
random variables, distributed according to exp(λ) and exp(μ),
respectively. The two threshold values qg and
qb are positive numbers satisfying
qg ≥ qb. Let
Y(t) be the accumulated output at time
t ≥ 0 and set Y = {Y(t) :
t ≥ 0}. Let
for n ≥ 1 and S0 =
S0′ = 0. In terms of these sequences,
Y(t) is defined by
The sample paths of Y increase strictly and continuously to
infinity. The status process I = {I(t)
: t ≥ 0} is given by
The sample paths of Y increase strictly and continuously to
infinity. Now we define recursively a sequence of level hitting times
τ1 < τ2 < ···
[searr] ∞ and the clearing process W =
{W(t) : t ≥ 0}. Let τ1 =
inf{t > 0|Y(t) =
qg} and W(t) =
Y(t) for t ∈ [0,τ1).
If τn and W(t), t
∈ [0,τn), have already been defined,
let
and set
By this construction, W(t) is defined for all
t ≥ 0. The content level process X =
{X(t) : t ≥ 0} is given by
We will show that W, and thus X, has an absolutely
continuous stationary distribution. The associated densities are denoted
by fX and fW,
respectively. Clearly,
The processes W and X are regenerative. Typical
realizations are depicted in Figures 1a and 1b.
W starts to increase at time 0 from W(0) = 0 at slope
a. Slopes alternate between a and b. Each time
level qg is reached by W. The
process has a negative jump leading to 0 if the current slope is
a or to qg −
qb if it is b. Cycles are defined to
be the intervals between jumps to zero. Let N be the number of
jumps to qg −
qb before the first return to zero. Then
T = τN+1 is that value
τn for which
I(τn) = 1 and
I(τj) = 0 for j < n.
[0,T) is the first cycle of W. The time intervals
[0,τ1),[τ1,τ2),…,[τn,τn+1)
are called subcycles.
Typical realizations of W and X.
Let θ1(x) be the probability that level
x is upcrossed by Y while growing at rate a;
similarly, let θ0(x) be the probability that
x is upcrossed by {Y(t −
U1) − aU1 : t ≥
U1} while the process grows at rate a. Note
that {Y(t − U1) −
aU1 : t ≥ U1} is the
accumulated output process starting in the bad state (i.e., growing at
rate b). These probabilities will appear throughout our
derivations.
The following facts are consequences of the strong Markov property of
W. For any n ≥ 1, the subcycle lengths
τ1,τ2 −
τ1,…,τn+1 −
τn are conditionally independent given that
N = n. Moreover, τ2 −
τ1,…,τn −
τn−1 are conditionally i.i.d. given that
N = n, provided n > 1. The number N
+ 1 of subcycles has a modified geometric distribution with parameter
θ0(qb): We have
P(N = 0) =
θ1(qg) and
3. STEADY-STATE ANALYSIS
In this section we derive ET and
For both, we need explicit formulas for θ1(x)
and θ0(x).
Lemma 1:
and
where λa = λ/a
and μb = μ/b.
Proof. Consider the process Z =
{Z(t) : t ≥ 0} obtained from Y by
replacing its linear pieces of slope b by independent
exp(μb)-distributed upward jumps. Formally,
let
where N(t) = max{n ≥
0|Sn ≤ t}. Z is
the sum of a compound Poisson process and the linear function at
and thus has the exponent
It is well known that the process M = {M(t)
: t ≥ 0} defined by
is a martingale (see [12]). Fix
x > 0 and define the stopping time
Applying the martingale stopping theorem to M and
Tx, we find that
where the third step of (3.5) follows from the lack-of-memory property
of the jump size distribution: Given the event
{Z(Tx) > x}, the random
variable Z(Tx) − x is
exp(μb) distributed for any x so that the
conditional Laplace transform of
Z(Tx) is
e−αxμb
/(μb + α).
The right-hand side of (3.3) has the two roots 0 and
If we can set α = α* in (3.5), we get
However, the event {Z(Tx) =
x} is equivalent to the event that Y upcrosses level
x during a good period. This means that
Solving for θ1(x) in (3.7), we obtain
(3.1).
There is a difficulty with inserting the root α* in (3.5): α*
< −μb is negative whereas the Laplace
transforms of Z(1) and of the
exp(μb)-distribution are only defined for
nonnegative arguments α and the defining integrals are not finite for
α < −μb. This problem can be solved
as follows. Being the logarithm of a Laplace transform, φ(α) is
only defined for
.
However, φ(α) can be analytically extended to
by identity (3.3), taking the right-hand side of (3.3) as definition of
φ. Next, the integral
is easily seen to be an analytic function of α for all
;
note that 0 ≤ Z(s) < x for s
∈ [0,Tx). Now let
Then G is analytic on
,
H is analytic on
,
and, by (3.5), G and H coincide on
.
By the identity theorem for analytic functions, we obtain
G(α) = H(α) for all
.
Since G(α*) = 0, it follows that H(α*) = 0,
i.e., (3.7).
To prove (3.2), let us make the dependence of
θ1(x) on λa and
μb explicit by writing it as
θ1(x,λa,μb).
Then a little reflection shows that
and (3.2) follows from (3.1). █
The following explicit formula for the steady-state density of
W is the main analytical result of this article. It is derived by
means of an intricate use of the LCT. The constant factor
fW(0+) can be determined by the normalizing
condition for the density; using another argument, an explicit formula for
fW(0+) is given in Theorem 1.
Theorem 1: The stationary density
fW of W is given
by
where θ0(·) and
θ1(·) have been computed in Lemma 1,
and fW(0+) can be determined from the
normalizing condition
.
Proof: By the limit theorem for regenerative processes
[1, p.170], the stationary
distribution function FW of W is
given by
where 1A denotes the indicator function of the set
A. Fix 0 < x < qg.
Then 0 < x < qg −
ε for sufficiently small ε > 0. We now use the basic
arguments of the level crossing approach. Consider the integral
which is equal to the amount of time that W spends in the
interval [x,x + ε) during the cycle
[0,T). Clearly,
(where an empty sum is defined to be zero). For x ∈
(0,qg − qb
− ε), the sum on the right-hand side of (3.12) vanishes since
all of the indicator variables in it are equal to zero. For x
∈ [qg −
qb,qg −
ε), each of them is equal to unity during an interval of length
either ε/a or ε/b or between these
numbers, because during each subcycle, the sample path of W is
increasing and runs through [x,x + ε)
linearly at alternating slopes a and b. It follows
that
for all ε ∈ (0,qg −
x). As W(t) is piecewise linear, the
derivative
exists, and (3.13) and a similar estimate for negative ε
together yield
Since N has a modified geometric distribution (see Section
2), we have
Thus, we can use dominated convergence to conclude that
where the third equality follows from (3.11). Hence,
FW is an absolutely continuous distribution
whose density can be computed by taking the derivative of
.
Let us compute
.
If x ∈ (0,qg −
qb), the interval [x,x
+ ε) is crossed exactly once in the cycle [0,T), and
during this crossing, the sample path has slope a with
probability θ1(x) + O(ε) and slope
b with probability 1 − θ1(x) +
O(ε) as ε [nearr ] 0. (The terms O(ε)
are due to the possibility that [x,x + ε) can
also be crossed by a piece of sample path in which the slope changes.)
Hence, the expected amount of time spent in [x,x +
ε) is
Relation (3.19) accounts for the term in square brackets in the
corresponding assertion of (3.8) for x ∈
(0,qg −
qb).
Now let x ∈ [qg
− qb,qg).
In this case we can write
where εY denotes the sojourn time of W in
[x,x + ε) during those subcycles of
[0,T) that start and end in the bad state. Indeed, the
expected amount of time spent in [x,x + ε)
during the first subcycle is the same as in the case x ∈
(0,qg −
qb); that is, it is given by (3.19). The
series in (3.20) gives the contribution of the other subcycles of
[0,T). With probability 1 −
θ1(qg), the first subcycle is
followed by n ≥ 0 consecutive subcycles starting and ending in
the bad state and a final subcycle starting in the bad state and ending in
the good state. The counting variable of the upcrossings of
[x,x + ε) at slope a (i.e., while
being in the good state) during the n subcycles starting and
ending in the bad state is a binomial random variable with success
probability
and the probability that the upcrossing of [x,x
+ ε) during the final subcycle (starting in the bad state and ending
in the good state) occurs while growing at slope a is
The corresponding values for slope b are 1 −
γ(x) + O(ε) and 1 − ν(x) +
O(ε), respectively. The probability of having n
consecutive subcycles leading from bad to bad followed by one leading from
bad to good is (1 −
θ0(qb))nθ0(qb).
These arguments explain (3.20). The series on the right-hand side of
(3.20) can easily be calculated in closed form, leading to the term in
square brackets in assertion (3.8) for x ∈
[qg −
qb,qg).
To complete the analysis, it remains to show that
During a cycle level, x is upcrossed exactly once for every
x ∈ (0,qg −
qb). Since the cycle [0,T)
starts in the good state, the initial slope is a, so that the
slope while crossing [x,x + ε) is a
with probability 1 − O(x) as x → 0.
Thus, by (3.18),
yielding (3.21).
Since fW is a probability density on
(0,qg), the value
fW(0+) can be obtained by the normalizing
condition ∫0qg
fW(x) dx = 1. The proof is
complete. █
An explicit formula for ET and thus
fW(0+) is given in the following theorem.
Theorem 2:
Proof: From (1.1) and (1.3), we can conclude that
(3.22) now follows from EN = (1 −
θ1(qg))/θ0(qb).
█
4. NUMERICAL EXAMPLES AND SENSITIVITY
ANALYSIS
We start with several asymptotic formulas. These results are quite
intuitive and their proofs, while not difficult, are quite tedious in some
cases and therefore omitted.
Lemma 2:
Notice that parts (a), (c), and (g) provide the classical EOQ results
for a system that is always in a good period. However, when a
→ ∞, the expected cycle time is zero (part (e)), whereas when
λ → 0 or μ → ∞, it is positive and finite. When
a → ∞, an order of size
qg /2 is made an infinite number of times
(part (h)). Parts (b) and (g) provide the classical EOQ result of expected
inventory and expected number of orders for a system that is always in a
bad period. However, when μ → 0 or λ → ∞, the
length of the cycle tends to infinity (part (d)).
Next we investigate the effect of parameter changes on the optimal
solution. The following example is typical of many other examples studied.
Consider a system with λ = 0.4/unit of time, μ = 0.6/unit
of time, a = 10/unit of time, and b = 5/unit of
time. We let h vary in {0.5,1,1.5,2,…,5} and K
vary in {5,10,15,20,…,50}. The optimal values of
qg,qb,EX,ET,EN,
and
C(qg,qb)
are given in Table 1, parts a–f,
respectively. The following observations can be made from Table 1:
1. As expected, as K increases, the optimal values
of
qg,qb,EX,ET,
and
C(qg,qb)
increase and the optimal value of EN decreases.
2. As expected, as h increases, the optimal
qg,qb,EX,ET,
and
C(qg,qb)
decrease and the optimal EN increases.
3. The expected number of orders (1 +
EN)/ET decreases when K increases, and it
increases when h increases.
4. When K increases (h increases), both
qg and qb
increase (decrease), as noted in observations 1 and 2, but their
difference qg −
qb increases (decreases). In other words,
qg is more sensitive to changes in K
than qb, and qb
is more sensitive to changes in h than
qg.
Optimal Values as a Function of h and
K
In Table 2, parts a–f, we set
h = $1/(unit × unit of time) and K =
$10/order and let a vary in {12,13,14,…,22} and
b vary in {1,2,3,…,11}. The following observations can be
made:
5. As a increases,
qg,EX, and
C(qg,qb)
increase and ET decreases. Both qb
and EN decrease as a increases.
6. As b increases,
qb,EX, and
C(qg,qb)
increase. Also, qg,ET, and
EN increase as b increases.
7. When a increases (b increases), the
difference between qg and
qb increases (decreases).
8. When either a or b increase, the
expected number of orders increases.
Optimal Values as a Function of a and
b
