Remaining useful lifetime estimation and noisy gamma deterioration process

Reliability Engineering and System Safety 149 (2016) 76–87

Contents lists available at ScienceDirect

Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress

Remaining useful lifetime estimation and noisy gamma deterioration process Khanh Le Son a, Mitra Fouladirad a,n, Anne Barros b a b

Institut Charles Delaunay, Université de Technologie de Troyes, UMR CNRS 6279 STMR, 12 rue Marie Curie, 10010 Troyes, France Norwegian university of science and technology, NTNU NO-7491 Trondheim, Norway

art ic l e i nf o

a b s t r a c t

Article history: Received 4 March 2015 Received in revised form 17 December 2015 Accepted 19 December 2015 Available online 4 January 2016

In many industrial issues where safety, reliability, and availability are considered of first importance, the lifetime prediction is a basic requirement. In this paper, by developing a prognostic probabilistic approach, a remaining lifetime distribution is associated to the system or component under consideration. More particularly, the system's deterioration is modelled by a non-homogeneous gamma process. The model considers a noisy observed degradation data and by using the Gibbs sampling technique, the hidden degradation states are approximated and afterwards the system's remaining useful lifetime distribution is estimated. Our proposed prognosis method is applied to the Prognostic and Health Management (PHM) 2008 conference challenge data and the interest of our probabilistic model is highlighted. To point out the interest of the prognostic, a maintenance decision rule based on the remaining lifetime estimation results is proposed. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Remaining useful life estimation Non-homogeneous gamma process Gibbs sampling Monitoring Condition-based maintenance Preventive maintenance

1. Introduction Prognosis and lifetime prediction have been largely developed in theory and industry. The quantity of interest is often the Remaining Useful Lifetime (RUL) and it is estimated based on the degradation measures. These measures often provide information such as the failure time data, failure sensors measurements to assess the reliability of systems. The variety and number of sources of uncertainty (e.g. forecast, complex system, unknown degradation process) encourage a probabilistic approach to model the degradation phenomenon. The probabilistic tools which permits the deterioration modelling are the stochastic processes. In many cases such as accumulation of damage in a mechanical systems or the deterioration due to defective products in case of a production line or due to corrosion/erosion in a structural systems, the problem is to deal with a system with monotone increasing deterioration. The appropriate deterioration process to model such deteriorations is the Gamma process because the paths of this latter are discontinuous and they can be thought of as the accumulation of an infinite number of small shocks which it is often how the degradation occurs, refer to [1,2]. Another interest of this process is that it allows feasible mathematical developments. A survey on the Gamma process and its application in degradation n

Corresponding author. Tel.: þ 33325718072. E-mail address: [email protected] (M. Fouladirad).

http://dx.doi.org/10.1016/j.ress.2015.12.016 0951-8320/& 2016 Elsevier Ltd. All rights reserved.

modelling is given in [2]. The Gamma process is largely used in deterioration modelling (eg. [3–7]). In practice, the available data on the health of a system are often noisy due to the presence of different sensors or measurement errors. Usually the sensor measurements on a monotone deteriorating system are contaminated by noises which can totally conceal the monotonicity of the deterioration process. To model such observations and not to neglect the monotonicity of the real deterioration phenomenon in this paper we propose to use a noisy gamma process where the gamma distributed deteriorating states are contaminated with Gaussian noise. This process gather both the noisy observations and the monotone deterioration phenomenon. To deal with this process first the observations are filtered to separate the noise and the deterioration states. Once we have access to the monotonically deteriorating states, we can make life time prediction based on the gamma process properties and the corresponding feasible mathematical calculations. In this framework, it is possible to give a probability distribution for the RUL and hence to obtain a confidence interval for its estimation and measure the precision of our prognosis method. As a case study for the deterioration model, we consider the degradation indicator obtained from the 2008 Prognostic Health Management Challenge (PHM) data [8] and explored in [9]. The proposed indicator in [9] has a monotonous trend with a large noise. Therefore it seems sensible to model the deterioration indicator evolution by a noisy stochastic process. In this aim a non-homogeneous Gamma process with a Gaussian noise is used to model the deterioration. Hence, a stochastic filtering is

K. Le Son et al. / Reliability Engineering and System Safety 149 (2016) 76–87

theoretically proposed to find the hidden degradation states of the gamma process. Gamma process with additive Gaussian noise has been already considered in the literature, refer to [10,11]. In both papers, the parameter estimation through the maximum likelihood calculation in the presence of unobserved deterioration is carried out by numerical approximations. Authors in [10] consider a Monte Carlo method whereas in Lu et al. [11] use the Quasi Monte Carlo method to reduce the calculation time. The Quasi Monte Carlo method proposed in [11] approximates the likelihood function by using the Genz transform. In comparison with [10] the proposed Quasi Monte Carlo method leads to a smaller error and reduces the computation time. In a general framework to estimate the non-observable system state, the stochastic filtering approaches are frequently used, see [12,13]. The stochastic filtering approaches give an estimation of the system state based on all collected measures history, therefore these methods can give reliable performances for the degradation states estimation. In the framework of stochastic filtering methods, the well-known Kalman filter is essentially used on linear models with Gaussian noises. In [14], the authors introduce a RUL estimation model based on Kalman filter for the aircraft engine. In industry, we usually deal with a system with non-linear behavior and non-necessarily Gaussian noisy observations. Following [15], the non-linear systems are considered in many areas of science and engineering. One of most popular non-linear filtering techniques is the particle filtering, a sequential Monte Carlo technique. We can find the application of the particle filtering in different domains such as signal processing [16], biology, biostatistics [17], etc. In this paper, we propose another Monte Carlo Markov Chain (MCMC) technique, Gibbs sampling, for stochastic filtering that can be applied to non-linear models and non-Gaussian noises. Recently, this technique is used in [18] to model the survival data (semi-competing risks model) via a gamma frailty. Gibbs sampling can be applied to the models with unknown parameters and can estimate stably the target function (unobserved data) based on the observed data. This sampling method which requires heavy calculations give a precise estimation of the real state and also its associated uncertainty. The use of the Gibbs sampling in the framework of the non-homogeneous gamma process is not a trivial task and this method is not applied before in prognosis. In this paper the use of Gibbs sampling technique is to extract the Gamma degradation states from noisy degradation indicator. The likelihood function is complex due to the hidden degradation states, therefore the estimation of the model parameters by maximizing the likelihood directly is infeasible. As proposed in [19], the stochastic expectation–maximization (SEM) algorithm is matching to solve this problem. The algorithm can iteratively take the approximated values of hidden data to the complete likelihood function and then estimate parametric inference of model. After the states estimation step, based on the filtered data a probability distribution is associated to the remaining useful life distribution. In previous works on the same data set (PHM data) a Wiener process was proposed to model the deterioration [9]. In the case of the Wiener process it is supposed that the large variability of data is not due to the presence of noise but to the properties of the system under consideration. Since rarely deteriorating systems have a non-monotone tendency this modeling does not remain valid for a big range of observation data. Another weakness of the use of Wiener process in the deterioration modeling is the unfeasible probability and mathematical calculations which lead us to solve our prognosis problem by Monte Carlo simulations and the prognosis is based on the mean value of the RUL (see [9]). In this paper we propose a more efficient and precise prognosis method which is accordingly more time consuming. The preliminary results of this method are presented in [20,21]. In this

