Renewal and Renewal-Intensity Functions with Minimal Repair

The renewal and renewal-intensity functions with minimal repair are explored for the Normal, Gamma, Uniform, and Weibull underlying lifetime distributions. The Normal, Gamma, and Uniform renewal, and renewal-intensity functions are derived by the convolution method. Unlike these last three failure distributions, the Weibull except at shape β = 1 does not have a closedform function for the n-fold convolution. Since the Weibull is the most important failure distribution in reliability analyses, the approximate renewal and renewal-intensity functions of Weibull were obtained by the time-discretizing method using theMeanValue Theorem for Integrals. A Matlab program outputs all reliability and renewal measures.


Introduction
Renewal (RNL) functions give the expected number of failures of a system (or a component) during a time interval and this is used to determine the optimal preventive maintenance schedule of a system [1].Renewal functions (RNFs) have particular importance in analysis of warranty [2][3][4][5][6].They have wide variety of applications in decision making such as inventory theory [7], supply chain planning [8,9], continuous sampling plans [10,11], insurance application, and sequential analysis [2,8,12,13].
Since RNFs play an important role in many applications, it is important to obtain them analytically, if possible.Based on analytical approach, the RNF () is the inverse of its Laplace transform (), where Laplace transforms will be defined later.Blischke and Murthy [6] state that "the advantage of analytical method is that one can carry out parametric studies of the RNF, that is, the behavior of () as a function of the parameters of the distribution." However, for most distribution functions, obtaining the RNF analytically is complicated and even impossible [8].Therefore, development of computational techniques and approximations for RNFs has attracted researchers [14].
One of the well-known approximations given by Täcklind [15] is () ≃ /  1 +   2 /2 2 1 − 1, where   1 and   1 are the first and second raw moments, which is generally known as the asymptotic approximation and is also cited in numerous papers such as Smith [16].The asymptotic approximation has a closed-form expression and thus it is easy to apply to optimization problems that involve a RNL process [8].However since the asymptotic approximation is not accurate for small values of , Parsa and Jin [8] propose better approximation by keeping the positive features of asymptotic approximations such as simplicity, closed-form expression, and independence from higher moments of an underlying distribution.Jiang [17] proposes an approximation for the RNF with an increasing failure rate (IFR) which is also useful in areas such as optimization where a RNF needs to be evaluated.
There are series methods available in the literature to approximate RNFs such as Smith and Leadbetter [18] who developed a method to compute the RNF for the Weibull by using power series expansion of   where  is the shape parameter of the Weibull.On the other hand, instead of using power series expansion, Lomnicki [19] proposes another method by using the infinite series of appropriate Poissonian functions of   .There are also many other approximations available such as Xie [20], Smeitink and Dekker [21], Baxter et al. [22], Garg and Kalagnanam [23], and From [24].There are usually three criteria: model simplicity, applicability, 2 Journal of Quality and Reliability Engineering and accuracy to evaluate the value of RNF approximation [17].Increasing the complexity may lead to more accurate approximation but may make the process complicated and difficult to implement in practice [25].
Furthermore, in the literature, lower and upper bounds on RNFs have been discussed, which can be used to obtain upper and lower bounds on warranty costs such as Blischke and Murthy [6].Marshall [26] provides lower and upper linear bounds on the RNF of an ordinary RNL process.Ayhan et al. [27] provide tight lower and upper bounds for the RNF which are based on Riemann-Stieljes integration.There are also many other studies conducted about bounds on RNF such as Barlow [28], Leadbetter [29], Ozbaykal [30], Xie [31], Ran et al. [32], and Politis and Koutras [33].
Finally, simulation can be considered as an alternate approach to estimate the value of a RNF. Brown et al. [34] use the Monte Carlo simulation to estimate the RNF for a RNL process with known interarrival time distribution.
This paper proposes a convolution method to obtain the Gamma and Normal renewal and renewal-intensity functions; throughout, we are assuming that repair time is negligible (or minimal) relative to Time to Failure (TTF).As a result of this assumption, the replacement time of a failed component in a system is minimal.A further objective is to obtain the renewal and RNL-intensity functions of the uniform distribution by using  = 2 to  = 12 convolutions and applying the normal approximation for convolutions beyond 12.Because the Weibull distribution, except at shape  = 1, does not have a closed-form function for its -fold convolution, the last objective is to approximate its renewal and RNL-intensity functions by discretizing time using the Mean-Value Theorem for Integrals.We will also highlight the differences between the renewal-intensity function () and the hazard function ℎ().

