CC BY
78
On Solution of Renewal Equation in the Weibull-
Gnedenko Model
Vladimir Rusev, Alexander Skorikov
Gubkin Russian State University of Oil and Gas (National Research University) mailto: [email protected], [email protected]@yandex.ru, [email protected]
Abstract
Renewal density of restorable systems and their components which depends on statistical estimates based on real operational data is studied. It is assumed that objects' entire life cycle is described by the Weibull-Gnedenko distribution. Analytical and discrete approaches for the solution to the renewal equation are proposed. New calculation schemes of the renewal density of restorable systems and their components are presented. Equivalence of suggested approaches is illustrated by numerical examples.
Keywords: the Weibull-Gnedenko distribution, reliability theory, renewal density (intensity), numerical methods, collocation knots, moments generating function
1 Introduction
All mechanisms, engineering constructions, operational systems are subjected to the processes of aging, degradation or failures in work. The renewal of normal operation mode possesses doubtless economical and sometimes vital importance. Construction of suitable for applications mathematical models of the renewal processes is thus far an actual challenge in the reliability theory, as existing models involve cumbersome calculations, while analytical solutions are not available in general case. This paper is devoted to the development of analytical and simple discrete schemes for the solution of the renewal density equation. Renewal functions have wide variety of applications in warranty analysis, inventory theory, supplies planning [1]. Examples of processing of real statistical data on refusals of technologically active elements of gas supply systems are considered in [2].
The scheme of a simple renewal process is the following. Let component (or system) failures occur at time moments tt, t2, ... , tn,... and it is assumed that replacement time is negligible relative to the operational time. Then tn represents the operational time until the n-th failure takes place. And it is supposed the time intervals between failures Tn = tn — tn-1 are independent and identically distributed. In this case Tn is the random life time of the n-th item with cumulative distribution function F(t) and probability density function (PDF) f(t), and N(t) is the number of renewals in the time interval (0, t). The renewal function H(t) is the expected value of renewals in that interval H(t) = E(N(t)). The renewal density (intensity) by definition is given by the equality h(t) = H'(t). The fundamental renewal density equation has the form (e.g. [3]):
h(t)=f(t) + f*h(T)f(t-T) dr. (1)
Solution of this equation does not have explicit form, except for some cases when the renewal process is driven by the exponential and the Erlang distributions. In this paper new
analytical and discrete methods of calculating of h(t) are presented for the Weibull-Gnedenko probability density function which depends on two parameters: a, named "scale" parameter and p, called "shape" parameter:
{aPpxP-ie-(ax)P, x > 0
/00 =-jo, x<0. (2)
This distribution was chosen for the study, because it allows to capture all life cycle of systems investigated in the reliability theory, that makes it one of key distributions. By results of many studies the typical curve of the hazard rate ([4]) usually is U-shaped, thus, it contains three main periods of life cycle: initial burn-in, normal operation and degradation. All these periods of the system functioning can be modelled by the Weibull-Gnedenko distribution with different shape parameter ([4], [5]). In particular, the first period corresponds to the Weibull-Gnedenko distribution with parameter p £ (0; 1), the second period — with parameter p « 1 and the third period — with parameter p > 2. It should be noticed that upon transition from the second to the third stage the value of shape parameter jumps from 1 to the value more than 2. And this property thus far remains the opened question for discussion.
For large t it is well-known the asymptotic result for the renewal density function
h(t)~±(t^»),
where n = E(Tn). But the values of h(t) can oscillate (see Fig. 3) about the asymptotic value, thus it is important to have opportunity to calculate values of h(t) more accurately.
W.L. Smith and M.R. Leadbetter [3] developed the method for computation of the renewal function for the Weibull-Gnedenko distribution by using power series expansion of t^, where p is the shape parameter of the Weibull-Gnedenko PDF. However, for p > 1, the numerical computation of this series is limited to the small range of t: 0 < t < 2,5. A. G. Constantine and N.I. Robinson [6] presented estimation method of H(t) (and automatically h(t)) by residue calculations of the Laplace transform of the renewal integral equation to form uniformly convergent series of damped exponential terms. There are also many other approximations of the renewal function explored for the Normal, Gamma, Uniform underlying lifetime distributions developed by L. Cui and M. Xie [7], E. Smeltink and R. Dekker [8], S. Maghsoodloo and D. Helvaci [9].
2 Methods
2.1 The moment problem
In the present paper the analytical solution of the renewal density equation is closely connected with the moment problem or the problem of unique determination of a distribution of a nonnegative random variable by its moments. Consider this problem for basic distributions used in the reliability theory. As for the background the problem of unique determination of a distribution by the sequence of its moments was first investigated by P. L. Chebyshev back in 1874 in connection with research on limit theorems of probability theory.
It can be shown that such distributions as exponential, normal, truncated normal, gamma distribution and the Weibull-Gnedenko distribution (with the shape parameter p > -) are uniquely determined by their moments, and the log-normal distribution, Student's t-distribution and the Pareto distribution can not be determined uniquely. The solution of this problem can be verified by the necessary and sufficient criterion of Krein M. (1944) (or the Krein condition [10]). The distribution with the PDF f(x) is determined uniquely if
r + ro lnf(x2) j
L o dx = ro.
J0 1+X2
2.2 Analytical solution of the integral renewal equation
To solve equation (1), we use the method of moments generating function [11] under the assumption of the two-parameter Weibull-Gnedenko distribution. The Laplace transform (or moments generating function) for the given distribution has the form:
№ = L
m dt = r l
n=0
(St)"
f(t) dt = rn=0(-l)nS-Vn
where vn is the n-th initial moment of a random variable % with PDF f(t):
Vn = c '"AO dt,
which for the Weibull-Gnedenko distribution has the form
and T(s) = f xs 1e Xdx is the Euler gamma function, s £ C. It should be noted that this series is absolutely convergent only for p > 1. Equation (1) in the Laplace transform domain has the form
h(s) = f(s) + h(s)f(s).
Consequently
h(s) =
f(s) 1-f(s)
¿rn=i
n!
nin-
№
(~i)n rf1 n\ (S
Z+=i(-1ï'
nvn„n n!;
1 Vn n
Applying the well-known in calculus technique of dividing infinite series, one can obtain the following expansion:
h(t) = n=0ik-F(k\t), (3)
vi
where F(0\t) = F(t) is the cumulative distribution function and F(k\t) designates the k-th derivative of F(t). The coefficients of the expansion have the following form:
Co = -1,
Ci =
mo,
ck =
1
-mo m1 m2
— r
0 1
-mo m1
(-1) mk-2 (-1)kmk-3
mn
m-,
m--,
where
It was proved that
consequently
mk
lr
_ VlVkn2 (kn2)l
0 m3
1 ...
-mo (-1)kmk_
, k = 0, 1, 2,.
, k = 2, 3,
^ Ck = -1nti?1(-1)k™k < +^
lim ck = 0.
k^rn
Thus, the solution (3) of equation (1) is obtained in terms of the probability moments of the initial distribution of a nonnegative random variable. Moreover, one can notice that the expansion (3) is true not only for the Weibull-Gnedenko distribution (2) under some conditions, not discussed in the present paper.
A good estimate of exact solution is found (Fig. 2) by taking only 7 leading terms of the expansion (3) unlike the series of damped exponential terms in [6] when it is necessary to compute 500 or more coefficients of the given series [4].
The offered solution (3) is applicable only for the values of P > 1, characteristic for the degradation processes [5]. The start of system operation corresponds to the case of 0 < p <1 and its description is sometimes actual. Numerical methods can work with any values of p. So, we proceed to the description of discrete methods.
n
n
0
0
0
1
3.3 The discrete methods
A numerical method which generates a cubic spline approximation of the renewal function by the Galerkin technique for solving the renewal equation was proposed by Z.S. Deligonul and S. Bilgen [12]. The discretizing time method has been used by M.Xie [13] to approximate the renewal equation. Numerical algorithm in papers of M.Xie was based on the definition of the Riemann -Stieltjes integral (RS-method). T.K. Boehme, W.Preuss, V. van der Wall also used the similar method [14]. M. Tortorella [15] presented a paper describing analysis of the method based on quadrature schemes for Stieltjes integrals. Some numerical procedures, when the time scale is discrete, can be found in the books on reliability theory by E.A. Elsayed [1], A. K. S. Jardine, A. H. C. Tsang [16].
In this paper the function h(t) is approximated by step functions or linear functions. The accuracy of this discretization was checked in three different ways. It is obvious that approximation error can be diminished by increasing of the number of collocation knots. All calculations presented below were done by using Wolfram Mathematica software including calculation of integrals (5), (8) without algorithm of numerical quadrature schemes.
The first step in the discretization scheme is the division of the specified time interval [0, t] into n equal-length subintervals by points (collocation knots)
t0 = 0, t1 = t0+ A,... , tn = t0 + nA, where A = t/n is the length of each subinterval.
The method of right knots. Approximate solution un(t) for function h(t) can be found in the form of a linear combination un(t) = YJk=± ukh(t), where uk = h(tk) and the so-called coordinate functions Ik(t) which are equal to zero outside the interval (tk-1, tk] and Ik(t) = 1, tk-1 < t < tk. Thus, the approximate solution is determined from the conditions in right collocation knots. Similar approximation was proposed by A. Brezavscek in [17].
The midpoint method. Let the value of the approximate solution in the k-th segment be the average value (Fig. 1)
uk
uk-i+uk 2 '
uk = h(tk).
Figure 1: The discrete methods
The approximate solution, respectively, has the form un(t) = Yk=i ukIk(t). The Line Spline Finite Element solution. The function h(t) in [tk-approximated by the Lagrange polynomials as follows
tk] can be
i
tk-tk-i tk-tk-i
where uk= h(tk).
Let us consider calculation schemes for the approximate solution of the renewal equation (1). The first method solution is defined by the following recurrent formulas:
(Uo = fo,
Ui = fi/(1 - Fi),
u2 = (f2+ UiF2)/(l - Fi),
u = (f + Y™-1 u
"-n \Jn ' ¿-IJ = 1 n—j
Fj+i)/(l - Fi)
where
fo=f(0),"■ ,fn=f(tn),
(4)
Fi = C f(r)dr, ■■■, Fn = S^ f(r)dr
(5)
The second method gives the next formulas: fuo = ^
= (fl + ^l)/(l-!2),
= {f2 + U-1F2 + Ui0F2 + U-1Fl)/{l-2),
= (fn +~~Fn + Y]—! ^ (Fj+i + Fj)) / (l - 2),
(6)
where fk and Fk,k = 0,1, ...,n are defined by formulas (5)
The recurrent formulas for approximate solutions obtained by linear splines have the form:
(u0 = fo,
Ui = (fi-UoGio)/(1 - Gll),
^2 = (Ï2- UoG2l + Ui(G22 - Gl0))/(1 - Gll)
un = tfn u0Gnn-l + Yj = l un-j(
(Gj+lj+l - Gjj-l))/(1 - Gll),
where
Gkj = S!^ f(r)lj(r)dr,
(7)
(8)
k( )=
tk-1 _ tk-t tk-tk-i A
It should be noted that the first and the third methods were considered earlier [18]. And large computing capacities were required for the solution of the corresponding linear systems. Formulas (4), (7) introduced in the present paper dramatically diminished the computational complexity of the algorithms. So, mentioned algorithms are important if the calculation involves large number of knots.
The second method of discretization, that is, the midpoints method (6) gives the best agreement with analytical solution (3). The third method is time consuming for large values of n > 40.
The discrete methods presented here can be used where the underlying lifetime distribution has the PDF with a singularity at the origin, such as the Weibull - Gnedenko distribution (2) with shape parameter less than unity. It is applied in the modelling of gas supply systems [19].
Fig. 2 illustrates the results of numerical computations of the renewal density performed according to all four methods of this paper. The solutions of equation (1) are presented for the Weibull-Gnedenko distribution with a = 1, P = 2 which coincides the Rayleigh distribution. Red curve for the analytical solution was plotted using seven terms of the expansion (3). There were used twenty collocation knots (n = 20) in calculations by the discrete methods. It is worth noticing
v
u
1
u
2
u
V
V
here that presented methods allow to carry out accurate calculation of the renewal function oscillations unlike the asymptotic formula mentioned above. The corresponding asymptotic value is shown by the straight green line. The curve of the PDF for the Rayleigh distribution is shown with green dotted line.
Figure 2: The Rayleigh renewal density
1.5
£ The Weibull-Gnedenko PDF , J3=4
6 1.0 T3 1 a « Hi oi ■ * ■ i i 1 i i i * — The method of moments generating function The midpoint method - The method of the right knots
1 2 3 4 Time 5 6
Figure 3: The Weibull-Gnedenko renewal density
In (Fig. 3) the curves for the analytical and all discrete solutions are shown for the Weibull-Gnedenko distribution with a = 1, p = 4 and the number of knots n = 100. This increase of collocation knots gives a practical coincidence of approximate solutions found using all the above methods. While for less amount of collocation knots (Fig. 2) the agreement between the results is not ideal. One can see, that all presented methods give oscillations of the renewal density which tends to the asymptotics for big values of t. This can be explained by the specifics of the renewal process, when the underlying lifetime distribution has rather big shape parameter p. The higher the value of p, the more often failures of the system take place. So, the operation and renewal process become unstable. The curve of the corresponding PDF is also shown in Fig. 3. As it can be seen from the formulas (4), (6), (7) the first member of the series representing the solution is equal to the value of the PDF for t = 0. It explains the proximity of the presented solutions to the PDF curve in the vicinity of t = 0. So, if one needs to know the system behavior in the initial period of time, it is sufficient to evaluate the PDF of the underlying lifetime distribution. The parameters of this distribution can be obtained from the statistical data processing.
3 Results & Discussion
The research techniques of restorable systems and their components renewal density are presented. They are based both on analytical and discrete methods and they consider the dependence of the renewal density on time. Using the Weibull-Gnedenko distribution some peculiarities of the renewal density behavior can be investigated.
The variation of the shape parameter p allows observation of the restorable system properties. As it was mentioned the higher the value of p the more often failures of the system take place, consequently, the less stable its behavior. For relatively small values of p the renewal density does not oscillate at all, that is the system operates in normal mode. This can be seen in Fig. 2. For larger values of p (Fig. 3) the oscillations are observed for some period of time, after which the renewal density goes to the asymptotics, corresponding to normal operation mode of the system. In [6] it was shown that for even higher values of p the oscillations period is longer. And this fact can be easily checked using the approaches of the present paper. Also the influence of the scale parameter a on the solution can be studied. Its variation gives the expansion or contraction of the curves over the time axis. So, a sort of criterion for restorable systems analysis can be developed if the failures statistics is approximated by the Weilbull-Gnedenko distribution. Unstable oscillating period is for sure undesirable for any application. The length and shape of this period are regulated by the parameters of the Weibull-Gnedenko distribution. The knowledge of these parameters values for a given restorable system can give a recommendation on its exploitation. It worth mentioning also, that the time scale in the renewal function dependences shown in Fig. 2 and Fig. 3 is conditional and should be adopted for each application separately. Practical significance of research results was demonstrated on several examples of processing of real statistical data on technologically active elements in gas supply systems failures [2, 19].
The advantages of the considered methods include their simplicity of algorithms and calculations. Nevertheless, one should keep in mind that the analytical approach (3) is valid only for the values of p higher than unity, which correspond to the most actual degradation period of restorable systems life cycles. Moreover, expansion (3) represents the asymptotic series, which does not converge uniformly. The summation of series (3) should be restricted by several members, giving the least relative error. This property was taken into account in curves plotting for Fig. 2 and Fig. 3. The presented discrete methods are more universal, thought their application is more time consuming.
References
[1] Elsayed E. A. Reliability Engineering. — John Wiley and Sons, Hoboken, 2012.
[2] Rusev V. N., Skorikov A. V. Analysis of elements of gas supply systems by method of moment generating functions. — Proceeding of Gubkin Russian State University of Oil and Gaz. — 2016. — no. 1(282). — P. 68-79. (In Russian).
[3] Smith W. L., Leadbetter M. R. On the Renewal Function for the Weibull Distribution. -1963. -Technometrics. - Vol. 5. — P. 393-396.
[4] Rinne H. The Weibull distribution. A handbook. — London, New York: CRC press, Taylor and Francis Group, 2009.
[5] Grigoriev L., Kucheryavy V., Rusev V., Sedyh I. Formation of estimates of reliability indicators for active elements in gas transport systems on the basis of refusals statistics. — Journal of Polish Safety and Reliability Association, Summer Safety and Reliability Seminars. — 2014. — Vol. 5, no. 1-2. — P. 41-47
[6] Constantine A. G., Robinson N. I. The Weibull renewal function for moderate to large arguments. — Computational Statistics and Data Analysis. — 1997. — Vol. 24. — P. 9-27.
[7] Cui L., Xie M. Some normal approximations for renewal function if large Weibull shape parameter. — Communications in statistics - -2003. — Vol. 32, no.1. —P. 1-16.
[8] Smeitink E., Dekker R. A Simple Approximation to the Renewal Function. — IEEE Transactions on Reliability — 1990. — Vol. 39, no.1. — P. 71-75.
[9] Maghsoodloo S., Helvaci D. Renewal and Renewal-Intensity Functions with Minimal Repair — Hindawi Publishing Corporation. Journal of Quality and Reliability Engineering — 2014.
— Vol. 2014, Article ID 857437. —10 Pages. http://dx.doi.org/10.1155/2014/857437.
[10] Krein M. G. On one extrapolation problem of A. N. Kolmogorov. —Doklady Akad. Nauk SSSR. — 1944. — no. 46(8). — P. 339-342. (In Russian).
[11] Sukharev M. G. Modeli nadezhnosti markovskogo tipa s prilozheniyami k neftegazovomu delu. — Moscow, Izdatel'skiy tsentr RGU nefti i gaza imeni I.M. Gubkina, 2012,
— 132 p. (In Russian).
[12] Deligonul Z. S., Bilgen S. Solution of the Volterra Equation of Renewal Theory with the Galerkin Technique using Cubic Spline. -Journal of Statistical Computation and Simulation — 1984. — Vol. 20, no.1 —P. 37-45.
[13] Xie M. On the solution of the renewal-type integral equations. —Communications in Statistics. —1989. — Vol. 18, no.1. — P. 281-283.
[14] Boehme T. K., Preuss W., van der Wall V. On a simple numerical method for computing Stielltjes integrals in reliability theory. — Probability in the Engineering and Informational Sciences — 1991. — Vol. 5, no.1. — P. 113-128.
[15] Tortorella M. Numerical solutions of renewal-type integral equations. —Informs Journal on Computing. — 2005. — Vol. 17, no.1. — P. 66-74.
[16] Jardine A. K. S., Tsang A. H. C. Maintenance, Replacement, and Reliability: Theory and Applications. Second Edition. — Boca Raton, CRC/Taylor & Francis. — 2013. — 330 P.
[17] Brezavscek A. A simple discrete approximation for the renewal function. —Business systems research. — 2013. — Vol. 4, no.1. — P. 65-75.
[18] Strelayev Y. M., Klimenko M. I. Primenenie metoda konechnyh elementov k resheniyu integralnyh uravneniy Volterra vtorogo roda. — Visnik ZNU. Fiziko - matematichni nauki. —2011.
— Vol. 18, no.2. —P. 131-135. (In Russian).
[19] Rusev V. N., Skorikov A. V. Analytical and discrete methods of research of failure intensity in gas transport. — Proceeding of Gubkin Russian State University of Oil and Gaz. — 2016. — no. 3(284). — P.104-117. (In Russian)