77

paper once the RUL is estimated a predictive maintenance is proposed to prevent system's failure and take advantage of the prognosis procedure. The remainder of the paper is organized as follows. Section 2 presents the deterioration model using a non-homogeneous gamma process with Gaussian noise and how to filter the hidden degradation states. In section 3 a method for the remaining useful life estimation is proposed and the impact of the observations vectors size on the RUL estimation is studied. Section 4 is devoted to a case study obtained from the 2008 Prognostic Health Management Challenge data. In Section 5, a maintenance decision rule based on the remaining useful life estimation is proposed. Finally, the conclusions of the results as well as the furthers works are given in Section 6.

2. Deterioration modelling and state estimation 2.1. Deterioration model Let us consider a gradually deteriorating system and we suppose that a scalar random variable can summarize the system state at time t. The system is supposed to have a monotone nondecreasing degradation which is the behavior observed in many physical deterioration processes. Denote by X(t) the system state at time t which is monotone non-decreasing. A survey on Gamma process given in [2] pointed out this process to model a monotone increasing deterioration. The paths of the Gamma process are discontinuous and represent the accumulation of an infinite number of small shocks which is often how the degradation occurs, refer to [1,2]. Let X(t) be a gamma distributed random variable with a scale parameter β 4 0 and shape function vðtÞ ¼ αt b which is a nondecreasing, right continuous, real-valued function for t 4 0 with α 4 0, XðtÞ  Γ ðvðtÞ; βÞ. The probability density function is given by: f ðxj v; β Þ ¼

βvðtÞ vðtÞ  1 x expð  βxÞI½0; þ 1½ ðxÞ; Γ ðvðtÞÞ

R1 where Γ ðxÞ ¼ z ¼ 0 zx  1 e  z dz is the Euler gamma function and I½a;b ðxÞ ¼ 1 if x A ½a; b and I½a;b ðxÞ ¼ 0 otherwise. The gamma process with shape function vðtÞ 4 0 and scale parameter β has the following properties:

 Xð0Þ ¼ 0,  XðsÞ  XðtÞ  Γ ðvðsÞ  vðtÞ; βÞ for all s 4t Z0,  X(t) has independent increments. It is supposed that the deterioration level is not directly observable. Let Y j ¼ Yðt j Þ; j ¼ 1; …; n, be the degradation indicators observed at inspection times 0 o t 1 o … ot n , and X j ¼ Xðt j Þ; j ¼ 1; …; n is the non-observable state at time tj modelled by a nonhomogeneous gamma process, thus the relation between Xj and Yj can be expressed as follows: Y j ¼ f ðX j ; ϵj Þ ¼ X j þ ϵj where ϵj are independent Gaussian random variables with standard deviation σj and mean equal to zero and it can be also expressed as a function of Xj and Yj as follows: ϵj ¼ gðX j ; Y j Þ ¼ Y j X j . For an efficient prognosis and maintenance planning it is necessary to estimate the non-observable state of the system. In this aim, the hidden degradation states vector X ¼ ðX 1 ; …; X n Þ should be evaluated from the observations Y ¼ ðY 1 ; …; Y n Þ. Due to uncertainties related to the vector X and to the observations Y, to propose a relevant estimation of X and to be able to analyse the estimation performances we should first calculate the conditional

78

K. Le Son et al. / Reliability Engineering and System Safety 149 (2016) 76–87

density of the non-observable states X 1 ; …; X n with respect to the observations Y 1 ; …; Y n defined as follows:

μX=Y 1 ¼ y1 ;…;Y n ¼ yn ðx1 ; …; xn Þ ¼ R

μX;Y ðx1 ; …; xn ; y1 ; …; yn Þ R ⋯ μX;Y ðx1 ; …; xn ; y1 ; …; yn Þdx1 …dxn ð1Þ

with ðx1 ; …; xn Þ and ðy1 ; …; yn Þ a realisation of X and Y, respectively. The density μX;Y ðx1 ; …; xn ; y1 ; …; yn Þ can be deduced as follows:

μX;Y ðx1 ; …; xn ; y1 ; …; yn Þ ¼ Ce  βxn n

 ∏ ðxj xj  1 Þαðt j

b

  

  

 t bj 1 Þ  1  ðg 2 ðxj ;yj ÞÞ=ð2σ 2 Þ ∂gðxj ; yj Þ i

e

j¼1

∂y

ð2Þ

where C is a constant and ∂gð.,yÞ=∂y is the partial derivative of g with respect to the second variable y. Hence, the conditional density of the hidden degradation states from the observations can be rewritten by:

μX=Y 1 ¼ y1 ;…;Y n ¼ yn ðx1 ; …; xn Þ ¼ K 1 e  βxn n

 ∏ ðxj xj  1 Þαðt j

b

  

  

 t bj 1 Þ  1  ðg 2 ðxj ;yj ÞÞ=ð2σ 2 Þ ∂gðxj ; yj Þ i

e

j¼1

∂y

ð3Þ

with respect to the observations μX=Y ðx1 ; …; xn Þ can be approximated and the unfeasible constant K1 calculation can be avoided. 2.2. Filtering of the hidden degradation states Gibbs sampler presented in [22,23], is a Monte Carlo Markov Chain (MCMC) algorithm used for generating random variables from a (marginal) distribution without having to calculate the density by using elementary properties of Markov chains. It is often used in statistical problems to approximate the likelihood function when this latter is difficult to be calculated in the closed form. In the Bayesian estimation problems, the Gibbs sampler is often used to generate non-conjugated posterior distributions. In this study, an ergodic Markov chain is generated whose invariant distribution is μX=Y . At each step q, q A N, the Gibbs algorithm gives the successive values of a Markov chain Zq ¼ ðZ q1 ; …; Z qn Þ, in which the hidden gamma degradation states X 1 ; …; X n can be approximated. In step q þ 1, we generate successively Z q1 þ 1 ; …; Z qj þ 1 ; …; Z qn þ 1 according to the following probability distribution:

 Generating the value zq1 þ 1 of Z q1 þ 1 according to the law of z1:

where K1 is the coefficient defined as follows: 1 ¼ K1

Z

Z

þ1

⋯

0

þ1

e



∂y

j¼1

0

ð4Þ Due to the presence of large number of integrals in Eq. (4) the calculation of the constant K1 is tricky and some times even unfeasible. To avoid this difficult calculation, the Gibbs sampler algorithm is proposed to estimate the conditional density μX=Y [22]. This algorithm is based on the marginal conditional densities for each component Xj of X given the other components. In our case we get:

 For j ¼ 1; α