Preliminaries
Suppose that component (or system) failures occur at times   ( = 1, 2, 3, 4, . ..) measured from zero and assuming that replacement (or restoration time) is negligible relative to operational time, then   represents the operating time (measured from zero) until the th failure.Because the probability density function (pdf) of  1 (or  1 ) may be different from the intervening times  2 = ( 2 −  1 ),  3 = ( 3 − 2 ),  4 = ( 4 − 3 ), . .., we consider only the simpler case of pdf of time to first failure  1 () being identical to those of intervening times  2 ,  3 ,  4 , . ... That is, all   's are assumed iid (independently and identically distributed).Therefore,   = ∑  =1   = sum of the times to the first failure from zero plus the intervening times of 2nd failure until the th failure.If  > 60, then the central limit theorem (CLT) states that the distribution of   approaches normality with mean , where  (the expected time between successive renewals) = (  ),  = 1, 2, 3, 4, . . .and with variance of   , denoted by (  ), equal to  2 .However, if the pdf of   , (), with  being the parent variable, is highly skewed and/or  is not sufficiently large, then the exact pdf of   = ∑  =1   is given by the -fold convolution of () with itself denoted by    () =  () ().Thus, by renewal we mean either the replacement of a failed component with a brand-new one, or the case when the failed component can be repaired in a negligibly short time to an almost brand-new condition.It has been proven by many authors both in stochastic processes and reliability engineering that the RNF for the duration [0, ] is given by where the random variable () represents the number (or count) of renewals that occur during the time interval [0, ] and  () () is the cumulative distribution function (cdf) of the -fold convolution of () with itself; that is,  () () = ∫  0  () (), where  () () is the -fold convolution density of ().It is also very well known that the RNI (renewal intensity) of a distribution function () by definition is given by () = ()/.
Authors in stochastic processes refer to () as the renewal density because () × Δ describes the (unconditional) probability element of a renewal during the interval (,  + Δ) [35]; further explanation is forthcoming in Section 5.2, while a few authors in reliability engineering, such as Ebeling [36] and Leemis [37], refer to () as intensity function.Because () is never a pdf over the support set of () for all failure distributions, throughout this paper we will refer to it as the renewal intensity function (RNIF).
The simplest and most common renewal process is the homogeneous Poisson process (HPP), where the intervening times are exponentially distributed at the constant interrenewal (or failure) rate .Because  is a constant and intervening times are iid, a Poisson process is also referred to as a homogeneous renewal process.It is also very well known that for a HPP, the pdf of interarrivals   's is given by () =  − , and that of the time to th failure (or renewal), measured from zero, is given by    () =  () () = (/Γ())() −1  − (the Gamma density with shape  and scale  = 1/).As a result, the use of (1a) for the interval [0, ] leads to the RNF: = , a fact that has been known for well more than a century.Further, the RNIF for a HPP is also a constant and is given by () = ()/ = ()/ = .Further, it is also widely known that for a HPP the [()] = .process variance  2 .Then, statistical theory indicates that time to the th failure (measured from zero) is the -fold convolution of (,  2 ) with itself, that is, time-to-the-thfailure   = TTF  ∼ (, 2 ).Hence, in the case of minimal repair, from (1a), the RNF is given by where Φ universally stands for the cumulative distribution function (cdf) of the standardized normal deviate (0, 1), and  () () gives the probability of at least  renewals by time .Cui and Xie [38] gave the same exact expression for the Normal RNF as in the above equation ( 2), which they used as an approximation for the Weibull renewal with shape  ≥ 3. It should be kept in mind that the normal failure law is approximately applicable in reliability analyses only if the coefficient of variation CV T ≤ 15% because the support for the normal density is (−∞, ∞), while TTF can never be negative (this assures that the size of left-tail below zero is less than 1 × 10 −10 ).If the CV T is not sufficiently small, then the truncated Normal can qualify as a failure distribution.[24] discussed the RNF for the truncated Normal.From a practical standpoint, the restriction on CV T can be relaxed to less than 20%.

