Remaining Useful Life Estimation Considering Prior Accelerated Degradation Data and Bayesian Inference for Multi-Stress Operating Conditions

Prediction of remaining useful life using the field monitored performance data provides a more realistic estimate of life and helps develop a better asset management plan. The field performance can be monitored (indirectly) by observing the degradation of the quality characteristics of a product. This paper considers the gamma process to model the degradation behavior of the product characteristics. An integrated Bayesian approach is proposed to estimate the remaining useful life that considers accelerated degradation data to model degradation behavior first. The proposed approach also considers interaction effects in a multi-stress scenario impacting the degradation process. To reduces the computational complexity, posterior distributions are estimated using the MCMC simulation technique. The proposed method has been demonstrated with an LED case example and results show the superiority of Bayesian-based RUL estimation. KeywordsRemaining useful life, Degradation data, Bayesian inference, Accelerated degradation test, Gamma process, MCMC simulation.


Introduction
Most of the recently designed engineering products and systems are highly reliable especially in the automotive, aerospace, electronics, and military industries. These highly reliable and high valued products are very critical and therefore, require extensive maintenance to ensure safer operations and reduce downtime. The remaining useful life (RUL) is a critical estimate that significantly influences the maintenance scheduling, spare parts inventory management, and overall prognostic operations (Pecht and Jaai, 2010). The RUL can be defined by the expected time remaining before a product fails (Si et al., 2011).
The majority of engineering products deteriorate over time due to external stress factors (Pan et al., 2016). This deterioration of physical systems also reflects on product performance. With the appropriate observation of the product's quality characteristic, the performance deterioration can be indirectly measured. The rapid development of sensor technologies and the emergence of the internet of things (IoT) have made the degradation data availability very convenient and feasible (Kwon et al., 2016). As a result of this, the prognostic and health management (PHM) of the high valued asset is getting more attention in the last few decades. Estimating the RUL of interested product or system is the main target of the prognostic analysis (Jardine et al., 2006;Pan et al., 2016).
In the existing literature, RUL estimations can be divided into three categories: physics-based, datadriven based, and hybrid types (Javed et al., 2015). However, the data-driven based stochastic process approach has become the most popular RUL estimation method due to well-developed mathematical models and the capability to capture the uncertainty (Son et al., 2013;Ye and Xie, 2015). For example, Wang et al. (2016) presented a RUL prediction model based on Wiener stochastic process for axial piston pump and Li et al. (2019) proposed RUL prediction based on Wiener stochastic process considering unit-to-unit variability. One of the drawbacks of the Wiener stochastic process is that assumption of degradation increment distribution can result in a negative value of degradation increment. The majority of the physical degradation processes such as wear, fatigue, corrosion are in fact nonnegative and monotonic degradation. This makes the gamma stochastic process a natural choice to model the degradation behavior for RUL prediction. In a gamma process, the shape and scale parameters describe the degradation rate and its variability.
Several works have been attempted to adopt the gamma process in maintenance decision making analysis (Dieulle et al., 2003;Crowder and Lawless, 2007;Pandey et al., 2009). Bagdonavicius and Nikulin (2001) presented a model with the gamma degradation process including the stress factor. Lawless and Crowder (2004) considered the random effect in addition to stress factors in their degradation model. Park and Padgett (2005) derived an approximation lifetime distribution for the gamma degradation model. Wei and Xu (2014) investigated the RUL estimation considering the measurement error of online monitored degradation data. Haowei et al. (2015) proposed a gammabased lifetime prediction model considering prior accelerated degradation information. Pang and Jia (2016) studied the gamma process-based RUL estimation to apply it for maintenance optimization. Recently, Ling et al. (2019) presented a RUL estimation work considering different phases of degradation behavior with Bayesian inference.
Besides the fact that the gamma process is well fitted to model degradation behavior, there has not been much work considering the gamma process for degradation modeling and RUL prediction especially considering the multi-stress environment. As technology becomes more complex every day, generally a product experiences more than one external stress factors (Zhu and Elsayed, 2013). Multi-stress factors also bring the issue of interaction between stresses that could lead a heavy nonlinear degradation behavior. Additionally, most of the existing literature considered only the shape parameter is stress-dependent and the scale parameter is independent of stress factors, which are questioned by recent works (Bayel and Mettas, 2010;Balakrishnan and Ling, 2014;Limon et al., 2020).
To address the research issues discussed earlier, this paper proposes a gamma process-based RUL model considering multi-stress factors with the interaction effect and stress-depended gamma parameters. During the product design and development stages, manufacturers test their products under harsher than the normal condition to expedite the failure process known as an accelerated life test (ALT) (Limon et al., 2017). If one is able to measure the degradation of performance characteristics during ALT without waiting for failure, it is called accelerated degradation testing (ADT). The degradation data collected during ADT could be very useful to predict the RUL of the product during normal operations. In the proposed approach, the Bayesian inference method has been used to predict the RUL. In the proposed Bayesian method, integrated degradation data collected during ADT, treated as prior information, and prognostic field degradation data used for the posterior update. The acceleration factor concept has been applied to convert the gamma model parameters from accelerated to normal operating conditions. The best-fitted distribution is considered for prior distribution. The Markov Chain Monte Carlo (MCMC) simulation technique has been implemented to evaluate the posterior parameters of Bayesian inference. Once the posterior parameters are updated for each new degradation observation, the RUL and its probability density function are estimated.
The rest of the paper is organized in the following manner. Section 2 provides the gamma degradation modeling and RUL prediction process using Bayesian inference. Section 3 provides the methodology of obtaining prior distributions for parameters. A case example with LED degradation data is demonstrated in section 4. Finally, the summary of the work presented in section 5.