μX=Y 1 ¼ y1 ;…;Y n ¼ yn ðx1 =x2 ; …; xn Þ ¼ K 2;1 x1

ðt b1 Þ  1

  2 2 ∂gðx1 ; y1 Þ 1ð0 o x o x Þ e  g ðX 1 ;y1 Þ=2σ 1  1 2  ∂y

where K 2;1 ¼ R

¼ y1 ;…;Y n ¼ yn ðx1 =x2 ;…;xn Þ

ðx2  x1 Þαðt 2  t 1 Þ  1 b

b

ð5Þ

1

μX=Y 1

dx1

.

 For 2 r jr n 1, μX=Y 1 ¼ y1 ;…;Y n ¼ yn ðxj =x1 ; …; xj  1 ; xj þ 1 ; …; xn Þ ¼ K 2;j ðxj  xj  1 Þαðtj  t j  1 Þ  1 b

ðxj þ 1  xj Þαðt j þ 1  t j Þ  1 eð  g b

b

where K 2;j ¼ R

μX=Y 1

 For j¼n

2

b

  ∂gðxj ; yj Þ   ∂y 1ðxj  1 o xj o xj þ 1 Þ

ðxj ;yj Þ=2σ 2 Þ 

1 ¼ y1 ;…;Y n ¼ yn ðxj =x1 ;…;xj  1 ;xj þ 1 ;…;xn Þdxj

  

  1ðxn  1 o xn Þ

ðxn ;yn Þ=2σ 2 ∂gðxn ; yn Þ

where K 2;n ¼ R

∂y

μX=Y 1

1 ¼ y1 ;…;Y n ¼ yn ðxn =x1 ;…;xn  1 Þdxn

μX=Y ðzj j zq1 þ 1 ; …; zqj þ11 ; zqjþ 1 ; zqn Þ; 2 r j rn  1

 Generating the value zqn þ 1 of Z qn þ 1 according to the law of zn: μX=Y ðzn j zq1 þ 1 ; …; zqn þ 11 Þ

Eventually, a Markov chain with the invariant distribution approximating μX=Y is obtained. The number of Gibbs iterations Q0 should be cautiously chosen in order to obtain a good approximation of ðX 1 ; …; X n Þ by the vectors Zq ¼ ðZ q1 ; …; Z qn Þ, q ¼ 1; …; Q 0 . In the case of our model, the approximations of hidden degradation states are obtained by using the conditional densities defined in (5)–(7) according to the two following probability distribution functions:

 For 0 rc r z rd 

b

ð7Þ .

The Gibbs sampler algorithm is presented more in detail in Section 2.2. By using the Gibbs sampler the problem of calculation of a joint density function is reduced to the calculation of several uni-dimensional distributions. As it can be noticed in Eqs. (5)–(7) the constants K 2;j , j¼1,…,n depend only on a simple integral calculation which is much easier to implement than the multiple integral calculation required for the constant K1 in Eq. (4). By calculating the marginal conditional distributions (5)–(7) in the framework of a Gibbs algorithm, the conditional density function of the non-observable states

 g 2 ðz; YÞ ð8Þ j g 0 ðz; YÞj 2 2σ Rd where C1 is the normalisation constant, c μ1 ðzÞdz ¼ 1=C 1 , γ and δ correspond to αðt bj  t bj 1 Þ  1 and αðt bjþ 1 t bj Þ  1 in the Eq. (6). For 0 rc r z   g 2 ðz; YÞ j g 0 ðz; YÞj μ2 ðzÞ ¼ C 2 ðz  cÞγ expð  βzÞexp  ð9Þ 2 2σ R þ1 where C2 is the normalisation constant, c μ2 ðzÞ dz ¼ 1=C 2 . The distribution μ1 defined in Eq. (8) is derived from Eqs. (5) and (6). When the aim is to draw Z1, μ1 is equal to Eq. (5) and therefore C 1 ¼ K 2;1 , c ¼0 and d ¼ z2 . When the aim is to draw Zj, j ¼ 2; …; n  1, μ1, is equal to Eq. (6) and therefore C 1 ¼ K 2;j , c ¼ zj  1 and d ¼ zj þ 1 . The distribution μ2, defined in Eq. (9) corresponds to the random variable Zn and is equal to the distribution proposed in Eq. (7). In this case C 2 ¼ K 2;n and c ¼ zn  1 . The distribution functions μ1 or μ2 are not known distribution functions and generating random variables according to these distribution functions is not a easy task. Therefore a rejection method as in [24] is used an detailed as follows.

μ1 ðzÞ ¼ C 1 ðz  cÞγ ðd  zÞδ exp 

. b

2

 Generating the value zqj þ 1 of Z qj þ 1 according to the law of zj:

ð6Þ

μX=Y 1 ¼ y1 ;…;Y n ¼ yn ðxn =x1 ; …; xn  1 Þ ¼ K 2;n e  βxn ðxn  xn  1 Þαðtn  tn  1 Þ  1 eg

μX=Y ðz1 j zq2 ; …; zqn Þ



b b 2 2 ∂gðxj ; yj Þ dx1 …dxn ∏ ðxj xj  1 Þαðtj  t j  1 Þ  1 e  ðg ðxj ;yj ÞÞ=ð2σ Þ  

n

 βxn



K. Le Son et al. / Reliability Engineering and System Safety 149 (2016) 76–87

 Generation of a random variable with the probability distribution μ1, with 0 r c r z r d: we realize a chain of drawing of lots

ðU 1 ðiÞ; U 2 ðiÞÞ independent with the uniform distribution on ½c; d  ½0; ðγ ðd  cÞ=δ þ γ Þγ ðδðd  cÞ=δ þ γ Þδ . If we suppose that I is the first index such as   g 2 ðU 1 ðIÞ; YÞ j g 0 ðU 1 ðIÞ; YÞj ; U 2 ðIÞ r ðU 1 ðIÞ  cÞγ ðd  U 1 ðIÞÞδ exp  2 2σ

so U 1 ðIÞ follows the probability distribution μ1.  Generation of a random variable with the probability distribution μ2, with z Z c note that   g 2 ðz; YÞ j g 0 ðz; YÞj r K 3 expð  β zÞ ðz  cÞγ expð  β zÞexp  2 2σ

ð10Þ

For where K 3 ¼ supz Z c ðz  cÞγ expð  g 2 ðz; YÞ=2σ 2 Þj g 0 ðz; YÞj . gðx; yÞ ¼ y  x, we have j g 0 ðz; YÞj ¼ 1 and K3 can be calculated as follows: maximizing the function ðz  cÞγ expð  g 2 ðz; YÞ=2σ 2 Þj g 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðz; YÞj , we get the solution zn ¼ y þ c þ ðy  cÞ2 þ 4γσ 2 =2 and therefore, K 3 ¼ ðzn  cÞγ expð g 2 ðzn ; YÞ=2σ 2 Þ. We then realize a chain of independent drawing of lots ðV 1 ðiÞ; V 2 ðiÞÞ with V 1 ðiÞ