The Renewal Function for a Gamma Baseline Failure Distribution
In order to obtain the RNF and RNIF of a Gamma nonhomogeneous Poisson process (NHPP), we first resort to Laplace transformation because it is the common procedure in stochastic processes.It has been proven (see [40, pages 277-280]) in theory of stochastic processes that the Laplace transforms (LTs) of () and () are, respectively, given by Further, it is also widely known that the LT of the Gamma density However, these last closed-form expressions for LTs are valid only if the shape  is an exact positive integer because the integration by parts does not terminate for noninteger values of .As a result, when the shape parameter  is not a positive integer, there exists no closed-form inverse-Laplace transform for the Gamma (); however, when the underlying failure distribution is Erlang (i.e., Gamma with positive integer shape), then there exists a closed-form inverse Laplace-transform for  = 2, 3, 4, . ... In fact, for the specific Erlang density with shape  ≡ 2 and scale  = 1/, it is very well known that () = ∫ ∞ 0  − () =  2 / 2 ( + 2) = −1/4 + /2 2 + 1/4( + 2).Upon inversion to the t-space, we obtain the well-known We used Matlab's ilaplace function as an aid in order to obtain the RNF and RNIF for the Erlang at shapes  ≡ 3 and 4 which are, respectively, given below The Gamma HZF at  ≡ 3 is given by ℎ , which exceeds the above () at  ≡ 3 for all  > 0.

