Remaining useful life estimation based on stochastic deterioration models: A comparative study

Reliability Engineering and System Safety 112 (2013) 165–175

Contents lists available at SciVerse ScienceDirect

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

Remaining useful life estimation based on stochastic deterioration models: A comparative study Khanh Le Son a, Mitra Fouladirad a,n, Anne Barros a, Eric Levrat b, Benoˆıt Iung b a b

Universite´ de Technologie de Troyes, Institut Charles Delaunay, CNRS UMR 6279 STMR,12 rue Marie Curie, 10010 Troyes, France, Universite´ Henri Poincare´, CRAN, CNRS UMR 7039, BP 239, 54506 Vandoeuvre, France

a r t i c l e i n f o

abstract

Article history: Received 20 January 2012 Received in revised form 23 November 2012 Accepted 28 November 2012 Available online 7 December 2012

Prognostic of system lifetime is a basic requirement for condition-based maintenance in many application domains where safety, reliability, and availability are considered of first importance. This paper presents a probabilistic method for prognostic applied to the 2008 PHM Conference Challenge data. A stochastic process (Wiener process) combined with a data analysis method (Principal Component Analysis) is proposed to model the deterioration of the components and to estimate the RUL on a case study. The advantages of our probabilistic approach are pointed out and a comparison with existing results on the same data is made. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Deterioration modelling Principal component analysis Stochastic process Wiener process with drift Prognostic Remaining useful life time estimation

1. Introduction Prognostic is a major scientific challenge for industrial implementation of maintenance strategies in which the remaining useful life estimation is an important task. For environmental, economical and operational purposes, prognostic and remaining useful lifetime (RUL) prediction arouse a big interest, (see [2]) and is developed in many domains as aerospace industry [29,35] electronics [6], medicine [1], nuclear [4,5], finance [8], weather forecast [17], etc. In the framework of prognostic and health management (PHM), many prognostic techniques exist and they are basely classified into three principal classes: data-driven approaches, model-based approaches, experience-based approaches [7,22]. However, in this paper, we consider a mixture of these three classes and divide mainly the tools for prognostic methods in two groups: non-probabilistic methods and probabilistic methods. The objective is to compare the efficiency of RUL prediction model on these two methods types with a study case (2008 Prognostic Health Management (PHM) Challenge data [29]). In non-probabilistic methods the deterioration phenomenon is not random and at most the observations on the deterioration can be noisy. The prognosis in these cases is based on the mean deterioration. Many applications of non-probabilistic

n

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

0951-8320/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ress.2012.11.022

methods on prognostics are realised. For example in [13], the authors used the multiple linear regression analysis to estimate RUL on HUMS (Health and Usage Monitoring Systems) conditions indicators extracted from vibration data of bearings. In [24], based on high-fidelity physics-based model, a prognostic RUL estimation tool is constructed for VDPs (Variable Displacement Pumps). In the particular case of PHM challenge data, the neural networks model is applied to RUL prediction of components in [26,14]. In [34] based on training data, the authors constructed a similaritybased model and applied to testing data in order to calculate the RUL of components. These above non-probabilistic methods are popularly used in prognostic applications. On the contrary, with the probabilistic methods, the deterioration phenomenon is considered to be random and with stochastic tools we deal with its random behaviour. In this case the prognosis is based on the future behaviour of the stochastic deterioration process and can give results in terms of probabilities. Recently many works on probabilistic prognosis methods have been carried out. In[30], a review on probabilistic prognostic methods is given. In these latter the RUL estimation is based on probabilistic tools and the deterioration phenomenon is modelled by a stochastic process. The Levy processes (specially the Wiener process and Gamma process) are the most common stochastic processes used for deterioration modelling. For 2008 IEEE PHM challenge, authors had to address the problem of RUL estimation for a given data set. Only nonprobabilistic methods have been proposed and the best results

166

K. Le Son et al. / Reliability Engineering and System Safety 112 (2013) 165–175

have been obtained by using respectively neural networks [26,14], Bayesian networks [10], and regression [34]. These results are considered here as a reference and a starting point. Our main objective is to propose a probabilistic prognostic model and to assess its performance. First Principal Component Analysis is used for the data analysis step in order to construct a degradation indicator. Due to the general trend of this indicator, the Wiener process has been selected. Actually, the use of Wiener process for the deterioration modelling is very popular where the deterioration process is noisy and varies non-monotonically. In [11,9], authors studied the statistical properties of the failure time in the case of a Wiener deterioration model and their results have been applied in reliability and lifetime analysis widely since 1970s. In [12], the authors applied Brownian motion with drift to a variable-stress accelerated life testing experiment. In [39,40] the impact of measurement errors on the Wiener degradation model of self-regulating heating cables is analysed. Authors in [23] also focused on the stopping time (failure time) of Wiener deterioration models. Recently, in order to have more realistic deterioration models, Wiener process with non-linear drift has attracted many interest, see for example [38,32,30,31]. As it is mentioned in [30] there are increasingly interesting results on probabilistic residual life estimation problem and some of them such as [36,37] focus on residual life estimation problem in the framework of Wiener deterioration model is developed. The remainder of the paper is organised as follows. In Section 2 a preliminary presentation of the RUL and the IEEE PHM challenge is carried out. In Section 3, a degradation indicator based on Principal Component Analysis is proposed. The prognostic method is applied with the new degradation indicator by using stochastic processes and the RUL estimation method is presented in Section 4. Finally, conclusions of the results as well as the further works are given in Section 5.

2. Layout of the study

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. 2.3. Assessment of RUL estimation quality Assessment of each method performance is based on a function of the error in predicting the number of remaining time cycles. The score function S is asymmetric and exponential, such that late predictions are penalised more heavily than early predictions and it is defined as follows: S¼

N X

Si

ð1Þ

i¼1