following an exponential distribution of parameter β and V 2 ðiÞ having the conditional uniform distribution on ½0; K 3 expð  β V 1 ðiÞÞ. If wesuppose that I is the first index such as   g 2 ðV 1 ðIÞþ c; YÞ j g 0 ðV 1 ðIÞ þ c; YÞj : V 2 ðIÞ r ðV 1 ðIÞÞγ expð  βV 1 ðIÞÞexp  2σ 2 therefore V 1 ðIÞ þc follows the probability distribution

μ2.

2.3. An example of state estimation In Fig. 1, an example of Gibbs sampling applied to a Gamma process is given. Observations trajectory are generated by a Gamma process with parameters ðα ¼ 1; β ¼ 1; b ¼ 1:5Þ and an additive Gaussian noise parameter σ ¼ 5. In this case, we consider a Gamma process with simple parameters in order to evaluate the impact of Gibbs sampling on RUL estimation. Two different iteration numbers of Gibbs sampler are depicted in Fig. 1. A small number of iterations does not permit the approximation of the Gamma process. As the iteration number grows the Gibbs approximation get closer to the real state. In [25,26], methods to chose the iteration number of the Gibbs sampler are proposed. To study the impact of the noise on the state estimation, a gamma process with the parameters ðα ¼ 1; β ¼ 0:5; b ¼ 1:5Þ and a

79

Gaussian noise parameter with three different variances is considered σ A f2; 5; 10g and the state estimations are depicted in Fig. 2. Let the Root Mean Square Error (RMSE) of the degradation states approximations be defined as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pn ^ 2 j ¼ 1 ðxj  x j Þ RMSE ¼ : ð11Þ n where xj and x^ j are respectively the real and estimated degradation states at time tj. For three different cases of Gaussian noise, the RMSE of the degradation states approximations given by the Gibbs algorithm are presented in Table 1. The RMSE increases as the noise get larger. However, the RMSEs in all of the three cases seem to be not very large and could be acceptable. 2.4. Parameters estimation As the deterioration states are hidden to estimate the unknown model parameters, the maximum likelihood method cannot be directly used. Based on an approximation of the deterioration states obtained by the Gibbs sampling, the unknown parameters are estimated through the Stochastic Expectation–Maximization (SEM) algorithm [27,28]. The Nelder–Mead algorithm is proposed (see [29]) to maximize the log-likelihood function. Recall that the likelihood function of gamma distributed deteðiÞ rioration states XðiÞ ¼ ðX ðiÞ 1 ; …; X ni Þ; i ¼ 1; …; N (N: number of comðiÞ ponents) given the observations Y ðiÞ ¼ ðY ðiÞ 1 ; …; Y ni Þ; i ¼ 1; …; N is defined as follows:

Lðα; b; β; σ Þ ¼

ni  N X X i¼1j¼1





ðiÞ αðt bj  t bj 1 Þ  1 ln X ðiÞ j  Xj  1



ðiÞ g 2 ðX ðiÞ j ; Yj Þ

  ðiÞ þ ln j g 0 ðX ðiÞ j ; Y j Þj    pffiffiffiffiffiffi  lnðσ 2π Þ þ αðt bj  t bj 1 Þln β  ln Γ αt bj  αt bj 1

ðiÞ  βðX ðiÞ j  Xj  1Þ 

2σ 2

ð12Þ The random variables

X ðiÞ j ;j¼

1; …; ni in the log-likelihood are

approximated by the Gibbs sampling and replaced by Z ðiÞ j ; j ¼ 1; …; ni in the likelihood function to estimate iteratively the parametric inference of model by the SEM algorithm. The algorithm is described as the following steps: 1. Choose the initial parameters ðα0 ; β 0 ; b0 ; σ 0 Þ and then simulate Z1 ¼ ðZ 1i;1 ; …; Z 1i;ri Þ by Gibbs sampling. 2. At the q-th iteration, by the samplers vector Zq ¼ ðZ qi;1 ; …; Z qi;ri Þ, we estimate: PN Pri g 2 ðZ qi;j ; Di;j Þ σ 2q ¼ i ¼ 1 Pj ¼N 1 ð13Þ i ¼ 1 ri and ðαq ; βq ; bq Þ are obtained by maximizing the likelihood function L using Nelder–Mead algorithm, 3. At ðq þ 1Þ-th iteration, use ðαq ; βq ; bq ; σ q Þ to simulate Zq þ 1 ¼ ðZ qi;1þ 1 ; …; Z qi;rþi 1 Þ. 4. Parameters of Gamma process are updated to a sufficient number of iterations Q, a parameters set ðαq ; βq ; bq Þ; q ¼ 1; …; Q is obtained, 5. The estimated parameters are finally calculated as follows:

αb ¼ Fig. 1. Gamma filtering of a simulated time series by Gibbs algorithm.

Q Q Q 1 X 1 X 1 X αq ; βb ¼ βq ; bb ¼ bq : Q q¼1 Q q¼1 Q q¼1

80

K. Le Son et al. / Reliability Engineering and System Safety 149 (2016) 76–87

Fig. 2. Illustration of a gamma process with noise σ ¼ 2 (left), σ ¼ 5 (middle), σ ¼ 10 (right) and its Gibbs filtration.

Z

Table 1 Roots mean squared error of the Gibbs approximations

L

¼ 0

σ

2

5

10

RMSE

1.102

1.811

2.7121

F αððt n þ hÞb  tb Þ;β ðL  xÞ:μX n =Y 1 ;…;Y n dx n

ð14Þ

where F αððt n þ hÞb  tb Þ;β ðL xÞ ¼ 1  F αððtn þ hÞb  t b Þ;β ðL  xÞ; n

n

3. Remaining useful lifetime estimation 3.1. RUL formulation In this paper, we define by the RUL the during time between the observing time and the failure time with respect to the noisy observation. The RUL can be considered as a random variable with a conditional probability distribution. The system is considered to be failed when the stochastic deterioration process cross a failure threshold L. The threshold L can be random or non-random. We shall denote by RULðt n Þ the random variable corresponding to the remaining useful life of the observations trajectory at time tn, the probability distribution function of RULðt n Þ is defined as follows: F RULðtn Þ ðhÞ ¼ PðRULðt n Þ o hÞ ¼ PðX t n þ h 4Lj X n o L; Y 1 ; …; Y n Þ ¼ PðX t n þ h  X n 4 L  X n j X n o L; Y 1 ; …; Y n Þ Z L ¼ PðX t n þ h  X n 4 L  xj X n ¼ xÞ:μX n =Y 1 ;…;Y n dx 0

and F αððt n þ hÞb  tb Þ;β ðL  xÞ is the distribution function of Gamma n