Examining Results of the Convolution Method. Jin and
Gonigunta [25] propose the use of generalized exponential functions to approximate the underlying Weibull and Gamma distributions and solve RNFs using Laplace transforms.Table 1 shows a comparison of the RNF from (7) with their results.They refer to () as the actual RNF [obtained from Xie's [20] numerical integration] and   () is their approximated RNF.Table 1 shows that their () becomes more accurate as compared to () from (7) for larger values of .
Further, it is well known from statistical theory that the skewness of Gamma density is given by  3 = 2/√ and its kurtosis is  4 = 6/, both of these clearly showing that their limiting values, in terms of shape , is zero, which are those of the corresponding Laplace-Gaussian (/, / 2 ).We compared our results, using our Matlab program, for the Gamma at  = 70,  = 15, and  = 5000 which yielded (5000) = 4.186155 (Matlab gives 15 decimal accuracy), while the corresponding Normal yielded (5000) ≅ 4.185793 expected renewals.

The Renewal and RNI Functions When the
Underlying TTF Distribution Is Uniform whereas () represents the cdf of  = TTF.However, Hildebrand [41] proves that () × () = L{∫
In order to obtain the RNF for the Uniform density, we substitute into (9a), for the specific Uniform (0, ) baseline distribution, for which minimum life  = 0, in order to obtain () = / + ∫  0 ( − )/; letting  −  =  in this last equation yields The above equation (10) shows that the RNIF is given by () = ()/ = 1/ + ()/.Solving this last differential equation and applying the boundary condition ( = 0) = 0, we obtain () =  / − 1, where time must start at zero (or the last RNL); that is, this last expression is valid only for 0 ≤  ≤ ,  > 0. Note that Ross [39, Problem 3.7, page 154] gives the same identical () only for the standard (0, 1) underlying failure density, whose base is equal to 1. (It is worth mentioning the fact that Ross [39, page 7] uses the common notation in stochastic processes of () = 1 − () for the reliability function at lifetime ; the authors of this paper are indebted to Ross's text that has provided so much help and guidance in writing this paper).
Because the Uniform renewal function () =  (−)/ − 1 is valid only for the interval [, ], we will obtain the -fold convolution of the (, )-density which in turn will enable us to obtain () for  >  by making use of (1a) that uses the infinite-sum of convolution cdfs,  () ().As stated by Olds [43], the convolutions of Uniform density of equal bases, , have been known since Laplace.However there are other articles on this topic since 1950.The specific convolutions of the Uniform density with itself over the interval [−1/2, 1/2] were obtained by Maghsoodloo and Hool [44] only for  = 2-fold, 3, 4, 5, 6, and 8-folds.There are other articles on the Uniform convolutions such as Shiu [45], Killmann and Collani [46] and Kang et al. [47].We used the procedure in Maghsoodloo and Hool [44] but re-developed each of the  = 2 through  = 8 convolutions of (, ) by a geometrical mathematical statistics method (the (1983) article is available from the authors, although the 8th convolution has some typos).Further, we programmed this last geometric method in Matlab in order to obtain the exact 9 through 12-fold convolutions of the (, ) with itself.For example, our 7fold convolution-density of (, ) with itself is given below, where  =  −  > 0  (7)

The Uniform Approximation Results
. The method we propose uses the exact  = 2 through 12 convolutions  () () and then applies the Normal approximation for convolutions beyond 12.The question now arises how accurate is the normal approximation for  = 13, 14, 15, . ..? We used our 12-fold convolution of the standard Uniform (0, 1) to determine the accuracy.Clearly, the partial sum  12 = ∑ 12 =1   , each   ∼ (0, 1) and mutually independent, has a mean of 6 and variance 12( − ) 2 /12 = 1, where  = 0, and maximum life  = 1.The summary in Table 2 shows the normal approximation to  (12) () for intervals of length 0.50− .Table 2 clearly shows that the worst relative error occurs at /2 and that the normal approximation improves as  moves toward the right tail.The accuracy is within 2 decimals up to one- and 3 decimals beyond 1.49.Therefore, we conclude that the normal approximation to each of Uniform  () (),  = 13, 14, 15, . .., due to the CLT, should not have a relative error at  = 0.50 exceeding 0.002960.
It should also be noted that the normal approximation to the Uniform  () (),  = 13, 14, 15, . . .must be very accurate from the standpoint of the first 4 moments.Because the skewness of -fold convolution of (, ) with itself is identically zero, which is identical to the Laplace-Gaussian (( + )/2, ( − ) 2 /12), and hence a perfect match between the first 3 moments of    () =  () () with those of the corresponding Normal.It can be proven (the proof is available on request) that the skewness of the partial sum   = ∑  =1   ,   's being iid like , is given by Journal of Quality and Reliability Engineering 7 Further, the kurtosis of   is given by where   ( = 2, 3, 4) are the universal notation for the th central moments.Equation (12b) clearly shows that the kurtosis of the uniform  () (), for  = 13, 14, and 15, respectively, is  4 = −1.20/13= −0.09231,−1.20/14 = −0.08571,and − 1.20/15 = −0.08000,with the amount −1.20 being the kurtosis of a (, ) underlying failure density.Thus, an  = 120 is needed in order for the kurtosis of  () () to be within 0.01 of Laplace-Gaussian (( + )/2, ( − ) 2 /12).Fortunately, summary in Table 1 clearly indicates that the Normal approximation is superior at the tails, where kurtosis may play a more important role than the middle of  () () density [48, page 86].

Approximating the Renewal Function for Unknown Convolutions
6.1.The Three-Parameter Weibull Renewal Function.Unlike the Gamma, Normal, and Uniform underlying failure densities, the Weibull baseline distribution (except when the shape parameter  = 1) does not have a closed-form expression for its -fold cdf convolution  () (), and hence (1a) cannot directly be used to obtain the renewal function () for all  > 0. When minimum life = 0 and shape  = 2, the Weibull specifically is called the Rayleigh pdf; we do have a closedform function for the Rayleigh () but it cannot be inverted to yield a closed-form expression for its ().
It must be highlighted that there have been several articles on approximating the RNF such as Deligonul [49], and also on numerical solutions of () by Tortorella [50], From [24], and other notables.The online supplement by Tortorella provides extensive references on renewal theory and applications.Note that Murthy et al. [42] provide an extensive treatise on Weibull Models, referring to the Weibull with zero minimum life as the standard model.Murthy et al. [42] state, at the bottom of their page 37, the confusion and misconception that had resulted from the plethora of terminologies for intensity and hazard functions.We will use the time-discretizing method used by Elsayed [1], and others such as Xie [20], with the aid of Matlab to obtain another approximation for ().

Discretizing Time in order to Approximate the Renewal Equation.
Because () = () + ∫  0 ( − )() and the underling failure distributions are herein specified, the first term on the right-hand side (RHS) of (), (), can be easily computed.However, the convolution integral on the RHS, ∫  0 ( − )(), except for rare cases, cannot in general be computed and has to be approximated.The discretization method was first applied by Xie [20], where he called his procedure "THE RS-METHOD, " RS for Riemann-Stieltjes.However, Xie [20] used the renewal type (9b) in his RS-METHOD, while we are discretizing (9a).
The first step in the discretization is to divide the specified interval (0, ) into equal-length subintervals, and only for the sake of illustration, we consider the interval (0,  = 5 weeks) and divide it into 10 subintervals (0, 0.50 week), (0.50, 1),. .., (4.5, 5).Note that Xie's method does not require equal-length subintervals.Thus, the length of each subinterval in this example is Δ = 0.50 weeks.As a result, ∫ The above equation ( 13) is similar to that of (7.10) of Elsayed [1], where his subintervals are of length Δ = 1.It should be highlighted that using the interval midpoints (instead of endpoints) inorder to attain more accuracy leads to more computational complexity.We first used the information (0) ≡ 0 at  = 10 to calculate the last term of (13);  (−1)/2 ()}, where (0.50) has been approximated.Continuing in this manner, we solved (13) backward-recursively in order to approximate ().The smaller Δ always leads to a better approximation of ().

Renewal-Intensity Function Approximation for the Weibull
Distribution.Next, after approximating the Weibull RNF, how do we use its approximate () to obtain a fairly accurate value of Weibull RNIF ()?Because () ≡ ()/ ≡ Lim Δ → 0 (( + Δ) − ())/Δ, then, for sufficiently small Δ > 0 the approximate () ≈ (( + Δ) − ())/Δ, which uses the right-hand derivative, and () ≈ (() − ( − Δ))/Δ, using the left-hand derivative.Because the RNF is not linear but strictly increasing, our Matlab program computes both the left-and right-hand expressions and approximates () by averaging the two, where  and Δ are inputted by the user.It is recommended that the user inputs 0 < Δ < 0.01 such that the probability of 2 or more failures during an interval of length Δ is almost zero.Further, our program shows that () < ℎ() if and only if the shape  > 1.
6.4.Time-Discretization Accuracy.The accuracy of above discretization method was checked in two different ways.
The above time-discretization-method, using the Mean-Value Theorem for Integrals, can similarly be applied to any lifetime distribution and should give fairly accurate results for sufficiently small subintervals.However, Table 3 clearly shows that even at Δ = 0.01, the computational time exceeds one minute and the corresponding error may not be acceptable; further, it is not a closed-form approximation.

Conclusions
This paper provided the exact RNF and RNIF for the Normal, Gamma, and Uniform underlying failure densities.We have devised a Matlab program that outputs nearly all the renewal and reliability measures of a 3-parameter Weibull, Normal, Gamma, and Uniform.We have quantified the differences between () and ℎ() for  > 0, except in the only case of CFR for which () ≡ ℎ() ≡ .Further, when ℎ() is an IFR, then () < ℎ() so that () <  −(0,) ,  > 0, for the four baseline failure distributions that we have studied.However, this is not quite consistent with () =  −(0,) given by some authors in reliability engineering, such as Ebeling [36], for a NHPP.Further, some authors such as Leemis [37] use the same notation () for the Weibull RNIF and also use () for the Weibull HZF (see his Example 6.5, page 165).
Similar works for other underlying failure densities are in the immediate future.Work is also in progress that incorporates nonnegligible repair-time requiring some knowledge ofconvolution of TTF distribution with that of time-torestore (TTR).Such a stochastic process is referred to as an alternating renewal process [51, page 350].

Table 1 :
Comparison of   () with () for the Gamma renewal function.