where N is the number of units (N ¼218), the prediction error di and the penalty score Si for unit i are computed as follows: di ¼ estimated RULi actual RULi ( Si ¼

edi =13 1,

di r 0

edi =10 1,

di 4 0

i ¼ 1, . . . ,N

ð2Þ ð3Þ

This score is originally used in 2008 PHM Conference Challenge. This score is asymmetrical and penalise more the late failure estimations where the error di presented in Eq. (2) is positive. It seems sensible to give advantage to methods which are more conservative and propose a failure prediction before the effective failure. Some papers use another criterion to compare the accuracy of methods, the root squared error (RSE) is defined as follows: vffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N uX 2 RSE ¼ t ð4Þ di i¼1

2.1. RUL definition Generally speaking, the RUL is defined at any time t as the remaining lifetime of a unit given that it is running at time t and given all the available information related to the unit at time t. Depending on the method used for estimating this quantity and depending on the available information, the formalisation of the RUL definition can change. In the following we propose a formalism connected to the probabilistic models we used. 2.2. Data presentation 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 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 [29]). 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

And finally some papers use another criterion the mean squared error of the prediction (MSE) defined as follows: MSE ¼

2 N X di 218 i¼1

ð5Þ

The lower is the score the better is the proposed RUL estimation algorithm. In the PHM conference challenge, the following papers obtained the best results. Wang in [34], obtained the total score of S¼5636 with a similarity-based approach for estimating the RUL. Felix in [14] with a recurrent neural networks method achieved the prediction error of 519.8. In the paper [26] an artificial neural networks mixed with a Kalman filter is used and the obtained mean square error is MSE¼984.

3. Degradation indicator 3.1. Lifetime distribution As a starting point before using stochastic processes, we apply the most basic probabilistic method for the RUL prediction. This method consists of a simple lifetime approach. This means that we only consider the failure times of the units, without taking into account the progressive degradation phenomena that leads from the initial state to the failure. Failure times are considered as independent and identically distributed random variables. Let T be the failure time (lifetime) of a unit. While treating with failure times, the average RUL at time t of each unit can be considered as the mean time to failure

K. Le Son et al. / Reliability Engineering and System Safety 112 (2013) 165–175

40

1

35

0.9

167

0.8

Distribution function

30 25 20 15

0.7 0.6 0.5 0.4 0.3

10 0.2 5 0 100

Accumulated distribution function Fitted Weibull distribution function

0.1 150 150

200

250

300

350

200

250

400

Cycles Fig. 1. Histogram of failure times for 218 units.

l

To estimate the Weibull distribution parameters, the empirical distribution function obtained from the failure times of the training data set is fitted to the cumulative distribution function of a three parameters Weibull distribution see Fig. 2. By a nonlinear regression the estimated parameters are obtained respectively ðm, l, yÞ ¼ ð2:606, 119:2, 95:02Þ and based on these estimated parameters the RUL is estimated by using the mean residual life given by Eq. (6) on the testing data set. The RUL of each unit is estimated and is taken into the score criterion, thus we get the score as S ¼9450. As expected, this score is not very convincing in comparison with other non-probabilistic methods. Therefore, the use of a more sophisticated prognostic model is necessary to improve these results and in the framework of probabilistic approaches it seems natural to propose a stochastic process to model the degradation phenomena. 3.2. Preliminary data analysis When using stochastic processes, it is first necessary to identify in the data a ‘‘degradation indicator’’, i.e. many realisations of a scalar indicator over a given time horizon (or a

350

Fig. 2. The empirical distribution function obtained from the failure times of the training data set and its fitted cumulative distribution function.

at time t considering that the unit is not failed before t presented as follows: Z 1 1 EðTt9T 4 tÞ ¼ RðxÞ dxt ð6Þ RðtÞ t

100 80 60

OP3

where RðxÞ ¼ PðT 4 xÞ is the reliability function. The Weibull distribution is one of the traditional distributions associated to failure times. According to the histogram of the failure times of 218 units presented in Fig. 1, it seems reasonable to propose a Weibull distribution to model the life time distribution, see for example [10] for more details. Since the failure time does not start from zero it can be sensible to consider a non-zero location parameter for the Weibull distribution. Therefore for the lifetime distribution modelling and the mean residual life calculation a three parameters Weibull distribution is considered. The Weibull probability density function with three parameters ðm, l, yÞ where m is the shape parameter, l the scale parameter and y location parameter is defined as follows:   m  xy m f ðx, m, l, yÞ ¼ ml ðxyÞm1 exp  ð7Þ

300

Cycle

40 20

0 60 40

OP1

20 0

1

0.8

0.6

0.4

0.2

0

OP2

Fig. 3. Three operational setting variables of all units.

multivariate indicator with reasonable size) representative of an evolution from a running state to a failed state. After a first study of the training data the following features were pointed out:

 The three measurements connected to operational settings can

 

be represented by six clusters in a 3-dimensional space and are illustrated in Fig. 3, where the variables OP1, OP2, OP3 correspond to the measures of the three sensors devoted to operational settings. The six observed dots are actually six highly concentrated clusters. Hereafter it is supposed that these clusters indicate six discrete operating conditions of the system. Hence, each cycle of a unit can be labeled by a regime ID from 1 to 6, replacing the original three variables of operational settings. A degradation indicator cannot be directly deduced from the 21 sensor paths (see Fig. 4), If we select only the measurements corresponding to the same regime ID, some of the 21 sensor measurements connected to units state exhibit prominent trend along the life of a unit. In [34], for each mode, the plot of the collected data from one sensor is used to show the trend of the system’s degradation (as an example refer to Fig. 5 given by [34]). Some sensors do not show a clear trend due to high noise or their low sensitivity to the degradation. Including them in the analysis

168

K. Le Son et al. / Reliability Engineering and System Safety 112 (2013) 165–175