process with shape function αððt n þ hÞb  t bn Þ and scale parameter β, and μX n =Y 1 ;…;Y n is the conditional density of Xn. The conditional distribution function of Xn, μX n =Y 1 ;…;Y n , above can be estimated by the Gibbs algorithm in the previous section. Eventually, the conditional distribution function presented in Eq. (14) is estimated by the Gibbs algorithm as follows: Q0 þ Q 1 X Fb RULðt n Þ ðhÞ ¼ F ðL  Z n ðqÞÞ b b Q q ¼ Q þ 1 αððtn þ hÞ  tn Þ;β

ð15Þ

0

where Zn(q) is the hidden degradation state approximation at time tn at the qth Gibbs iteration, Q0 is the necessary iteration number of Gibbs simulation to get a convergence state and Q is the required Gibbs simulation iterations to get a robust estimation. These iterations numbers are chosen by the variance ratio method proposed in [26].

K. Le Son et al. / Reliability Engineering and System Safety 149 (2016) 76–87

3.2. Validation of the RUL calculation method The mean associated to the distribution function defined in Eq. (14) is not reachable and it is approximated by the mean assod mean ðt n Þ. ciated to the distribution function Eq. (15) denoted by RUL In order to evaluate our proposed RUL estimation method we propose to compare at each time t n o t L the following quantities:

 RULreal ðt n Þ ¼ t L  t n .  RULmean ðt n Þ ¼ EðRULX ðt n ÞÞ the mean of the random variable RULX ðt n Þ obtained from the non-noisy deterioration with the following probability distribution function: F RULX ðtn Þ ðhÞ ¼ PðRULX ðt n Þ o hÞ ¼ PðX tn þ h 4 Lj X tn ¼ xn Þ ¼ PðX t n þ h  X tn 4 L  xn Þ ¼

Γ ðvðt n þhÞ  vðt n Þ; ðL xn ÞβÞ Γ ðvðt n þ hÞ  vðt n ÞÞ ð16Þ

where vðtÞ ¼ αt and Γ ðv; uÞ is the incomplete gamma function for u Z0 and v 40. d median ðt n Þ is the estimated median value of RULðt n Þ given by RUL b



d median Þ ¼ 0:5 Fb RULðtn Þ ðRUL where the distribution function Fb is estimated by Eq. (15).

d mean ðt n Þ is the estimated mean value of RULðt n Þ using the  RUL distribution function Fb defined by Eq. (15).

The quantity RULreal is calculated by using the training data set from 0 to tL. The quantity RULmean is calculated by using the training data set from 0 to tn and considering the future state of the system as unknown and random but not noisy. The quantities d median are calculated by considering only the d mean and RUL RUL estimated deterioration states and the noisy data set. By comparing these quantities, the performances of the RUL estimation method and the filtering method are investigated (Fig. 3).

81

As it is first presented in [20], in Fig. 4 the estimated RUL distribution F^ defined by Eq. (15) and the probability distribution function F defined in Eq. (16) are depicted and compared. For a better lecture of Fig. 4 the 95% and 5% quantiles are pointed out. It can be noticed that the two intervals have a large intersection, therefore it is reasonable to estimate the RUL distribution of the system or a component by considering F^ defined by Eq. (15). This result can be confirmed by the following implementation where d mean are compared. the RULmean and RUL With 100 sampled observation paths corresponding to 100 identical components, we study the 95% confidence interval of the RUL estimation. By checking the 100 considered components, it can be noticed in Fig. 5 that for all 100 components the RULmean is d mean obtained through the inside the 95% confidence interval of RUL ^ d mean can be estimated distribution F . Hence, in this paper RUL considered a robust estimator of RULmean. To study the consistency of our results, the estimated value d mean and the real value of RUL in the testing data set are RUL d mean and RULreal of 20 compocompared. The comparison of RUL nents is depicted in Fig. 6. The RULreal of only 2 components are d mean and the RULmean outside the 95% confidence interval of the RUL d mean and RUL d median are all in the 95% confidence interval. The RUL depicted in Fig. 6 are very close to the RULmean. We can deduce that d mean or it is sensible to estimate the real value of RUL by the RUL d . RUL median 3.2.2. Validation with a case study The data set, provided by the 2008 Prognostic Health Management (PHM) Data Challenge Competition, consists of multivariate time series that are collected from 218 identical and independent units of an unspecified system (the sensor

3.2.1. Validation with simulated data In order to evaluate the proposed probabilistic method, we consider 100 trajectories simulated by a non-homogenous gamma process with parameters ðα ¼ 1; β ¼ 1; b ¼ 1:5Þ and the failure threshold L ¼50. Each trajectory is simulated in time and it is finished when the degradation state attains the failure threshold L at time tL. We use the simulated deterioration from the beginning to failure (from t¼0 to tL) as the testing data test. By adding the Gaussian noise with parameter σ ¼ 3 on the gamma process (testing data), we obtain noisy observations.

Fig. 4. RUL distribution function for the non-observable and estimated degradation states.

Fig. 3. Gamma filtering and estimated degradation process. The failure threshold is given L ¼ 100 and RUL distribution is estimated at the inspection time.

Fig. 5. Estimated RULs interval of components and mean value of RULs obtained by the real degradation state.

82

K. Le Son et al. / Reliability Engineering and System Safety 149 (2016) 76–87

measurements are generated via a thermo-dynamical simulation model for the engine as a function of variations of flow and efficiency of the modules of interest [8]). Each time series is from a different instance of the same type of unit, the data might be from a fleet of ships of the same type. There are three operational settings that have a substantial effect on units performance. The data for each cycle of each unit include the unit ID, time cycle index, 3 values for the operational settings and 21 values for 21 sensor measurements connected to the units state. The sensor data are contaminated with noise. Two data sets are available: a training data set and a testing one. The prognostic method is proposed and constructed based on the training data set. The testing data set is used to estimate the RUL of each unit and then to evaluate the accuracy of the proposed prognostic method. For the training data set, operational settings and sensor measurements are complete from the starting cycle until the failure, and for the testing data set, operational settings and sensor measurements begin at the starting cycle until a given cycle upon which the RUL must be estimated. Failure time is also given to calculate the performance of the prediction methods. Deterioration modelling: The construction of degradation indicator obtained from the sensors measurements for the 2008 PHM Challenge data is presented in [30]. Let us consider the degradation indicator Di;½1:ni  ; i ¼ 1; …; 218 introduced initially in [30] and depicted in Fig. 7. Let be n ¼ 218 the total number considered components in the training data set and ni the number of

measurements corresponding to component i in the training data ðiÞ set. Let be Y~ ¼ Di;½1:ni  the given ni measurements of component i in the training data set. It can be noticed in Fig. 7 that the indicators Di;½1:ni  ; i ¼ 1; …; 218 are decreasing and in order to be able to use the deterioration model proposed in this paper the observations Y ðiÞ corresponding to the component i are defined as follows: ~ ðiÞ ~ ðiÞ Y ðiÞ j ¼ Y 1 Y j ;

j ¼ 1; …; ni :

ð17Þ