Remaining Useful Life Model 2.1 Degradation Modeling
The degradation of quality characteristics behaves probabilistically over time. This uncertain behavior can be best suited to a model with well-known stochastic processes such as the Wiener process, gamma process, and inverse Gaussian process. Because of monotonic and nonnegative properties and meaningful physical interpretation, the gamma process is very appealing for modeling the product's degradation behavior. Mathematically, the gamma process has an independent and non-negative increment that follows a two-parameter gamma distribution. Considering a degradation process is described by the random variable Y at time t, the degradation increment by gamma process can be written as follow: Where, α is the shape parameter, β is the scale parameter, and Λ(t) is a monotonic increasing timescale function with Δ Λ( ) = Λ( + Δ ) − Λ( ). Now, the degradation path at any time t can be computed as: The degradation rate and its variability are defined by two parameters of the gamma process. In the condition that timescale function is linear (Λ( ) = ), the gamma process represents a stationary process otherwise it is a non-stationary process that is practical for most of the natural degradation processes. Considering historical empirical studies, the timescale function Λ(.) is modeled with the power-law model as Λ(t)=t c , where c represents the nonlinearity constant. While c=1 it represents a linear degradation process, c>1 represents a convex degradation process, and c<1 represents a concave degradation process.
Let, D be the constant degradation threshold value and thus the failure time ξ can be defined when cumulative degradation Y(t) reaches the threshold D as per the first passage of time concept of the stochastic process: For a gamma process, the failure time ξ is conveniently approximated by the Birnbaum-Saunders (BS) distribution (Park and Padgett, 2005). The probability density function (PDF) of failure time ξ with time transformation hence can be written as (Balakrishnan and Kundu, 2019): Where, ′ = 1 √ ′ ⁄ and ′ = ′ ⁄ are BS distribution parameters. ′ = and assume there is no initial degradation value. Now consider that the degradation path Y(tk) at time tk does not reach the threshold D. Then the remaining useful life RULtk at time tk is defined as (Zhang et al., 2018): Therefore, the PDF of the remaining useful life using Equation (4) can be expressed as: Here Dβ =(D-Y(t)) β. Using equation (6) the expected value of RUL can be written as:

Bayesian Inference Modeling
The RUL can be estimated once the estimates of model parameters α, β, and c are obtained. The abundance and credible real-time field data are generally not adequately available to obtain the parameter estimates. However, the manufacturer's prior accelerated degradation test (ADT) data can be used to estimate parameters. The Bayesian method can be utilized to get refined posterior estimates of α, β, and c in the presence of available individual product's field degradation data. Let, P (α, β, c) represents the joint prior distribution function where α, β, and c parameters are mutually independent. P(α), P(β), and P (c) represent the respective prior distribution of parameters α, β, and c, hence, by the mutual independent property, P (α, β, c) = P(α) P(β) P (c).
Assuming that field degradation data = [ 1 ( 1 ), 2 ( 2 ), … … … ( )] are collected from an individual product at corresponding time and (∆ | , , ) is the likelihood function of the data. Then the joint posterior distribution ( , , |∆ ) can be expressed by Bayesian theory as follows (Martz and Waller, 1982): The marginal distribution functions of |∆ , |∆ , and |∆ can be written as: The posterior expectation of |∆ , |∆ , and |∆ then can be estimated as: It is noted that the marginal distribution functions listed in Equations (9-11) are usually in unknown forms. Therefore, evaluation of the expected value of posterior distribution listed in Equation (12-14) with direct mathematical derivation is extremely complex or sometimes impossible. The Markov Chain Monte Carlo (MCMC) simulation is a good alternative to resolve this mathematical complexity (Gelman et al., 2004). Once we have the posterior expectation estimates of ̂,̂, , Equation (6) and (7) can be used to obtain the PDF and expected RUL, respectively. In the presence of new field degradation data, the expected posterior estimates will be updated to get a real-time RUL prediction.

