Investigation of operation of intermediate storages applying probability density functions satisfying linear differential equation

In this paper operation of a batch/continuous processing system connected by an intermediate storage is considered. The filling process is supposed to be random and the output process is deterministic, consequently the process in the storage is a stochastic process. We investigate the problem of determination of necessary initial amount of material to avoid emptying of the storage. We define a function which is able to handle together the probability and expected time to shortage. We set up an integral equation for it, and in a special case of density functions of inter-arrival times we transform to an integro-differential equation. We provide analytical formula for the solution. We compare the analytical solution to the results arising from MonteCarlo simulations. In some cases we use only the form of the analytical solution but coefficients are determined by parameter fitting. Finally we use the computed functions to determine the necessary initial amount of material to a given reliability level.


Introduction
Intermediate storage is often used in chemical industry, pharmaceutical factories, food logistic or in case of environmental investments.The material produced by a factory is collected in a storage which operates as a buffer.The material is withdrawn from this buffer and usually is processed by further industrial operations.The reason for buffering is the different characteristics of producing and processing or provide appropriate spare amount in case of faults.One of the main goals of investigations is to determine the necessary initial amount of material.
The operation of intermediate storages is rather stochastic than deterministic [1][2][3].One can realize stochastic features in producing both in time and in amount of material.In this paper we investigate a stochastic model, when the filling process is a batch process and the withdrawing process is a continuous one.Such processes were previously investigated in [4,5].
The characteristics of the random behaviour of the producing process are difficult to determine.The distribution function or the density function of the random variables are impossible to determine exactly on the basis of data, consequently it is worth applying such distribution functions which are close to the distribution functions in real world.Such type of distribution functions are those which have probability density function satisfying a linear differential equation with constant coefficients.Their linear combinations are probability density functions as well and this set is dense in the set of all density functions having support [0,∞).In this paper we turn our attention to the case of these probability density functions.
We introduce an auxiliary function which is able to handle both the probability and the time to shortage of material in a unified way and we set up an integral equation for it.After transforming this equation into an integro-differential equation we present its analytical solution.We compare the results for probability and time of emptying based on these analytical solutions and the results coming from Monte-Carlo simulations.The form of the analytical solution provides opportunities for determining the constant values by parameter fitting.We present that parameter fitting provides good approximation of the analytical results.Finally we apply our results for solving the problem of finding the necessary initial amount of material to a given reliability level.

The investigated model, problem formulation, notation and assumptions
The material filled into the intermediate storage is produced by a production system.The fillings happen at random time points and the amount of filled material is random as well.The process begins at time 0 and the i th filling is at time i k=1 t k , i = 1, 2, . . . .The inter-arrival time between the (i-1) th and i th filling is t i .We suppose that t i , i = 1, 2, . . .are independent, identically distributed nonnegative continuous random variables with distribution function F(t), density function f(t), their common expected value is µ f .Let N(t) be the number of fillings in the interval [0, t], that is The amount of material filled into the buffer during the i th filling (i = 1, 2, . . . ) are described by independent, identically distributed nonnegative continuous random variables with distribution function G(y), density function g(y), expected value µ g .The filled amount of material is denoted by Y i ,i = 1, 2, . . . .We suppose that the filling process N(t) and filled material Y i , i = 1, 2, . . .are independent as well.The withdrawal process is continuous and withdrawal happens at constant rate c.This means that the material withdrawn from the buffer from time 0 to time t equals ct.Our main question is how many initial amount of material should be in the intermediate storage at time 0 in order to avoid the shortage (i.e.emptying the storage) in a large time interval.
To solve this problem we investigate the change of material in the intermediate storage in the function of time.If the initial amount of material is x (0≤x), then the amount of material being in the storage is ( In order to avoid the shortage of material, the inequality 0≤V(t) has to be satisfied in the interval 0≤t≤T, if the process is investigated in the time interval [0,T], or for any 0≤t, if the process is analyzed in unbounded time interval.V(t) is random for any fixed value of t, therefore V(t) is a stochastic process as a function of time.If the inequality 0≤V(t) does not hold for any value of t, we can define the time to shortage which is the 'first' time point when V(t)<0 holds, that is ) Let us define the probability of the shortage as well.Let ψ be the function describing reliability, that is The probability of shortage equals 1 − Ψ(x).
One realization of the stochastic process can be seen in Fig. 1.Inter-arrival times and filled material are exponentially distributed random variables with parameters µ f = 2.5 and µ g =4, respectively.The initial amount of material is 5 and the constant rate of withdrawal equals 2. The decreasing linear lines show the withdrawal of material and the jumps represent fillings.Shortage of material means that V(t) is below zero.In the case presented in Fig. 1 it happens first at about 8.9 unit of time.
Actually let us define the function φ(x, δ)for 0≤ x and 0≤ δ as follows: where E denotes the expectation of the random variable being in parenthesis, 1 A is the characteristic random variable of event A and x is the initial amount of material.This function is often used in different risk models in insurance mathematics and was introduced by Gerber and Shiu in 1998 [6].In economic context δ is a discounting factor.Returning to our model, one can realize that φ(x, δ) is the Laplace transform of the finite shortage time's density function and δ denotes the variable of the Laplace transform.
One can easily check, that If the function φ(x, δ) is given, then on the basis of ( 6) and ( 7) the probability of the shortage and the expected time to shortage can be computed.Consequently, the function φ(x, δ) is appropriate for handling the probability and expected time to shortage in a unified way.Therefore we will investigate the function φ(x, δ) in the next sections.