PCA. Then, the degradation indicator is defined at each time t, for each unit, by the distance between the projection of the seven sensor measures at time t in the failure space, and the failure place. However, the structure of the failure space, the behaviour of the obtained degradation indicator and the previous analysis about the impact of the six operational modes lead us to adopt a more complex strategy: we consider a different failure space for each mode by dividing the preceding sub-data set of the failure times into six groups (looking at the operational mode at failure time). Then the degradation indicator at the each time t, for each unit, is defined as the distance between the projection in the failure space of current unit’s mode, and the failure place of the current mode. For one selected mode, there is only one failure place corresponding to the barycenter of the failure projections set. This indicator construction method is detailed in the following section.

may lower the accuracy of prediction. According to [34], using a subset of the available sensors can help to improve the accuracy of RUL estimation, therefore they propose to use a selection of seven sensors (2, 3, 4, 7, 11, 12 and 15).

3.3. Main idea The goal is to identify:

 a set of paths of one indicator which can be considered as the 

realisations of the same stochastic process to make statistical inference, a real trend of this indicator that can be interpreted as a convergence to a failed state and that allows probabilistic RUL prediction.

3.4. Degradation indicator construction

A first study of the data presented in the previous Section shows that this scalar should be found using the selection of seven sensor measures proposed by Wang [34] but none of the winners of the 2008 PHM challenge propose to construct such an indicator. Then our idea is to identify a ‘‘failure place’’ in a ‘‘failure space’’ of dimension 2 with a Principal Component Analysis (PCA) applied on the sub-data set of the failure times, see [18] for more details on

3.4.1. Study of the space of the failure times with PCA First, we show the results obtained for the failure space when applying PCA without taking into account the operational modes. The PCA is explained very clearly in [18], and an example of using PCA technique for prognostic is given in [20]. We consider the data subset of the failure times for 218 units. Let X be the original data set where its rows indicate the information on individuals (at the failure time) and its columns indicate the information on variables (seven sensors). For the PCA corresponding to the chosen data set, the principal component’s coefficients and variances are the eigenvectors and eigenvalues of the covariance matrix of the centering matrix of X. The criterion to choose the principal component is based on the histogram of the variances of principal component’s. It can be noticed in Fig. 6 (left), that the cumulative percentage of the total variation of the first two principal components is more than 98%. The Fig. 6 (right) shows that the projections of failure times measurements on the two dimensional principal component space are divided into six different groups. This result makes it difficult to define one failure place and the distance to it. However we have tried it, defining the indicator as the distance between the current projection of the seven sensors measures at time t and the closest failure place at this time or the barycenter of the six failure places. But the resulting paths of the scalar indicator show no trend and no evolution that could be used for probabilistic RUL prediction. To overcome this problem we separate the operational modes of the failure times to build six data sets and we apply a PCA to data corresponding to each operational mode. Therefore, we obtain 6 different principal component spaces noted P 1 , . . . ,P6 hereafter. In Table 1 1 the total

Measurements of sensor 1 and 2 in the terms of times for unit 1

650

550

500

450 SM2 SM1

0

50

100

150

200

250

T Fig. 4. Measurements of sensor 1 and 2 for the unit 1.

Measurements of sensor 2 in mode1

644.5

2214.5

644

2214

643.5

2213.5

643

2213

SM8

400

SM2

Sensor measurement

600

642.5

2212.5

642

2212

641.5

2211.5

641 −350

−300

−250

−200

−150

Cycle

−100

−50

0

Measurements of sensor 8 in mode 4

2211 −400

−300

−200

−100

Cycle

Fig. 5. Measurements of sensor 2 in mode 1 (left), measurements of sensor 8 in mode 4 (right).

0

K. Le Son et al. / Reliability Engineering and System Safety 112 (2013) 165–175

Projections of failure places on the 2-D space of the PCs

Percentages of the PCs −4900

100 90

−5000

70

−5100

60

−5200

PC2

Percentage (%)

80

50 40

−5300

30

−5400

20

−5500

10 0

169

0

5

10

15

20

25

−5600

9000

9500

10000

PC1

PCs

Fig. 6. Principal component analysis (left), the two dimensional principal component space (right).

Failure projections for the operational mode 1

Table 1 The total variation of the three first components for each of the sixth operational mode.

−1922

Mode 1

Mode 2

Mode 3

Mode 4

Mode 5

Mode 6

−1924

60.85 38.04 0.66

72.64 26.75 0.28

61.45 37.85 0.34

54.41 44.55 0.56

58.95 40.07 0.41

79.65 19.11 0.77

−1926

3.4.2. Degradation indicator defined as a distance The degradation indicator can then be calculated at each time, for each unit, as the euclidean distance between the seven selected sensors measurements projected in the failure space of the current unit’s mode, and the unique failure place of the current mode. In order to improve the accuracy of degradation indicator, the distance of each unique failure space is divided by the dispersion. The dispersion of each unique failure space is defined as the standard deviation of the barycenter of the failure places cluster in each mode see Figs. 7 and 8. The construction of the indicator is formalised as follows:

 let Nk be the number of individuals at the failure time in mode k,  let ðx k ,y k Þ be the barycenter of the projected measurement of the failure places in mode k on Pk, k ¼ 1, . . . ,6,

 let ðxki ,yki Þ, be the ith projected failure places in mode k on Pk, 

i ¼ 1, . . . ,Nk , k ¼ 1, . . . ,6, let ðxkij ,ykij Þ be the projected measurement at the arbitrary operating time j of the component i in mode k on Pk, k ¼ 1, . . . ,6, i ¼ 1, . . . ,Nk .

We proceed considering the following distance for unit i, at time cycle j, in mode k: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxkij x k Þ2 þðykij y k Þ2 k ð8Þ Di,j ¼ Disperk

−1930 −1932 −1934 −1936 −1938 −1940 915

920

925

930

935

PC1 Fig. 7. The projections of the failure states on the two principal axes for the first mode.

−1905

Projections of failure on the 2-D space of PCs in mode k (Pk) Pi,jk

(Projection at the cycle j of unit i on the failure space in mode k)

−1910

k

−1915

PC2