Prior Estimates Using Accelerated Degradation Data
For posterior estimates of parameters, it is first required to obtain the prior distributions of the parameters. As stated earlier, the available ADT data can be used to obtain the parameter's prior distributions. The ADT data is modeled with a monotonic gamma process to obtain parameter estimates at each accelerated condition. Then the acceleration factor concept has been used to convert those parameters into normal operating conditions. Finally, best-fitted distribution is obtained for parameter estimates of α, β, and c.

Parameter Estimates at Accelerated Conditions
Consider a constant stress ADT to obtain lifetime information for a highly reliable product. Let, represents the ℎ observation of the ℎ sample under the ℎ stress level at the time ; the degradation increment is expressed by ΔY = − ( −1) ; 0 , , and represents the normal stress, maximum stress, and the i th stress level. To obtain the parameter estimates of ̂, ̂, and ̂, the maximum likelihood function for gamma degradation process can be written as:

Parameter Distribution at Normal Operating Conditions
Equation (15) will provide parameter estimates at different accelerated stress level for each sample degradation data. As the posterior parameter analysis is based on field conditions; therefore, these accelerated parameter estimates need to convert into normal operating conditions. The acceleration factor (AF) concept will be used to convert the parameter estimates at normal operating conditions. The acceleration factor can be defined by the ratio of two lifetimes at different stress level as follows: The idea in a gamma process is that due to acceleration, the parameters are changes and subsequently the lifetime. In literature, it assumes that for a gamma process, only shaper parameter α is stress dependable and scale parameter β is independent of stresses (Lawless and Crowder, 2004;Park and Padgett, 2005;Tseng et al., 2009). However, recent studies show that both the gamma parameters are stress-dependent (Rathod et al., 2011;Balakrishnan and Ling, 2014;Limon et al., 2020).
The mathematical expression of gamma parameters and accelerated stresses can be modeled by reaction rate models and their transformed stress relations (Park and Yum, 1997). For multi-stress ADT considering the interaction effects among the stresses, the gamma parameter-stress relation can be established as (Limon et al., 2020): where, represents the r th standardized transformed stress factor at i th stress level, a and b represent the corresponding stress coefficients, and suffix r (1, 2, 3, …R) represents the main stress factors, and suffix q (R+1, R+2, R+3, …Q) represents the interaction effect terms between main stress factors. Now considering the relationship given in Equations (17-18) and the likelihood function for all degradation increments data can be written as: where, ̂= (̂0,̂1, … ,̂, … ,̂,̂, … ,̂,̂0,̂1, … ,̂, … ,̂,̂, … ,̂,). The nonlinear constant c does not depend on the stress level and also there is no meaningful physical explanation of the relation between stress and nonlinear constant. However, it changes due to the random effects of sample units and their degradation behavior. Based on both gamma parameter stress-dependent assumption and parameter stress relations, the parameters can be converted into normal operating conditions by the following equations: ; 0 = ( ,0 ) ; 0 = Once all the parameter estimates α, β, and c are converted at normal operating conditions, the next step is to obtain the best-fitted distribution of these parameter estimates that can be treated as the prior distribution for the Bayesian model. As nonnegative parameters, Exponential, Lognormal, Weibull, and Gamma distributions are considered as potential candidates for the prior distribution. The Anderson-Darling statistic will be used to find the best fitted prior distribution (Thas and Ottoy, 2003).

Numerical Case Example
The light-emitting diodes (LED) are taken as an example to demonstrated the proposed method. The accelerated degradation data of LEDs are listed in Table 1 (Eghbali, 2000;Liao and Elsayed, 2006). The light intensity of LEDs is considered as quality characteristics that indirectly measure the degradation of LED performance. The temperature and electric current are considered accelerating stresses. Four combinations of constant accelerating stresses with five sample units for stress combination data are available for analysis. The light intensity (lux) as a degradation measure is observed in a 250 hours interval and failure is considered while light intensity drops below 50% of its initial value (Wang and Chu, 2012). The normal operating temperature is set 40 0 C and electric current at 10 mA.  Figure 1 illustrates the degradation path of each sample unit at different stress levels. The LED's degradation path indicates the nonlinear nature and monotonicity of the degradation behavior. This justifies the use of the nonlinear gamma process for degradation modeling.

