untitled A publication of CCHHEEMMIICCAALL EENNGGIINNEEEERRIINNGG TTRRAANNSSAACCTTIIOONNSS VOL. 33, 2013 The Italian Association of Chemical Engineering Online at: www.aidic.it/cet Guest Editors: Enrico Zio, Piero Baraldi Copyright © 2013, AIDIC Servizi S.r.l., ISBN 978-88-95608-24-2; ISSN 1974-9791 Residual Useful Life Estimation Based on a Time-Dependent Ornstein-Uhlenbeck Process Yingjun Deng, Anne Barros, Antoine Grall Institut Charles Delaunay, Université de Technologie de Troyes, CNRS UMR STMR, 12 rue Marie Curie, CS 42060, 10004 Troyes Cedex , France antoine.grall@utt.fr Residual useful life plays a core role in the field of the prognostics and health management of systems. For deteriorating systems or devices it can be estimated by establishing a good-fitting stochastic process model and solving corresponding first passage problem with the given failure level. In this paper, providing that the mean and variance function of the degradation process are given or obtained directly from observed data by statistical/empirical techniques, we propose a time-dependent Ornstein-Uhlenbeck process whose mean and variance are equal to the given mean and variance respectively. Based on this diffusion process, when a failure level is considered, the first passage failure distribution can be solved by a Volterra integral equation of second class. Furthermore, with the help of numerical methods, the residual useful life estimation is considered for a real case. 1. Introduction Estimating the residual useful life (RUL) from degradation records is a core issue in many fields, e.g. nuclear plant, aeronautical industry, mechanical engineering etc. However in a complicated environment with unknown mechanisms, degradation records may appear largely fluctuating. Therefore it is hard to use usual models to explain everything, e.g. Paris-Erdogan model or Gamma process, and more and more people add fluctuating uncertainties into the degradation models (K. Sobczyk, 2004). In general, people consider the fluctuation in the degradation records from 2 ways. The first is to consider the uncertainties as random errors e.g. measurement errors to the original degradation mechanism, which should be independent with each other at different times. This is supposed by plenty of physical or engineering work (Robert E. Melchers, 1999). The other supposes that the internal uncertainty cause the fluctuation. However in this case, we should notice that because of the complicated internal mechanisms, the outside appearance of the internal uncertainties can be rather unpredictable. And also as the complexity of such models is beyond the previous consideration, it is hard to perform later work such as estimation of RUL and mean time to failure (MTTF). However this is an active field attracted by many authors (Si et al., 2012). Compared to existing stochastic degradation models, the model stated in this paper is different from 3 aspects. Firstly, the influence of the uncertainty at the initial state is considered and connected to later performance. Secondly, the variance of the process is almost a constant, which can provide a more appropriate way to describe the system development. Thirdly, this model has good properties for numerical solution: the first passage failure can be estimated based on efficient and stable numerical methods without the usage of Monte-Carlo simulation. 2. Stochastic Modeling of Degradation Processes 2.1 System Descriptions and Hypotheses In this paper, we concentrate on deteriorating systems subject to gradually degradation processes such as fatigue crack, corrosion, erosion etc. And we suppose there exist mainly 2 kinds of uncertainties in deterioration model: measurement error and deterioration mechanism. Both of them are considered here, however the measurement error is treated as the internal uncertainty but not as an additive noise. Two 325 DOI: 10.3303/CET1333055 Please cite this article as: Deng Y., Barros A., Grall A., 2013, Residual useful life estimation based on a time-dependent ornstein- uhlenbeck process, Chemical Engineering Transactions, 33, 325-330 DOI: 10.3303/CET1333055 reasons motivate this simplification. In one hand in stochastic deteriorating systems, uncertainties at different times generally are dependent with each other. In the other hand, the internal uncertainty is a common phenomenon in degradation cases. For instance in crack and related fatigue phenomenon, a very important phenomenon named crack closure can cause fluctuation in the crack growth records even with perfect measurement (K. Sobczyk, 2004). Also a recent work on the impact of lightning in corrosion failure can also support this system hypothesis (Necci et al., 2012). To summarize the above and for mathematical simplification, the system in consideration is supposed to have the following properties: 1) The system deteriorates gradually without the risk of sudden shocks e.g. earthquakes, tsunami etc. 2) The initial system state cannot be observed explicitly, or has some uncertainties. 3) The degradation records have obvious fluctuations caused by internal system mechanisms and errors. 2.2 Time-dependent Ornstein-Uhlenbeck Process Following the hypotheses of the system in previous section, noting tx as the degradation level at time t , here we suppose that tx satisfies the following time-dependent Ornstein-Uhlenbeck process: ,))()('( ttt dBdttrmtmrxdx σ+++−= (1) where σ,r are non-negative real constants, ),0[)( 1 +∞∈Ctm , tB is a standard Brownian motion. The initial value 0x is a random variable with cumulative distribution function )(0 xF and probability density function(pdf) )(0 xf , moreover it has variance 0v and mean 0m . The solution of the equation (1) can be given by: .))(exp()exp()()( 000 −+−−+= t ut dButrrtmxtmx σ (2) And therefore the mean and covariance of tx are given when st ≥ : )).2exp(2))((exp( 2 ),(),()( 2 2 0 2 srvstr r xxCovtmxE stt +−+−== σ σσ (3) Remarks: When rv 2/2 0 σ= , then from (3) the variance of tx is a constant. We should also notice that this process has the mean-reverting and gaussian properties. Therefore it can be a substitute of wildly used noise- added model, as it reserves consideration of fixed-variance uncertainty while it also includes the interaction of uncertainties at different times as the covariance is not zero. Moreover even 0x is arbitrary, the asymptotic distribution of tx is gaussian with mean )(tm and variance r2/2σ . Moreover, for the process tx , the transition probability density ),|,(:),( sytxptxp = satisfies the Fokker- Planck equation: ,),( 2 )])(,()('[),( 2 22 x txpmxtxrptm xt txp t ∂ ∂+−− ∂ ∂−= ∂ ∂ σ (4) where ),-()0,(,0),( yxxptp δ==±∞ )(⋅δ is the Dirac mesure concentrated at y. And the analytic answer to (4) is given by: ), ),(4 )),((exp( ),(4 ),|,( 2),(),( st ystxe st esytxp stst γ β πγ αα −+−= (5) where γβα ,, are given by 326 ).1))(2(exp( 2 ),( )),(exp()()(),( ),(),( 2 −−= −−= −= str r st strtmsmst strst σγ β α (6) Remarks: when the initial value is a random variable, the transition pdf from 0x is an average of (5) with the form .)()0,|,()0,|,( 00 +∞ ∞− = dyyfytxpxtxp This can be calculated numerically especially when 0x has a truncated distribution. 2.3 Parameter Estimation In this section we will state how to estimate the parameters in (1) with given degradation data-set by maximum likelihood method. Supposing the degradation records are recorded on n independent equipments, the records are noted by }0,1|){( , iijij mjnitx ≤≤≤≤ , where )( , ijij tx stands for the record at the j th time in the i th trajectory. Moreover there exist totally n trajectories with im records in the i th trajectory respectively, and we set .1,00 niti ≤≤= If we set )(tm as a 3-parameter real function, e.g. 0)1)1(()( mtatm b +−+= , then (1) becomes a 5- parameter model with parameter ),,,,( 0 σθ rmba= . When records for the i th trajectory )( , ijij tx are given, the conditional PDF of tx is based on the observation of whose PDF appears as a Dirac delta distribution, so from (5), we know that: ), ),(4 )),(( exp( ),(4 ),|,( )1( 2 )1( ),( )1( )1( ),( )1()1( )1()1( ijji ijijji tt ji ijji tt ijijjiji tt xttex tt etxtxp ijjiijji + ++ + ++ −+ −= ++ γ β πγ αα (7) where γβα ,, are given by (6). Then we can construct the log-likelihood function for the i th trajectory as follows: .)),,|,(log())|(log()(log 1 0 )1()1(00 − = +++= im j ijijjijiii txtxpxfL θθθ (8) Moreover, the log-likelihood function for the whole data-set is given: .)),,|,(log())|(log()(log 1 1 0 )1()1( 1 00 = − = ++ = += n i m j ijijjiji n i i i txtxpxfL θθθ (9) From (9), we can maximize )(log θL to get the optimal estimated parameters. 2.4 Residual Useful Life Calculation Various ways exist to describe a system failure. In our research, we concentrate on deterioration phenomenons, e.g. corrosion, erosion, crack et al., in which deterioration cumulates until it reaches the material strength. So we will consider the first passage failure throughout: given an alert level of deterioration, a system failure occurs once accumulating deterioration exceeds the alert level. The above question is essentially the classical first passage time (FPT) problem in a mathematical view. That is to say, if tX is a given stochastic process, Ω is a given safe area, we want to investigate the first time when tX leaves the area Ω based on the present observation ),( sy , i.e. }.|:{inf yXXt stRt =Ω∉= +∈ τ The residual useful life therefore can be expressed by s−τ . Also we can define the reliability of the system at time t simply by )( tP >τ . Here we just consider the one-dimension problem, given the boundary as ),(C)( 1 +∈ RtL where the safe area is supposed to be )}({ tLx ≤=Ω . Therefore to derive the distribution function ),|),(( 00 txttLg for the first passage time τ , based on the initial value ),( 00 tx , we refer to the numerical solution to a Volterra integral equation of second class (E. Pirozzi et al., 2001) as follows: ,)),(|),((),|),((2),|),((2),|),(( 0 000000 dsssLttLKtxssLgtxttLKtxttLg t t +−= (10) 327 where ).,|),((] ),(2 ),()( 22 )()(')()('[),|),(( ),(),(2 syttLp st ystetLetrmtmtrLtLsyttLK stst γ βσ αα −+−−−+= (11) with the same notations as in equations (5) and (6). Then equation (10) is with a non-singular kernel, and therefore can be numerically solved by several methods, such as compound trapezoid formula with the following form: )),(|),((2),|),((2 0000 1 1 0000 tjtjLtktkLKgtxtktkLKg k j jk +Δ+Δ+Δ+ΔΔ++Δ+Δ−= − = , (12) where .,),,|),((,0),,|),(( 0 0000000 Nk k tttxttLgggtxjtjtLgg kj ∈−=Δ==Δ+Δ+= (13) Remarks: when we consider the random initial value 0x , ),|),(( 00 txttLK has an integral form which is hard to be given analytically, therefore we just consider the pdf of 0x is truncated in the interval ],[ ba , and therefore we should calculate it with numerical integral scheme, e.g. compound trapezoid formula with the following form: )].(),|),(( 2 1)(),|),(( 2 1 )(),|),(([)( 00000000 0000 1 0 bftbtkhtkhLKaftatkhtkhLK jhafjhtjhatkhtkhLKhtK q j q ++++++ +++++= − = (14) 3. Study of a Real Degradation Case According to previous analysis, in the following we consider a real case study where the degradation data-set is chosen as 415 records of degradation in 159 independent equipments. And we do the following tests in a personal computer with R software (version 2.15.1). Figure 1 and Figure 2 express the degradation records in an intuitive way, where Figure 1 expresses the records of degradation, Figure 2 expresses the recorded trajectories on different equipments. Single points in Figure 2 mean that only one record is recorded on that equipment. The measurement scale of the time adopted here is 1000 hours per unit, and the degradation level here is an integrated index without real measurement scale. Figure 1: Degradation records. Figure 2: Recorded trajectories on each component. 328 Figure 3: Mean and standard variance. Figure 4: Several trajectories of the fitted process. Figure 5: Probability density of the first-passage failure time: red (12), blue(25), black(50). Figure 6: Cumulative density of the first-passage failure time: red (12), blue(25), black(50). We here suppose that the records is subject to the stochastic process defined in Eq (1), and a general description of the chosen model can be seen from Table 1. As the initial information is missing, we here suppose that initial value 0x is uniformly distributed in ]12,0[ 0v , where r v 2 2 0 σ= . We also suppose 03)1)1(()( vtatm b +−+= , therefore (1) now is a 4-parameter model with parameters ),,,( rba σ . From the remarks in section 2.2, the constructed model is a stochastic process with constant variance 0v , mean )(tm , and covariance between tx and sx is tsstrv ≤−− )),(exp(0 . Table 1: General description of the chosen model Mean Variance Covariance Initial distribution Unknown parameters 03)1)1(( )( vta tm b +−+ = r v 2 2 0 σ= ))(exp(0 strv −− uniformly distributed in ]12,0[ 0v ),,,( rba σ 329 Using the default Nelder-Mead algorithm in R with the searching start )1,6,35.1,6.0( ==== rba σ and the initial interval is divided to 1000 parts, the parameters are optimized and the results can be seen from Table 2, where 0v =10.36483 . In Figure 3, the red curve is the expression for the fitted )(tm , and the upper and lower blue curves express 0)( vtm ± respectively. The upper and lower black lines are for 03)( vtm ± . Here if we treat tx as a pure gaussian process (the asymptotic distribution of tx is gaussian distributed), then a single point has a probability of 0.692 to appear in the area between 2 blue curves, has a probability of 0.999997 to appear in the area between the 2 black curves. From this probability view, Figure 3 fits the hypothesis on the gaussian property of tx . Figure 4 shows several simulated trajectories of the fitted process. Table 2: Fitted parameters Parameter a B σ r r2/v 2 0 σ= Value 0.6044493 1.3966703 2.2810409 0.2505040 10.38536 Noticing the technical procedures for estimating MTTF and RUL are generally identical, and also the difficulty to estimate tens of elements’ RUL, we here just consider to estimate the MTTF instead. Moreover here instead of estimating one equipment’s MTTF, we consider the average MTTF for all equipments. We select 3 constant boundaries )(tL : 12, 15 and 20, and consider the corresponding first passage time as the element’s failure time. In Figure 5, the probability density and cumulative density function for the first passage failure are compared, where the red, blue, black curves stand for the boundaries 12, 15, 20 respectively. The average MTTF for all equipments is estimatedand the results can be seen from Table 3. Table 3: MTTF estimation based on first-passage failure time. Boundary L=12 L=25 L=50 MTTF 3.79315 10.75301 20.58544 4. Conclusion In this paper, we presented a time-dependent Ornstein-Uhlenbeck Process for modeling degradation processes, with which the RUL can be estimated. This model considers not only the uncertainty in the observations, but also the interaction of the uncertainties at different times. Another contribution is to add the uncertainty in initial system state into the modeling work. Acknowledgement We greatly appreciate the useful discussion with Prof. Michel Roussignol in University of Paris-Est Marne-le- Valley, and also the partial financial support from China Scholarship Council. References Amos Necci, Giacomo Antonioni, Elisabetta Renni, Valerio Cozzani, Alberto Borghettib, Carlo Alberto Nuccib and Elisabeth Krausmannc, 2012, Equipment Failure Probability due to the Impact of Lightning, CHEMICAL ENGINEERING TRANSACTIONS, 26, 129-134. Pirrozi E., Di Nardo E., Nobile A. G. and Ricciardi L. M., 2001, A computational approach to first-passage-time problems for Gauss-Markov processes, Advances in Applied Probability, 33(4), 453-482. Sobczyk K., 2004, Fatigue crack growth in random residual stresss, International Journal of Fatigue, 26(11), 1179-1187. Robert E. Melchers, 1999, Corrosion uncertainty modeling for steel structures, Journal of Constructional Steel Research, 52(1), 3-19. Sheng-Tsaing Tseng and Chien-Yu Peng, 2007, Stochastic diffusion modeling of degradation data, Journal of Data Science, 5, 315-333. Xiao-Sheng Si, Wenbin Wang, Chang-Hua Hu, Donghua Zhou and Micheael G. Pecht, 2012, Remaining useful life estimation based on a nonlinear diffusion degradation process, IEEE Transactions on Reliability, 60(1), 50-61. 330