Each trajectory Y ðiÞ ; i ¼ 1; …; 218 corresponds to the life cycle of a component in the training data set. Considering each Y ðiÞ as the ðiÞ noisy measurement of the real state vector XðiÞ ¼ ðX ðiÞ 1 ; …; X ni Þ of component i we can use a non-homogeneous gamma process with ðiÞ noise to model the measurements Y ðiÞ 1 ; …; Y ni . Using the convergence criteria proposed by [26] and the SEM method, the number of iterations to attain a stable state is fixed to Q 0 ¼ 15; 000 and the model parameters are estimated and presented in Fig. 8. The estimated parameters are presented in Table 2. These values are used to calculate the RUL of the system. RUL estimation – Training data set: In the training data set let us estimate the model parameters and the components state at failure times. Since we observe a large variation on the failure state we can consider that the failure threshold L is a random variable. Using a goodness of fit to a Weibull distribution, the following probability density function is associated to L: f ðl; μ; λ; θÞ ¼ μλ

μ

  μ  lθ ; ðl  θÞμ  1 exp 

λ

ð18Þ

with λ ¼ 1:1498, μ ¼ 3:3430 et θ ¼ 1:2693. RUL estimation – Testing data set: Based on the proposed deterioration model in the previous section the aim is to estimate the RUL of the components. Considering a random failure threshold the cumulative distribution function (14) of the remaining useful life time at time tn is estimated as follows: Q0 þ Q Z 1 X F αððt n þ hÞb  tb Þ;β ðl Z qn Þ:f L ðlÞ dl Fb RULðt n Þ ðhÞ ¼ n Q q ¼ Q þ1

ð19Þ

0

Fig. 6. Comparing the estimated RULs interval of components with the values of RULs.

where fL(l) is the probability density function of L et Znq is the output of the Gibbs algorithm at the q-th iteration. To measure the efficiency of the proposed RUL estimation method the criteria proposed in PHM 2008 conference is

Fig. 7. The degradation indicator corresponding to component 1 (left) and all units (right).

K. Le Son et al. / Reliability Engineering and System Safety 149 (2016) 76–87

83

Fig. 8. Parameter estimations with respect to the iteration number of the Gibbs algorithm.

Table 2 Parameter estimations and their variance.

Estimated parameter Variance

b α

b β

b b

σb2

0.8655 0.0006

82.1799 0.6777

1.0070 0.0001

1.7861 0.0001