Parameter's Prior Distributions
Equation (15) is used to obtain the maximum likelihood estimates (MLE) of model parameters for each sample unit as given in Table 2. These parameters are for accelerated stress levels which need to convert at the normal operating level. Equation (19) is used to estimate the overall coefficients considering all the stress level degradation data together and estimates are (̂0,̂1,̂2,̂0,̂1,̂2,) = (7.3698, -7.159, -0.1318, 7.4505, -7.8276, -1.7978, 0.4748). Now using these overall coefficients, parameter estimates of α, β, and c can be calculated at any stress levels.
This also provides the opportunity to calculate the parameters and total acceleration factors for different stress levels listed in Table 3.
Once the acceleration factor for each parameter is obtained, the parameter estimates are derived at normal operating conditions that are set as 40 0 C and 10 mA. The parameter estimates at normal operating conditions are listed in Table 4.    The Anderson-Darling (AD) statistic is used to choose the best fitting prior distribution. Table 5 provides AD statistic for four different continuous distributions. For our case example, it is concluded that all three parameter's prior distribution can be best fitted by the gamma distribution with smaller AD statistic. The prior distribution parameters are ̂~(0.9727, 0.0001564), ̂~( 2.890, 0.0006969), and ̂(9.996, 21.725). It is noted that our approach is equally applicable and remains the same even each parameter having different prior distributions.

Parameter's Posterior and RUL
In the presence of field operating degradation data, the Bayesian method is applied to update the model parameters. Since the field degradation is not available, we use simulated degradation data as given in Table 6. shown in Figure 2. For each parameter, 10,000 MCMC samples were used, however, for more reliable information and estimation accuracy, the first 500 samples are eliminated. The posterior expected estimates ( |∆ ), ( |∆ ), and ( |∆ ) for the fifth degradation data point are listed in Table 7.  In Bayesian inference, once new data is available, the previous posterior distribution is considered as prior information to provide posterior estimates for the next time period. The parameter's expected posterior mean for all field degradation data are illustrated in Figure 3. Once the posterior means are obtained, the RUL is estimated by using the Equation.
(7). From Figure 3, it is clear that posterior means of parameters are constantly updated based on the current degradation behavior. Therefore, the RUL estimates are also based on the updated degradation behavior of the product. Figure 4 shows the comparative RUL estimates by Bayesian updated method and traditional without Bayesian updated approach. In traditional without Bayesian RUL estimation, the same initial gamma parameters are used to estimate the RUL for each newly available degradation data.
On the other hand, in the Bayesian method, gamma parameters are first updated for each newly available degradation data and then RUL estimated using the updated parameters. This gamma parameter update ensures to capture the degradation behavior of products based on the most recent information. The Bayesian method also provides an opportunity to understand the individual product's degradation propagation and variation due to random effects. Therefore, the Bayesian method with an integrated ADT approach provides more realistic RUL estimates compared to without Bayesian updated method because updated parameters reflect the product's degradation behavior more accurately.
As the degradation behavior is probabilistic nature, thus degradation propagation is varying randomly. In Table 6, it is observed that degradation increment is varying even the time interval is the same. The increment for the last interval (fifth observation) is less than the previous degradation increment. The change of degradation behavior is well captured into the Bayesian method by changing the gamma parameters. The lower degradation increments were also indicative of a slower degradation rate and hence, the RUL estimate for fifth observation slightly higher than the previous RUL estimate (see Figure 4). As traditional without Bayesian method does not consider the parameter updates based on product's degradation behavior, thus it estimates usual lower RUL at fifth observation than the previous. To better understand the RUL estimates, the density distribution of each updated degradation point is illustrated in Figure 5. The density estimates will provide an insight about the range of values of RUL estimates. This RUL estimates and density for updated degradation behavior can be further utilized for condition-based maintenance planning.

Conclusion
In this work, the RUL prediction of highly reliable products has been proposed using the integration of ADT data and the Bayesian inference method. The monotonic gamma process has been adopted to model the degradation behavior because of its nonnegative increment property. We have also considered the multi-stress condition with the interaction between stresses. A nonstationary gamma process is considered in this work that is well equipped to capture the nonlinear degradation behavior. Ignoring this interaction and nonlinear degradation phenomena could lead to inaccurate degradation modeling and RUL estimation as well. Further, this approach is equally applicable for linear degradation behavior with a single stress scenario.
The acceleration factor has been used to convert the accelerated parameter into normal operating conditions. As it is considered that both gamma parameters are stress-dependent thus acceleration factor for each parameter is separately estimated. It is observed that the stress factors affect each parameter's acceleration differently. The random effect of each parameter also considered. The case example results indicate the superiority of the Bayesian-based RUL prediction compared to the non-Bayesian method. Because the Bayesian-based method provides more accurate product performance behavior in real-time by updating the model parameters based on the most recent degradation information. Also, as more real-time degradation data is available, the accuracy of the model parameters as well as RUL estimates increases with a tighter density of RUL (see Figure 5).
As the gamma stochastic process is considered for degradation modeling therefore, this approach is only applicable for physical systems, products, or components showing monotonic degradation behavior only. Any natural degradation type phenomena, for example, the water or glacier level of