Microsoft Word - TOC_R.doc
HUNGARIAN JOURNAL
OF INDUSTRIAL CHEMISTRY
VESZPRÉM
Vol. 37(2) pp. 89-94 (2009)
PROBABILITY AND EXPECTED TIME OF EMPTYING OF INTERMEDIATE
STORAGES UNDER STOCHASTIC OPERATION
É. ORBÁN-MIHÁLYKÓ1 , B. G. LAKATOS2, CS. MIHÁLYKÓ1, L. LUCZ1
1University of Pannonia, Department of Mathematics, H-8200 Veszprém, Egyetem u. 10. HUNGARY
E-mail: orbane@almos.uni-pannon.hu
2University of Pannonia, Department of Process Engineering, H-8200 Veszprém, Egyetem u. 10. HUNGARY
A mathematical model for the operation of an intermediate storage under stochastic operational conditions is presented.
The material input of the storage occurs at randomly, both in terms of time and amount of material, while the output is of
constant volumetric rate. The main practical question is the required initial amount of material in the storage allowing no
emptying of storage. In order to determine that we define a two variable function. We develop an equation satisfied by
this function and present solution for the case of Erlang(n) distributed inter-arrival time. We provide numerical examples
and apply this function for solving the original problem.
Keywords: Process system, random operation conditions, intermediate storage, reliability, integro-differential equation.
Introduction
Intermediate storage is frequently used in processing
systems. In chemical engineering, the material produced
by batch operational units is often collected in a storage
system from which it is withdrawn for further processing
in the subsequent plant units. The reasons of applying
storage may be the different operational characteristics
of the input and output subsystems, or necessity of
sparing the processed material in case of failures during
the operation. It seems to be an important question to
determine the necessary size of the storage in order to
avoid the overflow and the necessary initial amount of
material not to run out the material stored. This is to
ensure the work without failure of the continuously
operated subsystem supplied by material from the storage
system.
Investigating the operation of intermediate storage,
researchers realised that the operation is rather stochastic
than deterministic [2]. General stochastic model was set
up for batch-batch systems [3], and a method was given
for determining the size and the initial amount of
material required for operation of the system without
failures at a given reliability level.
In this paper we investigate such a model in which
the input of the storage system is formed by a pro-
cessing system of a large set of batch operational units
while the output of the storage occurs continuously, i.e.
the material from the storage is withdrawn at a constant
rate. In previous work [4], we investigated this problem
assuming Poisson input process. We have set up the
equations modelling the behaviour of the system which
were solved analytically in some special cases, and we
have given a method for solving the sizing and initial
amount problem based on these solutions or computer
simulation.
In this paper we deal with the problem of initial
amount of material under more general conditions. We
determine the probability of running out of material in
the storage as a function of the initial amount of material.
On the other hand it is interesting to also know the
expected time of this event. In order to handle these
problems together, we introduce a general function
which at the value zero provides the probability of
running out of material, and its derivative provides the
expected time of emptying as well. Applying the tools
of probability theory we set up an integro-differential
equation satisfied by this function. We show some
qualitative properties of the solution and in some cases
we determine analytical solutions for it. Furthermore,
we present numerical examples of solutions which are
applied for solving the practical problem.
The model: basic assumptions and notations
Consider a processing system consisting of a batch
and a continuous subsystem that are connected by means
of an intermediate storage system, shown schematically
in Fig. 1. The material to be processed is filled into
the storage by the input sub-systems. The filling
happens randomly, both in terms of time and amount
of material. The investigation starts at t0 = 0 and the
time intervals between the consecutive filling points
are described by random variables ...3,2,1,kt
Therefore, the ith filling is at time ∑
=
i
k
kt
1
.
90
Figure 1: Intermediate storage connecting the batch and continuous subsystems
Let the random variables ti be assumed to be
nonnegative, independent, identically distributed random
variables, hence the input process is a renewal process.
If the random variables are exponentially distributed,
then the input process is a Poisson process. In a general
case let their common distribution function be denoted by
F(t), while the density function be denoted by f(t). Let
μf be their finite expectation.
Let N(t) denote the (random) number of fillings in
the time interval [0, t]. In the ith filling event, the
amount of material filled into the storage is random
as well, and it is described by the nonnegative
random variable Yi, i = 1, 2, …. These are supposed
to be independent of each other. Their distribution
function is denoted by G(y), the density function is
g(y), while the finite expectation is μg. Furthermore,
we assume that the input process N(t) and Yi are
independent as well.
The output process is a deterministic process with
constant output rate c. The main goal is to determine
the initial amount of material necessary at a given
reliability level, which means that there occurs no
running out of the material with a given probability.
For this purpose it is worth investigating the change
of material in storage as a function of time.
The change of material is the difference of the
filled and the withdrawn material. If we denote the
initial amount of material by x, which is the amount
of material at t = 0, and we sum up the change to the
initial amount of material x, we get the function
describing the amount of material in the storage in
terms of time. This is a random function, which means
that considering it in terms of time provides a stochastic
process. In order not to run out of material it is necessary
to satisfy the following inequality
0)(
)(
1
≥−+= ∑
=
ctYxtV
tN
i
i (1)
for any value of t ≥ 0 where V(t) denotes the amount
of the material in the storage at time t. As V(t) is a
stochastic process, the inequality (1) holds for some
realization and does not hold for other ones, i.e. it
holds only with a given probability (reliability).
Let us introduce a function describing reliability,
i.e. the probability of not emptying as a function of
the initial amount of material.
If x is the initial amount of material then we define
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
⎭
⎬
⎫
⎩
⎨
⎧
≤∀−+≤= ∑
=
)(
1
0:0)(
tN
i
i ttctYxPxR .
The probability of running out of the material is
given by
ψ(x) = 1 – R(x),
and the time of emptying is expressed as
}{⎩
⎨
⎧
<≥<≥
≥≥∞
=
0)(:0existsthereif,0)(:0inf
0anyfor0)(if,
)(
tVttVt
ttV
xTV
Furthermore let the function φ(x,δ) be defined as
follows. For 0 ≤ x and 0 ≤ δ we have
)1(),( )(
)(
∞<
−= xT
xT
V
VeEx δδφ .
This function is, in essence, the Laplace transform of
the probability density function of the finite time of
running out of the material. It can be easily shown that
φ(x,δ) = ψ(x),
and
)1(
),(
)(
)(
0
∞<
−
=
=
∂
∂
− xT
xT
V
VeE
x δ
δδ
δφ
,
hence the function φ(x,δ) is a suitable tool for
determining the reliability and the expectation as well.
Consequently, we will deal with the function φ(x,δ).
We mention that similar function was introduced in
insurance mathematics to investigate the ruin probability
in [1] called Gerber-Shiu discounted penalty function.
Equations for the function φ(x,δ)
First we note that the function φ(x,δ) is monotonic
decreasing in both variables. It is continuous and the
inequality 0 ≤ φ(x,δ) ≤ 1 holds. Furthermore it can be
proved that:
Theorem 1: For any x ≥ 0 and δ ≥ 0 function φ(x,δ)
satisfies the following integro-differential equation:
+−+= ∫∫
∞
− dtdyygtfctyxex
c
z
t )()(),(),(
0 0
δφδφ δ
∫
∞
−
+
c
x
c
x
dttfe )(
δ
. (2)
Batch
units
Continuous
subsystem
Intermediate
storage
Input Output
1
K
q(t)=c
Y1i, t1i
YKi, tKi
Yi, ti
91
Theorem 2: Equation (2) has a unique solution for δ > 0
in the set of bounded functions. This solution tends to
zero in exponential order. Namely φ(x,δ) ≤ Ke–rx for
suitable constants 0 < K and 0 < r, depending on δ.
In the case of δ = 0 the bounded solution of Equation
(2) is not unique, as the Heaviside function, i.e. constant
of value one for all x ≥ 0 is a solution to Equation (2)
and, if c
f
g >
μ
μ
and ∫
∞
−≤
x
rxeKdttf 2)( , r > 0, then there
exists such a solution as well for which φ(x,0) ≤ K1e
–rx
with a suitable constant K1 > 0.
Remark: Condition c
f
g >
μ
μ
means that during a unit
time interval the expectation of the amount of filled
material is more than the amount of withdrawn material.
In the opposite case it can be proven that the solution of
the processing problem is
φ(x,0) ≡ 1, x ≥ 0.
Theorem 2 ensures that in the case c
f
g >
μ
μ
for any
0 < α < 1 there exists such an initial amount of material
x for which the material in the intermediate storage will
not run out with probability 1 – α hence the processing
problem can be solved theoretically. The solution will
be provided by taking the inverse function of
φ(x,0) = ψ(x) at point y = α. Consequently we want to
determine the solution of Equation (2).
Theorem 3: Let the random variables be the sum of n
independent, exponentially distributed random variables
with parameter λ, that is Erlang(n) distributed. Now, one
can prove that Equation (2) can be transformed into the
following form:
∫
∑
∞
−
=
+⎟
⎠
⎞
⎜
⎝
⎛
=
=⎟⎟
⎠
⎞
⎜⎜
⎝
⎛
⎟
⎠
⎞
⎜
⎝
⎛ +
∂
∂
0
0
)(),(
),(
dyygyx
c
i
n
cx
x
n
inn
i
i
i
δφ
λ
λδδφ
(3)
with initial conditions
i
x
i
i
cx
x
⎟
⎠
⎞
⎜
⎝
⎛
−=
∂
∂
=
δδφ
0
),(
1,...,1,0 −= ni (4)
If δ > 0, then the function
∑ ∑
=
−
=
−=
m
i
n
j
xkj
ij
i
iexcx
1
1
0
)()(),( δδδφ (5)
is the unique solution of Equation (3) where ki(δ) is the
root of the equation
∫
∞
− =⎟
⎠
⎞
⎜
⎝
⎛
−⎟
⎠
⎞
⎜
⎝
⎛
−
+
0
0 )( dyyge
c
k
c
ky
nn λλδ
, (6)
which is the characteristic equation of Equation (3). The
multiplicity of ki(δ) is ni, and ∑
=
=
m
i
i nn
1
. The constant
values cij(δ) can be uniquely determined from (4) by
solving a system of linear equations. If the limit
)(lim)0(
0
δ
δ ii
kk
+→
= (i=1, …, m)
exists, then
∑∑
=
−
=
−=
m
i
n
j
xkj
ij
i
iexcx
1
1
0
)0()0()(ψ
and
)(lim)0(
0
δ
δ
ijij cc
+→
= .
Corollary: If the assumptions of Theorem 3 are valid and
the roots of characteristic equation are of multiplicity 1,
then
∑
=
−=
n
i
xk
i
iecx
1
)()(),( δδδφ
and
∑
=
−=
n
i
xk
i
iecx
1
)0()0()(ψ .
Now
xk
n
i
iiixTV
i
V
ecxkcxTE )0(
1
''
)( ))0()0()0(()1)((
−
=
∞< ∑ −−=
where
1
0
)0(
1
'
)0()(
)0(
)0(
−∞
−
−
⎟
⎠
⎞
⎜
⎝
⎛
−−⎟
⎠
⎞
⎜
⎝
⎛
⎟
⎠
⎞
⎜
⎝
⎛
−
=
=
∫
n
i
yk
n
n
i
i
k
c
ndyyyge
c
k
cc
n
k
i
λλ
λ
Especially, if n = 1, that is the inter-arrival times are
exponentially distributed, then
xkex )(1),( δδφ −= (7)
xkex )0(1)( −=ψ (8)
xk
yk
TV e
dyyygec
x
xTE
i
V
)0(
0
)0(
1
)(
)1)(( −
∞
−
∞<
∫−
=
λ
. (9)
If Yk is exponentially distributed, that is
0,
1
)(
1
≥=
−
yeyg
y
g
gμ
μ
,
then
92
x
c
ccc
ggg
ex 2
4
)()( 2
),(
μ
δ
δλ
μ
δλ
μ
δφ
+−−+−−−
−
= (10)
x
c gex
)
1
(
)( μ
λ
ψ
−−
= (11)
x
c
g
g
TV
g
V
e
cc
x
xTE
)
1
(
)(
)1)(( μ
λ
λμ
λμ −−
∞< −
= . (12)
In the case when g(y) is a Dirac-delta function at y = 1,
that is Yk ≡ 1, then an equation corresponding to
Equation(2) can be set up by similar argumentation
and its solution expressed as
xkex )(1),( δδφ −= (13)
xkex )0(1)( −=ψ (14)
xkTV ekc
x
xTE
V
)0(
1
1
))0(1(
)1)(( −∞< −+
=
λ
, (15)
where k1(δ) is the unique root with positive real part of
the characteristic equation
0=−−
+ −ke
c
k
c
λλδ
(16)
for δ ≥ 0.
We note that as ψ(x) is exponential function in both
cases, consequently its inverse function can be
analytically determined easily.
Numerical results
We discuss some examples of functions presented in the
previous sections. Fig. 2 shows the function φ(x,δ)
determined by expression (13) in the case of exponentially
distributed time interval tk with expected value μf = 0.5
(λ = 2), constant amount of filled material Yk = 1 and
constant withdrawing rate c = 1. Here, x is the initial
amount of material and δ is a parameter. The exponent
k1(δ) was computed numerically from formula (16).
Figure 2: The form of function φ(x,δ) depending on the
initial amount of material x and parameter δ
Fig. 3 demonstrates φ(x,0) = ψ(x) using the previously
given parameters. The exponent in expression (14) is
actually 1.5936.
0 0.5 1 1.5 2 2.5 3
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figure 3: Diagram of function ψ(x) describing the
probability of running out of the material
depending on the initial amount of material
The inverse function of ψ(x), i.e. ψ–1(α), provides the
necessary initial amount of material in the function of
variable α, where α provides the probability of running
out of material in storage. Function ψ–1(α) is shown in
Fig. 4 depending on the parameter α.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
1
2
3
4
5
6
7
8
Figure 4: The initial amount of material as a function of
the probability of emptying the storage
The expectation of the time of emptying given by
formula (15) as a function of the initial amount of
material can be seen in Fig. 5.
The function φ(x,δ) given by expression (10) is
shown in Fig. 6 in the case of exponentially distributed
consecutive filling times with parameter μf = 0.5 and
exponentially distributed filled amount of material with
parameter μg = 1.
φ(x,δ)
x
δ
μf = 0.5
Yk ≡ 1
c = 1
ψ(x)
x
α
ψ–1(α)
μf = 0.5
Yk ≡ 1
c = 1
μf = 0.5
Yk ≡ 1
c = 1
93
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
Figure 5: The expectation of the time of emptying the
storage depending on the initial amount of material
Figure 6: Function φ(x,δ) in the case of Poisson filling
process and exponentially distributed filled
material into storage
The functions expressing the probability of emptying
can be seen in Fig. 7 using different values for the
expectation of the filled material.
Figure 7: The probability of emptying storage ψ(x) as a
function of initial amount of material
In Fig. 8 we represent the function φ(x,δ) in case of
Erlang(2) distributed consecutive input times with
λ = 0.3, that is
3.0
2
)( =itE , and lognormal distri-bution
of the filled amount of material. Parameter m equals 2,
the standard deviation of the normal distribution is σ = 1,
hence the expectation of the lognormal distribution is
18.12)( 2/
2
== +σmi eYE . The solution of the integro-
differential equation (3) is given by expression (5). The
exponents were determined numerically from the
equation (6). The integral for the case of lognormal
distribution of the filled amount of material, i.e. the
Laplace transform of the lognormal density function can
not be computed analytically. It was determined
numerically by using Gauss-Laguerre quadrature
formulas based on polynomial of degree 20. For this
case, the bivariate function φ(x,δ) is presented in Fig. 8.
Figure 8: Function φ(x,δ) in case of Erlang (2)
distributed inter-arrival times and
lognormal distribution of filled amount of material
Finally, let us see an example in details how to
determine the initial amount of material according to the
reliability level 1 – α = 0.95, i.e. we would like to
determine how much initial amount of material is needed
if we want to assure the operation of the system without
running out of material with probability 0.95. It means
that the initial amount of material in the storage x is to
be determined corresponding to R(x) = 1 – α = 0.95.
Let us suppose that the distribution of the consecutive
filling times is Erlang(2) distribution with parameters
λ = 2.1, c = 1 and Yi ≡ 1.
First we determine the roots of the equation ψ(x) = α.
As ψ(x) = φ(x,0), we deal with the function φ(x,δ), which
is of the form (5). Now g(y) is a Dirac-delta function at
y = 1, hence Equation (6) takes the form
0
22
=⎟
⎠
⎞
⎜
⎝
⎛
−⎟
⎠
⎞
⎜
⎝
⎛
−
+ −kye
c
k
c
λλδ
.
This equation has two positive roots, one lies in the
interval (0, λ + δ), and the another one is in the interval
(λ + δ, ∞) with multiplicities equal 1. These roots were
)1)(( ∞<VTV xTE
μf = 0.5
μg = 1
c = 1
δ
φ(x,δ)
x
μg = 1
μg = 2
μg = 3
μf = 0.5
c = 1
φ(x, δ)
x
λ = 0.3
n = 2
m = 2
σ = 1
c = 2
x
δ
x
μf = 0.5
Yk ≡ 1
c = 1
ψ(x)
94
determined numerically by the Newton method. The
solution is of form
)(
2
)(
1
21 )()(),( δδ δδδφ kk ececx −− +=
where the constants c1(δ) and c2(δ) are expressed from
the initial conditions for i = 0 and i = 1, that is from the
system of linear equations c1(δ) + c2(δ) = 1 and
c
ckck
δ
δδδδ −=−− )()()()( 2211 . The function φ(x,δ)
is shown in Fig. 9.
Figure 9: Function φ(x,δ) in the case of Erlang (2)
distributed inter-arrival times and constant amount of
filled material
The numerical computations provided the values
k1(0) = 0.1968, k2(0) = 2.6564, while the constants are
c1(0) = 1.0800, c2(0) = –0.0800, hence finally
xx eex 6564.21968.0 0800.00800.1)( −− −=ψ .
The function ψ(x) = 1 – R(x) = φ(x,0) can be seen in
Fig. 10. This function is shown together with the
function ψ(x) equals 0.05 for x = 15. We numerically
determined the roots of equation ψ(x) = 0.05, which is
x = 15.61 and ψ(x) = 0.01 which holds for x = 23.7945.
Therefore, if we want to ensure operation of the
continuous processing system without emptying of the
intermediate storage with probability 0.95 we have to
start to operate the system with initial amount of
material equal to 15.6154 units.
0 2 4 6 8 10 12 14 16 18 20
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figure 10: Function ψ(x) in case of Erlang (2) distributed
inter-arrival times and constant amount of material
However, in spite of having this amount of material the
storage may still become empty with probability 0.05.
The expected time of this event is 5398.12)1( =
∞<VTV
TE
units.
REFERENCES
1. GERBER H. U.,SHIU E. S. W.: On the time value of
ruin, North American Actuarial Journal 2 (1998),
48–78.
2. LEE E. S., REKLAITIS G. V.: Intermediate sto-rage
and operation of batch processes under batch failure.
Comp. Chem. Eng., 13 (1989), 411–498.
3. ODI T. O., KARIMI I. A.: A general stochastic model
for intermediate storage in non-continuous processes.
Chem. Eng. Sci., 45 (1990), 3533–3549.
4. ORBÁN-MIHÁLYKÓ É., LAKATOS G. B.: Inter-
mediate storage in batch/continuous processing sys-
tems under stochastic operation conditions. Comp.
Chem. Eng., 28 (2004), 2493–2508.
x
φ(x,δ)
x
λ = 2.1
n = 2
c = 1
Yi ≡ 1
x
ψ(x)
δ
λ = 2.1
n = 2
c = 1
Yi ≡ 1