variation of the three first components for each of the sixth operational mode are given. The percentages of histogram for the first two principal components in Table 1 shows that we can use the two-dimensional principal component space to build a degradation indicator for each unit in each mode. Moreover, the high concentration of the positions of different component failure states on the two principal components space illustrated by Fig. 7 leads us to consider only one failure place for each operational mode. We can consider the barycenter of the cluster of failure states on the two principal components space as the unique failure place.

Barycenter of failure places

−1928

PC2

PC1 PC2 PC3

−1920

Di,j =

¯ ) distance ( Pi,k j ; P k Disperk

−1920 Barycenter of the ¯k failure projectionsP

−1925

−1930

−1935 Failure space in mode k −1940 900

905

910

915

920

925

930

935

PC1 Fig. 8. Degradation indicators definition.

where Disperk the dispersion of the failure places in mode k is defined asvfollows: ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u Nk u 1 X t ððxk x k Þ2 þ ðyki y k Þ2 Þ ð9Þ Disperk ¼ Nk 1 i ¼ 1 i

170

K. Le Son et al. / Reliability Engineering and System Safety 112 (2013) 165–175

Degradation in dicators of all units

Degradation in dicator of component 1

1.25

2.5

Degradation in dicator

3

Degradation in dicator

1.5

1

0.75 0.5

0.25 0

2 1.5 1 0.5

0

50

100

150

200

0

250

0

50

100

150

Cycle

200

250

300

350

400

Cycle

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

Unit 2

2

1.5

1.5

1

1

0.5

0.5

0

0 0

50

100

150

200

Unit 188

2

Unit 100

2

0

50

100

150

200

250

Unit 218

1.4 1.2

1.5

1 0.8

1 0.6 0.4

0.5

0.2 0

0

50

100

150

200

250

0

0

20

40

60

80

100

120

140

Fig. 10. Some examples of degradation paths.

For each unit, the degradation indicator is sequentially calculated by applying Eq. (8) at each time cycle j and looking at the current mode k at this time. Then the degradation indicator over ½1, . . . ,ri  for ri cycles can be formalised by the vector kðr i Þ ,Dkð2Þ , . . . ,Di,r  Di,½1:ri  ¼ ½Dkð1Þ i,1 i,2 i

the degradation indicator and the stochastic process not only takes into account the non-linear tendency but also the presence of uncertainties. The two methods are compared and their performances are analysed.

ð10Þ

where k(j) denotes the mode at time j. As an example, in Fig. 9 (left), the degradation indicator for component 1 is depicted and in Fig. 9 (right) the degradation indicators of all components are illustrated. It can be noticed in Fig. 9 (right) and Fig. 10 that there is an important variation between the degradation indicators of different units and there is a general trend: non-linear and decreasing (at the beginning cycles, the indicator value takes around 1 and at the end cycles, at the failure times, it tends to zero). Therefore we model the trend of this indicator by considering these characteristics and then we propose the criterion of RUL estimation. In the following section, the degradation indicator is modelled once with a nonlinear deterministic function and once by a stochastic process. The non-linear function takes into account the non-linear trend of

4. RUL estimation In this section, based on the built indicator, a RUL estimation method is proposed using a stochastic process to model the degradation. First, the stochastic prediction model based on Wiener process is presented and its parameters are estimated with the training data set. Then the RUL is estimated for each unit of the testing data set and the score of our method is calculated. At last, in Section 4.2, a comparison with a non-probabilistic method allied on PHM 2008 conference challenge data is carried out. For the comparison, we did not select the Bayesian network ones because the a priori knowledge that is necessary to get good performance indicators seems to us really unrealistic. Then, neuron network and similarity-based models constitute the two main classes of models that are relevant to compare with

K. Le Son et al. / Reliability Engineering and System Safety 112 (2013) 165–175

stochastic ones. Considering neuron networks, the comparison is really quick, since we can only look at the obtained performance indicators. So we have mentioned it quickly in the paper to show that the performance of our model is good compared to the one of a neuron network. On the opposite the comparison with similarity-based models is much more interesting since there is a common step (the construction of the degradation indicator) and a non-common step (the choice of the model to fit with the degradation indicator). Hence we can separate the performance due to the construction of the degradation indicator and the performance due to the choice of a degradation model (nonstochastic for similarity-based models, and stochastic for us). That is why one complete section is devoted to this comparison. Hereby, the similarity-based prognostic method is used as a benchmark to investigate the relevance of our stochastic model. This method differs in the way the degradation indicator is modelled and the way the uncertainties are taken into account. More precisely, the aim of this section is to show that the performance of our method are not only due to the relevance of the indicator Di,½1:ri  , but also to the use of a stochastic process for the RUL estimation. With the similarity-based models, the main idea is to construct many off-line possible sub-models for degradation indicator and then to choose on-line which sub-models is the good one, given the monitoring information. For example in [34,41,42,19], the authors constructed the sub-models to estimated the RUL component and the estimated RUL is given by a similarity-weighted sum. Based on the advantages of similaritybased method, the authors [15,16,33] used the degradation indicator in [34] for RUL estimation by a similarity-based interpolation (SBI). Recently, in [3], the RUL of system is estimated by the ensemble of bootstrapped models and similarity-based method. Hence the similarity-based method for data-driven prognostic model is popular and many variations exist. We select and present in our paper the one that gives the best result for the PHM 2008 challenge data. 4.1. Stochastic prediction model 4.1.1. Estimation of Wiener process parameters—training data set In this study, the Wiener process is used for the degradation modelling corresponding to PHM data, refer to [27] for more details. The distances between the failure and the operational states classified in the principal component analysis (PCA) satisfies the properties of a Wiener process. Here, it is supposed that Di,½1:ri  ,r i Z 1 for unit i is one realisation over ½1,ri  of a unique Wiener process (non-stationary Brownian motion with drift) with three parameters ðZ, s, aÞ, that are the same for all the units. Hence the variance of the paths of Di,½1:ri  observed in Fig. 9 for i ¼ 1, . . . ,218 is modelled by the variance of the stochastic process. By definition this process noted D,j at time j satisfies the following conditions: 1. fD,j ,j Z1g has independent increments. 2. For every 0 os o j, D,j þ s D,j is normally distributed with mean Zððj þ sÞa ja Þ and variance s2 ððj þ sÞjÞ ¼ s2 s.

171

estimated as the optimal variables of the following loglikelihood function: l ¼ logðf Z, s, a ðDd1 , . . . , DdN ÞÞ ¼

¼

N X i¼1

0 log@

N X

logðg Z, s, a ðDdi,1 , . . . Ddi,ri ÞÞ

i¼1 ri Y j¼1

ðDdi,j Zððj þ 1Þa ja ÞÞ2 1 pffiffiffiffiffiffiffiffiffiffiffiffi exp 2s2 2ps2

!!1 A ð11Þ

where N is the number of units, f Z, s, a is the probability density function of all increments and g Z, s, a is the probability density function of all increments corresponding to each unit. The estimated parameters are obtained with the confidence interval of 95% as follows:

Z^ ¼ 1:2525e05, s^ ¼ 0:5086, a^ ¼ 2:0734

ð12Þ

These parameters will be used for the simulation of a nonstationary Wiener process from the last cycle of the testing data to the failure time in the order to find the estimated RUL for each unit. The next section will introduce the RUL estimation method by Monte Carlo simulations. 4.1.2. RUL estimation by simulation of a Wiener process—testing data set The testing data set is used to calculate the remaining useful life of each unit. The last cycle of the testing data set is the time upon which we calculate the RULs. Before setting the criterion of RUL estimation, we propose a method for the RUL simulation. Let kðr 0 Þ be r i0 the last cycle of the testing data set, Di0 ,ri0 the degradation i indicator at time r i0 of unit i0 , L a fixed degradation level and TL be the crossing time of the level L. If it is supposed that L is the failure threshold then L¼0. The Wiener process is nonmonotonous and therefore the process can cross several times the same threshold. This property can be noticed in Fig. 12 where a simulated path of the non-stationary Wiener process with the estimated parameters ðZ^ , s^ , a^ Þ obtained in the previous section shows many passage times from the degradation level L. In this framework considering the failure time as the first passage time or as the last passage time to L might overestimate or underestimate the RUL. The choice of the criterion for the RUL simulation is crucial and has a big influence on the accuracy of RUL estimation. First we consider the failure time as the first passage time. In [31], the authors introduce a mathematical theory of first hitting time with a diffusion non-linear process. Our proposed Wiener process with the nonlinear part : Zt a is a special case of the proposed model in [31]. Applying to our model with the negative drift Z and the failure threshold L, we obtain the density of first hitting time at the inspection time r i0 with the available degradation measurement Xn as follows: " # 1 ðZdðhÞLn Þ2 f RULðri0 Þ ðhÞ ¼ C pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðZ  ðdðhÞahðh þ r i0 Þa1 ÞLn Þexp  3 2hs2 2ph s2

ð13Þ

In this paper only Monte Carlo simulations are used, so the analytical formulation and the modeling of the first increment are postponed to future works. Let Di,j be the degradation indicator of component i at an arbitrary time j. The increments of degradation Ddi,j ¼ Di,j þ 1 Di,j , i ¼ 1, . . . ,N, for each unit follow a Gaussian distribution a N ðZððj þ 1Þa j Þ, s2 Þ. Let us consider the observations, Ddi ¼ ðDdi,1 , . . . Ddi,ri Þ, for i ¼ 1, . . . ,N where ri is the last sensors measurements before the failure of unit i. By using the maximum likelihood method [21], the Wiener parameters ðZ, s, aÞ are

a

r a0

and where C is the normalisation constant, dðhÞ ¼ ðh þr i Þ i Ln ¼ LX n . By testing this analytical function of first hitting time and Monte Carlo simulations, we obtain the same results. For instance in Fig. 11 we have two cases of Wiener parameters ðZ ¼ 1, a ¼ 1:5, s ¼ 0:2Þ and ðZ ¼ 1, a ¼ 1:5, s ¼ 2Þ where the FHT density in [31] are compared with distributions obtained by Monte Carlo simulations. By applying the proposed FHT density to our degradation indicator, we obtain the estimated RUL for each testing unit. However the obtained results S¼ 17 199 and MSE¼1341 are not 0

172

K. Le Son et al. / Reliability Engineering and System Safety 112 (2013) 165–175

0.5

3.5 Wiener : (η = 1, α = 1.5, σ = 0.2)

Wiener : (η = 1, α = 1.5, σ = 2)

0.45

3

0.4

2.5

0.35 0.3

2

0.25 1.5

0.2 0.15

1

0.1

0.5 0 1.4

0.05 1.6

1.8

2

2.2

2.4

2.6

0 0

2

4

6

8

10

12

Fig. 11. Comparison of the Monte Carlo simulation with the proposed FHT density in [31].

Table 2 The estimated remaining useful life time and its standard deviation.

Degradation Indicator

D−,t

Unit i0

RULi0estimated

si

1 2 3 4

97 78 84 104

0.76 0.53 0.50 0.14

0

min (t|D−,r+t < L)

L

to our model, we obtain Z 1 F T i0 max ðhÞ ¼ C½ZðhÞ=h2  þ

RUL = (TLmin+TLmax) / 2 − r

L

max(t|D −,r+t > L) TLmin

r

TLmax

Fig. 12. The proposed criterion of RUL calculation by Wiener simulation.

satisfying comparing to the results obtained by the winners of the PHM conference challenge. To improve our results, in this paper it is considered that the simulated RUL at time r i0 is the duration between the mean of the first passage time and the last passage time to degradation level L and r i0 , depicted in Fig. 12, and defined for unit i0 and one simulated path of Wiener process as follows: RULi0simulated ¼

ðT i0L minþ T i0L maxÞ r i0 2

ð14Þ

rLg and T i0L max ¼ inffj 4r i0 8s 4 0, where T i0L min ¼ inffj 4r i0 ,DkðjÞ i0 ,j kðj þ sÞ Di0 ,j þ s r Lg are respectively the first passage time and the last passage time to the degradation level L. Therefore, PðRULi0simulated o hÞ ¼ PðT i0L minþ T i0L maxÞ o 2ðhþ r i0 ÞÞ

ð15Þ

which can be calculated by the convolution of the distribution functions of Tmin and Tmax. The distribution function of FHT can be calculated from (13) and the distribution function of last hitting time (LHT) can be defined as follows: ! PðT i0L maxo MÞ ¼ P

sup fX h 4 Lg o M

ð16Þ

h 0

1

Ff½ZðhÞZðtÞ=ðhtÞ2 gF T i0 max ðtÞ dt L

ð17Þ pffiffiffiffiffiffi R x where FðxÞ ¼ ð1= 2pÞ 1 expðu2 =2Þ du and CðxÞ ¼ 1FðxÞ. This mathematical formula cannot be calculated directly and requires a heavy numerical computation method hence to bypass these complexity problems and avoid the time consuming methods the Monte Carlo simulations are used in our paper to evaluate the RUL distribution. By n ¼10 000 simulations of a non-stationary Wiener process from the last time of the testing data set unit’s, the estimated RUL, RULi0estimated , can be obtained as the empirical mean of the path simulated RULs. RULi0estimated ¼

n 1X RULi0,s simulated ns¼1

ð18Þ

0 where RULi0,s simulated is the sth simulated RUL of unit i given by Eq. (14) and n is the number of simulations. For example, the estimated RUL of four units are given in Table 2 with their 0 empirical standard deviation s i . By using the criterions of PHM Conference Challenge, the final obtained score S of the Wiener process plus Monte Carlo simulations method is 5575, the root mean squared error RSE ¼423 and the Mean Squared Error MSE¼ 819. Hence the performance indicator obtained is better than the one proposed by the winners of the 2008 PHM Data Challenge Competition. The study hereafter aims at showing that this good performance of our prognostic method is not only due to the relevance of the degradation indicator Di,½1:ri  , but also to the use of a stochastic process for the RUL estimation.

r i0 r h

The LHT distribution function for a standard Brownian motion is considered in [28]. However, in the case of Wiener process with nonlinear drift, the calculation of LHT is very complicated and time consuming. In [25], the authors evaluated the crossing time distribution of standard Wiener process to certain curves on various intervals to give the LHT distribution function. Applying

4.2. Similarity-based prognostic method For our proposed comparative study, a ‘‘mixture’’ of both methods is considered, taking our indicator Di,½1:ri  ,r i Z 1 to model the degradation and taking the similarity-based prognostic method proposed by Wang [34] for the RUL prediction. Hence, if

K. Le Son et al. / Reliability Engineering and System Safety 112 (2013) 165–175

the use of a stochastic process is really relevant, our complete method must get a better score than this one based on a ‘‘mixture’’.

173

Obsevations

1.6

Fitted curve

4.2.1. Construction of a collection of 218 models—training data set The main idea of the prediction method proposed by Wang [34] is to construct a model Mi for each unit of the training data set and then to calculate a RUL using the collection of these 218 models in the testing data set. To construct Mi ,i ¼ 1, . . . ,218, the degradation indicators are presented in the inverse time for all units. This operation is a way to concentrate the uncertainty on the state of the units at the starting time and not at the failure time. We consider 0 as the failure time, thus all the other times are negative. In Fig. 13, the degradation indicators in the inverse times are depicted. Afterwards, the degradation indicator trend of each component is fitted to an exponential (non-linear) regression model. For this model, the degradation indicator trend of unit i is fitted to Mi defined as follows: M i : f i ðtÞ ¼ ai ðebi t þ ci eci Þ, r i þ 1r t r0

Degradation indicator

1.4 1.2 1 0.8 0.6 0.4 0.2 0 −200

−150

−100

−50

0

Inverse cycle Fig. 14. The degradation indicator filtering by the regression model.

ð19Þ Move the block of degradation indicator along time, find the most probable position with regard to the curve of model Mi

Degradation indicator

where ri is the last sensors measurements before the failure of unit i in a non-inverse time scale and t ¼0 corresponds to the failure time in the inverse time scale. The fitted model of each unit i gives us a parameters set ðai ,bi ,ci Þ and a variance P s2i ¼ ð1=ri Þ 0t ¼ ri þ 1 ðf i ðtÞD i,t Þ2 , where Di,t is the degradation indicator of unit i at time t in the inverse time scale. For instance, in Fig. 14 the non-linear regression fitting for unit 1 is depicted, the estimated parameters are (ai ,bi ,ci )¼(  0.60, 0.02, 0.42). Finally, the estimated parameters set of the exponential regression model is used as a library in order to estimate the RULs of testing data set.

Model Mi RU Li

4.2.2. RUL estimation with the collection of 218 models—testing data set For the RUL estimation we calculate for a unit i0 of the testing data set, the degradation indicator Di0 ,½1,ri0  where r i0 is here the last

Degradation indicator of a test unit i

Cycle

date of sensors measurements before failure. Then we represent Di0 ,½1,ri0  by its inverse time scale transformation noted D i0 ,½1,ri0  . For one unit i0 of the testing data set, the distance between any model Mi of the training data set and the degradation indicator built kð1Þ

kð2Þ

kðr 0 Þ

from the testing data set D i0 ,½1,ri0  ¼ ½D i0 ,1 ,D i0 ,2 , . . . ,D i0 ,ri0  is defined

0

Fig. 15. RUL estimation criterion.

as the following Euclidean distance: dðt,D i0 ,½1,ri0  ,M i Þ ¼

i

r i0 X

kðjÞ

ðD i0 ,j f i ðtr i0 þ jÞÞ2 =s2i ,

0 r t r r i r i0

j¼1

ð20Þ

3 2 i

where s is the prediction variance given by the model Mi. In this case, dðt,D i0 ,½1,ri0  ,M i Þ is a function of t, the number of cycles that

2.5

Degradation indicator

the sequence D i0 ,½1,ri0  is shifted away from cycle zero of the model Mi. The distance measure tells the similarity level that a test unit behaves as model Mi at its history t. Closer D i0 ,½1,ri0  is to Mi higher

2

is the similarity. By considering this distance the RUL of unit i0 given by model Mi can be estimated as follows:

1.5

RULi0 ,i ¼ arg min dðt,D i0 ,½1,ri0  ,M i Þ

ð21Þ

t

1

The collection of models M i ,i ¼ 1, . . . ,218 of the training data set give us a sub-set of RULi,i0 ,i ¼ 1, . . . ,218 corresponding to the unit

0.5

0 −400

−350

−300

−250

−200

−150

−100

−50

Inverse cycle Fig. 13. The degradation indicator in the inverse times.

0

i0 depicted in Fig. 15. The final RUL of each unit i0 is estimated by the weighted sum of the obtained RULs. X X RULi0 ¼ wi  RULi,i0 , where wi ¼ 1 ð22Þ i

i

In the framework of PHM conference (see at the score formula (1)), we choose the weights wi as proposed by Wang [34], the

174

K. Le Son et al. / Reliability Engineering and System Safety 112 (2013) 165–175

following weighting method is used, which emphasises on the upper and lower boundaries only (minimum and maximum of RULi,i0 and considering the exponential penalty to prediction errors). RUL ¼ ð13=23Þ  minðRULi Þ þ ð10=23Þ  maxðRULi Þ i

i

ð23Þ

By using the criterion above with all units, we obtain the RUL of each unit. Finally, the score is S¼6690. Hence, the performance of our probabilistic model (S ¼5520) is also due to the use of a stochastic process for the RUL estimation and not only on the proposed degradation indicator. For the similarity based prognostic method, we also calculated the value of MSE¼809 and RSE ¼420 which are slightly better than the ones corresponding to the Wiener based method. 4.3. Discussion The scores of our method based on Wiener process are S¼ 5520, RSE ¼423 and MSE¼819. These results show the relevance of a probabilistic approach using stochastic process by comparison with a probabilistic approach with lifetime model (S¼9870) and by comparison with a similarity based approach using the same degradation indicator (S ¼6690, RSE ¼420, MSE¼809). The MSE and RSE of the Wiener based method are not better than the ones corresponding to the similarity based prognostic method but we obtain a better score S with the Wiener based method. The reason is that in the criterion of score, we focus on the early prediction, it means that our model obtains more early predictions than the similarity based prognostic method. These last results seem to justify that for the proposed degradation indicator, the use of our probabilistic method based on the stochastic process can be more appropriate in the framework of maintenance optimisation. We can also compare our results with best obtained results of the PHM 2008 conference challenge as follows:

 Non-probabilistic model with the similarity-based prognostic 

method of Wang in the 2008 PHM conference [34]: the score is S¼5636. Non-probabilistic model based on the neural networks of [26] and [14] show respectively MSE¼984 and an error of 519.8.

These results and our obtained score justify the accuracy of our probabilistic method. A probabilistic approach seems to be promising since the result are comparable to the best ones, and our degradation indicator makes sense. What is more, a Wiener process seems to be relevant compared to neural networks and similarity-based prognostic method. This first study encourage us to improve our model by using other stochastic processes for a better RUL estimation (Gamma process with noise for example) and to model more precisely the evolution of the operational modes (with covariates for example).

5. Conclusion This paper shows a probabilistic prognostic method on the 2008 PHM conference challenge data set by using a Principal Component Analysis in order to find a better degradation indicator and model it by a Wiener process. The obtained result is better than the one given by a simple lifetime modelling based on only the historical failure data [10] and is as good as the winning RUL estimations methods of the PHM Data Challenge. Hence the advantage of the use of stochastic processes for the prognostic is pointed out.

In further works in order to improve the obtained results on the case study, a more sophisticated stochastic model (gamma process with a Gaussian noise) can be considered. This model requires more complex probabilistic tools to deal with uncertainties but in this case, the RUL estimation seems easier and less time consuming. Another path to follow is to take into account the operational mode changes with more complex models. Once these ideas explored, it could be interesting to use our RUL estimation method in the framework of a prognosis-based maintenance policy and to study the properties of the policy under consideration. All these points constitute important perspectives and our future works will be devoted to their development.

Acknowledgments This work is part of the PhD research work of Khanh Le Son financially supported by Conseil Re´gional de ChampagneArdenne, France, and by Groupement d’Inte´rˆet Scientifique ‘‘Surveillance, safety and security of the big systems’’. References [1] Abu-Hanna A, Lucas P. Prognostic models in medicine. Methods of Information in Medicine-Methodik der Information in der Medizin 2001;40(1):1–5. [2] Banjevic D. Remaining useful life in theory and practice. Metrika 2009;69(2): 337–49. [3] Baraldi P, Mangili F, Zio E. Ensemble of bootstrapped models for the prediction of the remaining useful life of a creeping turbine blade. In: 2012 IEEE international conference on prognostics and health management, June 18–21, Denver, Colorado; 2010. [4] Bond L, Ramuhalli P, Tawfik M, Lybeck N. Prognostics and life beyond 60 years for nuclear power plants. In: 2011 IEEE international conference on prognostics and health management; 2011. p. 1–7. [5] Bond L, Taylor T, Doctor S, Hull A, Malik S. Proactive management of materials degradation for nuclear power plant systems. In: IEEE international conference on prognostics and health management, 2008. PHM 2008; 2008. p. 1–9. [6] Brown D, Kalgren P, Byington C, Roemer M. Electronic prognostics—a case study using global positioning system (GPS). Microelectronics Reliability 2007;47(12):1874–81. [7] Byington C. Prognostic enhancements to diagnostic systems (PEDS) applied to shipboard power generation systems. Technical report, DTIC Document; 2004. [8] Carbone R, Armstrong J. Note. Evaluation of extrapolative forecasting methods: results of a survey of academicians and practitioners. Journal of Forecasting 1982;1(2):215–7. [9] Chhikara R, Folks J. The inverse Gaussian distribution as a lifetime model. Technometrics 1977;19(4):461–8. [10] Coble J, Hines J. Prognostic algorithm categorization with PHM challenge application. In: International conference on prognostics and health management, 2008. PHM 2008; 2008. p. 1–11. [11] Cox D, Oakes D. Analysis of survival data. Chapman & Hall; 1984. [12] Doksum K, Ho´yland A. Models for variable-stress accelerated life testing experiments based on Wiener processes and the inverse Gaussian distribution. Technometrics 1992:74–82. [13] He D, Bechhoefer E. Development and validation of bearing diagnostic and prognostic tools using hums condition indicators. In: 2008 IEEE aerospace conference. IEEE; 2008. p. 1–8. [14] Heimes F. Recurrent neural networks for remaining useful life estimation. In: International conference on prognostics and health management, 2008. PHM 2008; 2008. p. 1–6. [15] Hu C, Youn B, Wang P. Ensemble of data-driven prognostic algorithms for robust prediction of remaining useful life. In: 2011 IEEE conference on prognostics and health management (PHM). IEEE; 2011. p. 1–10. [16] Hu C, Youn BD, Wang P, TaekYoon J. Ensemble of data-driven prognostic algorithms for robust prediction of remaining useful life. Reliability Engineering & System Safety 2012;11:22. [17] Hyndman R, Koehler A. Another look at measures of forecast accuracy. International Journal of Forecasting 2006;22(4):679–88. [18] Jolliffe I. Principal component analysis. Springer Verlag; 2002. [19] Liu J, Djurdjanovic D, Ni J, Casoetto N, Lee J. Similarity based method for manufacturing process performance prediction and diagnosis. Computers in Industry 2007;58(6):558–66. [20] Mazhar M, Kara S, Kaebernick H. Remaining life estimation of used components in consumer products: life cycle data analysis by Weibull and artificial neural networks. Journal of Operations Management 2007;25(6):1184–93. [21] Moler C. Numerical computing with MATLAB. Society for Industrial Mathematics; 2004.

K. Le Son et al. / Reliability Engineering and System Safety 112 (2013) 165–175

[22] Muller A, Suhner MC, Iung B. Formalisation of a new prognosis model for supporting proactive maintenance implementation on industrial system. Reliability Engineering & System Safety 2008;93(2):234–53. [23] Padgett W, Tomlinson M. Inference from accelerated degradation and failure data based on Gaussian process models. Lifetime Data Analysis 2004;10(2): 191–206. [24] Palazzolo J, Scheunemann L, Hartin J. Leakage fault detection method for axial-piston variable displacement pumps. In: 2008 IEEE aerospace conference. IEEE; 2008. p. 1–8. [25] Park C, Schuurmann F. Evaluations of barrier-crossing probabilities of wiener paths. Journal of Applied Probability 1976;11:267–75. [26] Peel L. Data driven prognostics using a Kalman filter ensemble of neural network models. In: International conference on prognostics and health management, 2008. PHM 2008; 2008. p. 1–6. [27] Revuz D, Yor M. Continuous Martingales and Brownian motion. New York: Springer-Verlag; 1999. [28] Salminen P. On the first hitting time and the last exit time for a Brownian motion to/from a moving boundary. Advances in Applied Probability 1988;11:411–26. [29] 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. [30] Si X, Wang W, Hu C, Zhou D. Remaining useful life estimation—a review on the statistical data driven approaches. European Journal of Operational Research 2011;213(1):1–14. [31] Si X, Wang W, Hu C, Zhou D, Pecht M. Remaining useful life estimation based on a nonlinear diffusion degradation process. IEEE Transactions on Reliability 2012;61(1):50–67.

175

[32] Tseng S, Peng C. Stochastic diffusion modeling of degradation data. Journal of Data Science 2007;5(3):315–33. [33] Wang P, Youn B, Hu C. A generic probabilistic framework for structural health prognostics and uncertainty management. Mechanical Systems and Signal Processing 2011;28:622–37. [34] Wang T, Yu J, Siegel D, Lee J. A similarity-based prognostics approach for remaining useful life estimation of engineered systems. In: International conference on prognostics and health management, 2008. PHM 2008. p. 1–6. [35] Wang X, Rabiei M, Hurtado J, Modarres M, Hoffman P. A probabilistic-based airframe integrity management model. Reliability Engineering & System Safety 2009;94(5):932–41. [36] Wang W. A two-stage prognosis model in condition based maintenance. European Journal of Operational Research 2007;182(3):1177–87. [37] Wang W, Carr M, Xu W, Kobbacy K. A model for residual life prediction based on Brownian motion with an adaptive drift. Microelectronics Reliability 2011;51(2):285–93. [38] Wang X. Wiener processes with random effects for degradation data. Journal of Multivariate Analysis 2010;101(2):340–51. [39] Whitmore G. Estimating degradation by a Wiener diffusion process subject to measurement error. Lifetime Data Analysis 1995;1(3):307–19. [40] Whitmore G, Schenkelberg F. Modelling accelerated degradation data using Wiener diffusion with a time scale transformation. Lifetime Data Analysis 1997;3(1):27–45. [41] Zio E, DiMaio F. A data-driven fuzzy approach for predicting the remaining useful life in dynamic failure scenarios of a nuclear system. Reliability Engineering & System Safety 2010;95(1):49–57. [42] Zio E, Peloni G. Particle filtering prognostic estimation of the remaining useful life of nonlinear components. Reliability Engineering & System Safety 2011;96(3):403–9.