considered. This criteria is presented as follows: ( K X e  dk =13  1; dk r 0 S¼ k ¼ 1; …; K Sk where Sk ¼ dk 4 0 edk =10 1; k¼1

ð20Þ

and dk ¼ estimated RUL of component kreal RUL of componentk ð21Þ Other criteria are also considered such as the mean square error (MSE) defined as follows: MSE ¼

K 1X ðd Þ2 Kk¼1 k

ð22Þ

In order to evaluate the performance of our method these criteria are applied to the estimated RUL of each component. The errors dk of each component are presented in Fig. 9. For most of the components the estimated RUL is smaller than the real RUL in

the training data set (dk Z 0). This result highlights the efficiency of our method and shows the interest of its use in the framework of a preventive maintenance where it is necessary to plan maintenance actions before the failure. By using the Gamma process with Gaussian noise model to model our degradation indicator, we obtain respectively the results of predicted score S¼4170 and the Mean Squared Error MSE¼ 864. These results improve the scores obtained for the same indicator in [9] where a Wiener process is used to model the deterioration, S ¼5520 and MSE ¼819. The latter highlights the interest of our deterioration model for the RUL estimation on PHM 2008 challenge data. We can also compare our results with best obtained results of the PHM Challenge [31,32] and the other models [9] in Table 3. Table 3 shows that our model gives the best result in terms of score criterion. Generally, in this case the probabilistic model show better performences than the non-probabilistic models. By comparing different modeling approach on the same degradation indicator, the efficiency of our filtering approach for RUL estimation by using Gibbs sampling is highlighted. These results point out the efficiency of our deterioration model and also the probabilistic RUL estimation method. After having studied the impact of probabilistic approach with filtering method on a study case, we try to explore the use and interest of RUL estimation in industrial applications. One of the

84

K. Le Son et al. / Reliability Engineering and System Safety 149 (2016) 76–87

with cost Cc. We suppose that after replacements the system is as good as new and the inspection cost is much lower than the replacement costs. The corrective replacement is carried out as soon as a failure is detected. The maintenance policies proposed in this paper differ in the way they handle the preventive replacements. These policies are developed in the following sections. To evaluate the maintenance policies, we focus on the asymptotic expected maintenance cost per unit over an infinite time span considered as a cost criterion: C 1 ¼ lim

CðtÞ ; t

t-1

where C(t) is the cumulated maintenance cost at time t where CðtÞ ¼ C i N i ðtÞ þC p N p ðtÞ þ C c N c ðtÞ þ C d dd ðtÞ;

Fig. 9. Error estimations dk for the testing data set.

Table 3 Evaluation of the different prognostic approaches. Prognostic model Type Probabilistic

Non-probabilistic

Criterion Approach Lifetime distribution [9] Wiener process [9] Gamma process with noise

S 9450 5575 4107

MSE 1011 823 864

Similarity-based [9] Similarity-based [31] Neuron network [32]

6690 5636 X

809 X 984

most common use of the RUL is in the framework of preventive maintenance. This problem is addressed in the next section.

with Np(t) the number of preventive replacements before t, Nc(t) the number of corrective replacements before t, dd(t) the cumulative unavailability duration of the system before t and Ni(t) the  number of inspections before t. Note that Ni ðtÞ ¼ t=T where ½x denotes the integer part of the real number x. In the following paragraphs the maintenance policies are presented and evaluated. 4.1. State-based maintenance model: (T,M) policy The system at the inspection time Tk is denoted by Xk. In the framework of the state-based maintenance policy a fixed threshold M o L is defined and at each inspection time Tk:

 if M o X k oL, the system is preventively replaced with a cost Cp.  if X k o M, the decision is postponed until the next inspection T k þ 1.

 if X k Z L, the system has already failed, a corrective replacement is carried out with a cost Cc.

4. Remaining useful life and maintenance In this section, to have an idea about how the RUL can be used in the framework of the maintenance in this section, two maintenance policies based on the remaining useful life of the system (or component) are considered and these policies are compared with a condition-based policy where only the state of the system is taken into account. In this section, we try to highlight the interest of a probabilistic prognosis method where in difference with the non-probabilistic method we have access to the RUL distribution. For this comparison one of the two RUL-based maintenance policies presented in this section is based only on the mean value of the RUL and the other policy is based on the quantile of the RUL. The comparison of these two policies can emphasize the probabilistic prognosis methods. In the policies under consideration in this section, the maintenance decision is taken on the estimated state after filtering and not on the noisy observation. Using these estimated states let us to proceed prognosis as if we had a non-noisy monotone increasing deterioration (Gamma process) and therefore use the gamma process properties to estimate the failure time. These calculations are possible and reliable because we have used a robust and efficient filtering method called Gibbs sampler. Throughout this section, it is supposed that the system is inspected periodically at inspection times T 1 , T 2 ; … where T k ¼ kT with k A N and T A R is the inspection interval. The deterioration state is hidden, hence at each inspection time the hidden degradation state is estimated and an inspection cost Ci is incurred. To simplify notations we shall denote Xk the filtered state at the inspection time Tk that we consider as the real state at this time. At each inspection time two maintenance operations are possible: a preventive replacement with cost Cp and a corrective replacement

The inspection period T and the preventive threshold M are the two decision variables to be optimized. The optimal values of inspection period Topt and of preventive threshold Mopt are obtained as follows: C 1 ðT opt ; M opt Þ ¼ min C 1 ðT; MÞ T;M

where T 4 0 and 0 o M o L. 4.2. E(RUL)-based maintenance model: (T,RULmin) policy Recall that the RUL of the system at the inspection time Tk is denoted by RULðT k Þ and the expectation of RULðT k Þ is defined as follows: Z dF RULðT k Þ ðhÞ h: dh ð23Þ EðRULðT k ÞÞ ¼ dh In the framework of the E(RUL)-based maintenance policy a fixed threshold RULmin is defined and at each inspection time Tk:

 if X k oL and EðRULðT k ÞÞ rRULmin , the system is preventively replaced with a cost Cp.

 if X k o L and EðRULðT k ÞÞ 4 RULmin , the estimated mean RUL is still 

higher than the threshold RULmin, the decision is postponed until the next inspection T k þ 1 . if X k Z L, the system has already failed, a corrective replacement is carried out with a cost Cc.

The inspection period T and the preventive threshold RULmin are the two decision variables to be optimized. The optimal values of inspection period Topt and of preventive threshold RULmin;opt are

K. Le Son et al. / Reliability Engineering and System Safety 149 (2016) 76–87

85

obtained as follows:

4.4. Numerical examples

C 1 ðT opt ; RULopt Þ ¼ min C 1 ðT; RULmin Þ

The aim of this study is to compare the performances of the maintenance policies. For the numerical example, the maintenance costs are chosen as C i ¼ 2, C p ¼ 50, C c ¼ 100, C d ¼ 150. In the order to compare the performance of RUL-based maintenance policies, we consider two cases of parameters for the nonhomogenous gamma degradation process:

T;RULmin

where T 4 0 and 0 o RULmin o MTTF.

4.3. FRUL-based maintenance mode: (T,Q) policy Recall that F RULðT k Þ and Xk are respectively the cumulative distribution function of remaining useful life RULðT k Þ and the degradation state at the inspection time T k . In the framework of the FRUL-based maintenance policy a fixed percentile Q, 0 oQ o 1, of the RUL distribution function is defined and at each inspection time Tk:

 If X kT oL and PðRULðT k Þ o TÞ 4 Q , the system is preventively replaced with a cost Cp.

 If X kT o L and PðRULðkTÞ o TÞ o Q , the decision is postponed until the next inspection.

 If X kT Z L, a corrective replacement is carried out with a cost Cc. The optimal values of T and Q are obtained by the following expression: C 1 ðT opt ; Q opt Þ ¼ minfC 1 ðT; Q Þg T;Q

where, 0 o Q o 1 and T 40. Table 4 Optimal mean cost per time unit for three different policies in two cases Policy Case 1

Case 2

(T,RULmin)

(T,Q)

(T,M)

C1 opt

9.9417

9.9398

10.3402

Optimal variables

(1.6, 2.2)

(1.6, 0.19)

(1.3, 12.1)

C1 opt

6.0085

6.0083

6.3706

Optimal variables

(4.8, 5.7)

(4.8, 0.18)

(2.4, 11.4)

 Case 1: α ¼ 1, β ¼ 1, b¼1.5, the failure threshold L¼ 20.  Case 2: α ¼ 1, β ¼ 2, b¼1.5, the failure threshold L¼ 20. In this paper, the gamma distributed deterioration increments are in the presence of a Gaussian noise with parameter σ ¼ 2. The optimal costs and the optimal parameters of the three policies in the two cases are presented in Table 4. It can be noticed that the optimal inspection periods are the same for both prognosis policies and the optimal long-run expected maintenance costs are very close. The state-based policy gives the least attractive results. Hence the prognosis policies seem to be more cost-efficient. At this point we can highlight the interest of using the RUL in the framework of a maintenance policy. It seems that the choice of the type of prognosis maintenance policy is not crucial. To have a better conclusion and comparison, in the following paragraphs other performances of the two prognosis-based policies are studied. In order to have a better comparison of the policies, the impact of the maintenance costs variations on the total maintenance cost should be studied. A comparison of optimal expected cost and inspection period by varying the maintenance costs is carried out as follows:

 Variable inspection cost: C i A ½2; 20, Cp ¼50, Cc ¼100, Cd ¼150.  Variable preventive replacement cost: Ci ¼ 2, C p A ½5; 100, Cc ¼100, Cd ¼150.

 Variable inactivity cost: Ci ¼2, Cp ¼50, Cc ¼100, C d A ½10; 200. Figs. 10 and 11 show the optimal inspection time and the optimal long-run expected maintenance cost by varying the

Fig. 10. Evaluation of the optimal inspection period in case 1.

86

K. Le Son et al. / Reliability Engineering and System Safety 149 (2016) 76–87

Fig. 11. Evaluation of the optimal long-run expected maintenance cost rate in case 1.

maintenance costs in case 1. It can be noticed that the ðT; RULmin Þ policy and (T,Q) policy have a similar behavior but the ðT; RULmin Þ policy is more sensitive to variations of the preventive cost. The results are similar in case 2 and for this reason we have skipped the figures in this case. Comparing the computation times, the (T,Q) policy is less time-consuming. However, in the case of the proposed parameters both policies have the same profile. 4.5. Performance assessment In Fig. 12 the optimal costs for different values of Q in the (T,Q) policy (up), for different values of RULmin in the ðT; RULmin Þ policy (down) in case 2 are depicted. The (T,Q) policy is more stable around the optimal value. In Fig. 13 the optimal maintenance costs at RULmin;opt and Qopt for different values of the inspection interval T are evaluated. From this figure, we can highlight the stability of the (T,Q) policy to the variations of T around Topt. The (T,Q) policy is more robust because the use of RUL distribution in the maintenance policy makes this latter more flexible than the policy based on a unique estimated RUL value. When for example we are not sure to be able to inspect exactly at the optimal inspection time, the maintenance decision should be made on the basis of the quantile instead of EðRULðT k ÞÞ. In sum, both maintenance policies (T,Q) and ðT; RULÞ have better performances than the traditional control-limit conditionbased maintenance policy (T,M). In this sense there are both more advantageous than (T,M) policy. Moreover, the two RUL based policies lead to similar optimal cost values and this latter remains

Fig. 12. Cost variation for different values of Q in the (T,Q) policy (up), for different values of RULmin in the ðT; RULmin Þ policy (down) in case 2.

true for different maintenance cost units. Therefore the advantage of one policy over the other is debatable. However, a small advantage of the (T,Q) can be pointed out as follows:

K. Le Son et al. / Reliability Engineering and System Safety 149 (2016) 76–87

Fig. 13. The optimal cost variation at RULmin;opt and Qopt for different T values in case 2.

 The (T,Q) policy optimisation is slightly less time-consuming.  In the (T,Q) policy for a large number of non-optimal values T, the long-run average maintenance cost is close to the optimal cost value. This result remains also true for the non-optimal values Q. This property ensures that even if the policy is applied with non-optimal parameters chosen reasonably in advance the long-run average maintenance cost remains close to its optimal value. Therefore the policy implementation is less sensitive to the optimisation procedure which makes its application to real problems much easier.

5. Conclusion This paper shows a probabilistic model for prognostic of a deteriorating system modeled by a noisy degradation Gamma process. This model is suitable for a degradation phenomenon not directly observable. A filtering of deterioration model is proposed using Gibbs sampling technique and the impact of this method is pointed out. Based on the probabilistic properties of the proposed model the remaining useful life of the system is estimated and its properties are analyzed. In the case of Gamma deterioration with noise, we can consider the proposed Gibbs sampling technique as a Benchmark for stochastic filtering in order to evaluate the other complicated filtering methods on simulating data or on our study case (2008 PHM conference challenge data). The RUL estimation method is applied to the real data and its performances are analyzed. Eventually to highlight the utility of the RUL estimation taken into account the maintenance decision-making, two prognosis based maintenance policies are proposed.

References [1] Abdel-Hameed M. Inspection and maintenance policies of devices subject to deterioration. Adv Appl Probab 1987:917–31. [2] Van Noortwijk J. A survey of the application of gamma processes in maintenance. Reliab Eng Syst Saf 2009;94(1):2–21.

87

[3] Lawless J, Crowder M. Covariates and random effects in a gamma process model with application to degradation and failure. Lifetime Data Anal 2004;10 (3):213–27. [4] Singpurwalla N. Gamma processes and their generalizations: an overview. Eng Probab Des Maint Flood Prot 1997:67–75. [5] Lee M, Whitmore G. Threshold regression for survival analysis: modeling event times by a stochastic process reaching a boundary. Stat Sci 2006; 21(4):501–13. [6] Schirru A, Pampuri S, De Nicolao G. Particle filtering of hidden Gamma processes for robust predictive maintenance in semiconductor manufacturing, in: 2010 IEEE Conference on Automation Science and Engineering (CASE), IEEE; 2010. p. 51–6. [7] Park C, Padgett W. Accelerated degradation models for failure based on geometric Brownian motion and gamma processes. Lifetime Data Anal 2005; 11(4):511–27. [8] Saxena A, Goebel K, Simon D, Eklund N. Damage propagation modeling for aircraft engine run-to-failure simulation. In: Proceedings of 2008 International Conference on Prognostics and Health Management, 2008. [9] Le Son K, Fouladirad M, Barros A, Levrat E, Iung B. Remaining useful life estimation based on stochastic deterioration models: a comparative study. Reliab Eng Syst Saf 2013;112:165–75. [10] Kallen M, van Noortwijk J. Optimal maintenance decisions under imperfect inspection. Reliab Eng Syst Saf 2005;90(2–3):177–85. [11] Lu D, Pandey MD, Xie W. An efficient method for the estimation of parameters of stochastic gamma process from noisy degradation measurements. Proc IMechE Part O: J Risk Reliab 2013;227:1–9. [12] Jazwinski A. Stochastic processes and filtering theory, vol. 63, Academic Pr; 1970. [13] Baruah P, Chinnam R. Hmms for diagnostics and prognostics in machining processes. Int J Prod Res 2005;43(6):1275–93. [14] Batzel T, Swanson D. Prognostic health management of aircraft power generators. IEEE Trans Aerosp Electron Syst 2009;45(2):473–82. [15] Cappé O, Godsill SJ, Moulines E. An overview of existing methods and recent advances in sequential monte carlo. Proc IEEE 2007;95(5):899–924. [16] Doucet A, De Freitas N, Gordon N. Sequential Monte Carlo methods in practice. New York: Springer Verlag; 2001. [17] Liu J. Monte Carlo strategies in scientific computing. New York: Springer Verlag; 2008. [18] Zhang Y, Chen M-H, Ibrahim JG, Zeng D, Chen Q, Pan Z, et al. Bayesian gamma frailty models for survival data with semi-competing risks and treatment switching. Lifetime Data Anal 2014;20(1):76–105. [19] Fox J-P. Stochastic em for estimating the parameters of a multilevel irt model. Br J Math Stat Psychol 2003;56(1):65–81. [20] Le Son K, Fouladirad M, Barros A. Remaining useful life estimation on the nonhomogenous gamma with noise deterioration based on gibbs filtering: A case study. in: 2012 IEEE international conference on prognostics and health management, Denver, Colorado; 2012. [21] Le Son K, Fouladirad M, Barros A. Deterioration model filtering by Gibbs algorithm and RUL estimation by using a gamma process with noise. in: 8th Safeprocess 2012, Mexico city, Mexico; 2012. [22] Casella G, George E. Explaining the Gibbs sampler. Am Stat 1992;46(3): 167–174. [23] Tierney L. Markov chains for exploring posterior distributions. Ann Stat 1994;22(4):1701–28. [24] Ripley B. Stochastic simulation, vol. 183, Wiley Online Library; 1987. [25] Brooks S, Roberts G. Convergence assessment techniques for Markov chain monte carlo. Stat Comput 1998;8(4):319–35. [26] Gelman A, Rubin D. Inference from iterative simulation using multiple sequences. Stat Sci 1992;7(4):457–72. [27] Nielsen F. The stochastic EM algorithm: estimation and asymptotic results. Bernoulli 2000;6(3):457–89. [28] Zhan H, Chen X, Xu S. A stochastic expectation and maximization algorithm for detecting quantitative trait-associated genes. Bioinformatics 2011; 27(1):63. [29] Lagarias J, Reeds J, Wright M, Wright P. Convergence properties of the Nelder– Mead simplex method in low dimensions. SIAM J Optim 1999;9(1):112–47. [30] Le Son K, Fouladirad M, Barros A, Levrat E, Iung B. Remaining useful life estimation based on probabilistic model. In: 17th ISSAT Int Conf Reliab Qual Des, Vancouver, Canada; 2011. [31] Wang T, Yu J, Siegel D, Lee J. A similarity-based prognostics approach for remaining useful life estimation of engineered systems. In: PHM 2008 international conference on prognostics and health management; 2008, 2008, p. 1–6. [32] Heimes F. Recurrent neural networks for remaining useful life estimation. In: PHM 2008 international conference on prognostics and health management, 2008; 2008. p. 1–6.