Integral equation for the investigated function
First we mention that investigation and determination of the function φ(x, δ) requires many complicated and sophisticated mathematical tools, therefore we will publish the details in a mathematical journal in the immediate future.Actually we state only the necessary theorems.

78
Éva Orbán-Mihálykó / Csaba Mihálykó Figure 1 A realization of the stochastic process in the function of time Comparing it to the results arising from Monte-Carlo simulation we get very good coincidences.The plotted points are simulated results while the continuous function is given by (35).Simulations were executed to T=3000 N=10000 times.
This equation has a unique bounded and continuous solution for any 0< δ.Moreover this solution is exponentially bounded, namely φ(x, δ) ≤ e −δ x c .
Therefore φ(x, δ) → 0, if x → ∞ for any 0< δ.If δ=0 then the solution of Eq.( 8) is not unique.If condition is satisfied then it can be proved that φ(x, 0) → 0 supposing x → ∞.This boundary condition assures that the solution of ( 8) is the function given by (??).As φ(0, 0) = 1 holds, for any 0 < α < 1 there exists a value x for which This value x is the necessary initial amount of material to the reliability level 1 − α.Consequently, determination of function φ(x, δ)is important for solving our problem.We note that condition (10) expresses that the material filled into the storage is more than the material withdrawn from the buffer, in average.
We note that if condition is satisfied, then, applying probabilistic argumentation, one can prove that φ(x, 0) ≡ 1.

Case of density functions satisfying linear differential equation
Explicit solution for the integral equation ( 8) can not be provided generally hence we restrict ourselves to a special set of density functions f (t).In order to be widespread applicable, this set of density functions should be dense in the set of all density functions continuous on [0,∞).This means that the continuous density functions can be well approximated by the elements of this set.We investigate the case when the density function of the inter-arrival times satisfies a linear differential equation with constant coefficients subject to some general initial conditions.This set of density functions is dense in the set of continuous density functions [8] and is investigated in risk models as well [9].Furthermore, it contains such famous distributions as Coxian distributions, K n -type distributions, Erlang(n) distributions, combination of exponential distributions.For this type of random variables we transform the integral equation ( 8) into an integro-differential equation and we provide its analytical solution.It can be proved that the solution of the integro-differential equation set up by the help of the approximate density function is close to the solution belonging to the real density function (which was approximated).
We also suppose that g(y) is a continuous function.
First we summarize what is known about a density function that satisfies a linear differential equation with constant coefficients.Let be a polynomial with coefficients a n = 1, a 0 0. Assume that the density function f (t) satisfies the linear differential equation subject to initial conditions From the theory of the differential equations we know the following: if the number of different roots of the polynomial p(z)is m, the different roots are denoted by ρ i , and the multiplicity of the root ρ i is s i , then the solution of the linear differential equation is where b i j are appropriate real numbers.f (t) is a density function, therefore in the case Re(ρ i ) ≥ 0 b i j = 0 hold for all j = 1, 2, , . . ., k i .After some computations we can see that holds if and only if Investigation of operation of intermediate storages applying probability density functions 79 Theorem 2 [7] Let p(z) = n i=0 a i z i , with a n = 1, a 0 0, and suppose that the density function f (t) satisfies Eq. ( 13) subject to initial conditions (14).Then, with the notation q 0 (z) ≡ 0, q i (z) function φ(x, δ) satisfies the following integro-differential equation subject to appropriate initial conditions, where I denotes the identity operator.The initial conditions can be determined by substituting 0 into the equations obtained after taking derivatives of Eq.( 8), but their determination can be difficult.We note that if A k = 0 for k = 0, 1, 2, . . ., n − 2 and a 0 = A n−1 0, then the initial conditions are Actually we present the solution of Eq. (20) in the following theorem: Theorem 3 [7] Suppose that the density function f (t) satisfies Eq. ( 14) with a n = 1, a 0 0 subject to the initial conditions (15).If the polynomial p has roots with nonnegative real parts as well, it can be simplified to a polynomial of lower degree and f satisfies the linear differential equation defined by the reduced polynomial as well.Consequently, we can assume that all the roots of the polynomial p(z) have negative real parts.Then where the coefficients c i j (δ) are appropriate real numbers and k i (δ) are the roots of with multiplicity n i (δ) for i = 1, 2, . . ., r, where Especially, if all the roots of Eq.( 23) has multiplicity 1, then Consequently, and where

Application of analytical results and comparison to the results arising from Monte-Carlo simulation
In this section we present how to determine analytically the probability of shortage on the basis of the previously presented equations.In order to be able to compare the exact solution and simulation results, furthermore for the sake of simplicity of computations first we suppose that t i , (i=1,2,. . . ) are Erlang(2) distributed with parameter λ, that is f (t) = λ 2 te −λt , and Y i are exponentially distributed with parameter λ y = 1/µ g .We note that µ f = 2/λ.Now, with p(z) =(z+λ) 2 , f(t) satisfies Eq.( 14) subject to the initial conditions A 0 =0, A 1 = λ 2 .After some computation we get Eq.( 23) looks like The roots of Eq. ( 30) can be determined analytically.If δ=0 then two of them is negative and one is zero.The negative ones are Both these roots have multiplicity 1, therefore from (25) The coefficients c 1 (0) and c 2 (0) can be easily determined applying (21) for j = 0, 1 and substituting x = 0: and Per. Pol.Chem.Eng.
If we would like to determine the expectation of the time to shortage we need the derivative of c i (δ) and k i (δ) at δ = 0, i=1,2, as well.k i (0), i=1,2 can be determined from (27).c i (0) can be calculated as follows: take the derivatives of (21) with respect to δ and substitute x = 0and δ = 0.After that we get a linear system of equation for c i (0) (i=1,2) to solve.In Fig. 2 we present the exact probability of shortage in the case of Erlang(2) distributed inter-arrival times with parameter λ = 1 (µ f = 2) and exponentially distributed amount of material with λ y = 0.4 (µ g = 2.5), furthermore c=1.The negative roots of Eq. ( 23) were k 1 (0) = −0.1367and k 2 (0) = −1.4633(see (31)).From (32), the function describing probability of emptying is Comparing it to the results arising from Monte-Carlo simulation we get very good coincidences.The plotted points are simulated results while the continuous function is given by (35).Simulations were executed to T =3000 N=10000 times.To determine the expected time of emptying we need the derivatives.Applying (27) and solving the mentioned system of equations we get The function describing the expected time of emptying is presented in Fig. 3.The parameters were the same as previously.
The simulation was repeated N=10000 times.One can realize that the results coming from simulations and arising from the analytical formula are close to each others.Determination of coefficients k i (0) is based on Eq. ( 23) but c i (0) depend on the explicit form of initial conditions.In such cases when the initial conditions are difficult to compute, we To determine the expected time of emptying we need the derivatives.Applying (27 solving the mentioned system of equations we get , , - can fit a curve of form (22) to simulated results of probabilities with known values k i (0) and unknown values c i (0).Applying the least square method for determining c f 1 (0) and c f 2 (0) we get c f 1 (0)=1.1016and c f 2 (0)= -0.1016, and the most difference between the fitted curve φ f (x, 0) and the exact one φ(x, 0) is 1.2 10 −3 .The difference is presented in Fig. 4. Parameters of the simulation were the same as previously.
After determining c f 1 (0) and c f 2 (0) by least square method, we computed the derivatives of them.We compared the exact form of E(T V (x)1 T V (x)<∞ and the formula using fitted coefficients E f (T V (x)1 T V (x)<∞ ).We realized that their difference is sufficiently small.The maximal difference was 0.018, which shows that the determination of unknown parameters by least square method applying the function type given by ( 26) is an appropriate method for approximation.Consequently, analytical result is very useful because it serves the form the fitted function and we have to apply least square method to determine the coefficients.
curve and the exact one The difference is presented in Parameters of the simulation were the same as previously.
After determining c f 1 (0) and c f 2 (0) by least square method, we computed the der of them.We compared the exact form of and the formula usin coefficients ( ).We realized that their difference is sufficiently The maximal difference was 0.018, which shows that the determination of un parameters by least square method applying the function type given by (26 appropriate method for approximation.Consequently, analytical result is very because it serves the form the fitted function and we have to apply least square me determine the constant values.Finally we present an example when the initial conditions (21) for j = 1 are not given explicitly.We compute the roots of Eq. ( 23) exactly then we determine the coefficients by fitting parameters by least square methods.Let Two roots of (38) can be given analytically at δ = 0 as follows: 2 , (39) They are negative numbers if (10) holds.We note that there is a third root as well, namely zero.Now (25) can be applied with n=2 and φ y, coefficients and can be determined easily by least square method.Actually, coefficients c 1 (0) and c 2 (0) can be determined easily by least square method.Fig. 5 shows the simulation results and the fitted curve in the case of parameters λ 1 = 1/6, λ 2 = 2/3, µ g = 5,c=5.The roots of (38) are k 1 (0) = −0.0804and k 2 (0) = −0.5530.The simulation was executed N=1000 times to T=1000.The fitted coefficients are c f 1 (0) = 0.7485 and c f 2 (0) = 0.2518.The fitted function looks φ f (x, 0) = 0.7485e −0.0804x + 0.2518e −0.5530x , (41 and it is very close to the points computed by Monte-Carlo simulation. Returning to the problem how much initial amount of material is needed to the reliability level 1 − α, according to (11) we have to solve the equation 0.7485e -0.0804 x + 0.2518e −0.5530x = α.It can be executed numerically.In case of α = 0.1 we get x ≈ 25.0, in case of α = 0.05 we get x ≈ 33.65, and in case of α = 0.01 we get x ≈ 53.62.

Summary
Operation of a batch/continuous processing system connected by an intermediate storage is investigated.The input process is supposed to be a random process and the output occurs at constant volumetric rate.The probability and the expected time to shortage are handled by the help of the Laplace transform of the density function of finite shortage time.Integral equation is set up for it, transformed and solved in special cases.The analytical solution is compared to the solution coming from simulation supporting the applicability of Monte-Carlo simulation.The form of the analytical solution gives possibility for parameter fitting and finding the solution numerically.Finally we use the computed functions to solve the physical problem, i.e. to determine the necessary initial amount of material to a given reliability level.

Figure 2 Fig. 1 .
Figure 2 Probability of shortage of material in the function of initial amount of material (-analytical result, * simulated results) x

Figure 2 Figure 3 Fig. 2 .
Figure 2 Probability of shortage of material in the function of initial amount of material

Figure 3 Fig. 3 .
Figure 3 Expectation of the time to shortage in the function of the initial amount material (-analytical result, * simulated result) Fig. 3. Expectation of the time to shortage in the function of the initial amount of material (-analytical result, * simulated result)

Figure 4 φ
Figure 4 Differences between the exact probabilities ) 0 , x ( φ and values of the function ) 0 , x ( f φ

5 Fig. 5 .
Fig. 5. Comparison of the fitted curve and the simulation results in the case of linear combination of exponential probability density functions (-fitted function given by (37), * simulated result)