Remaining Useful Life Prediction for Two-Phase Nonlinear Degrading Systems with Three-Source Variability

Recently, the estimation of remaining useful life (RUL) for two-phase nonlinear degrading devices has shown rising momentum for ensuring their safe and reliable operation. The degradation processes of such systems are influenced by the temporal variability, unit-to-unit variability, and measurement variability jointly. However, current studies only consider these three sources of variability partially. To this end, this paper presents a two-phase nonlinear degradation model with three-source variability based on the nonlinear Wiener process. Then, the approximate analytical solution of the RUL with three-source variability is derived under the concept of the first passage time (FPT). For better implementation, the offline model parameter estimation is conducted by the maximum likelihood estimation (MLE), and the Bayesian rule in conjunction with the Kalman filtering (KF) algorithm are utilized for the online model updating. Finally, the effectiveness of the proposed approach is validated through a numerical example and a practical case study of the capacitor degradation data. The results show that it is necessary to incorporate three-source variability simultaneously into the RUL prediction of the two-phase nonlinear degrading systems.


Introduction
With the rapid development of technology, the equipment in the fields of aerospace, high-speed rail, and ship manufacturing tend to be large-scale, complex, sophisticated, and long-life, which puts forward higher requirements for system reliability [1,2].Recently, owing to the important role in realizing predictive maintenance, improving system reliability, reducing maintenance costs, and avoiding catastrophic eventualities, prognostics and health management (PHM) has obtained extensive attention both in academia and industry [3][4][5].As an essential part of PHM, the remaining useful life (RUL) is defined as the length from the current time to failure [6].RUL prediction results provide sufficient information support for decision-making and maintenance scheduling.The performance of RUL prediction primarily relies on the degradation modeling approaches.Generally, the RUL prediction approaches can be systematically classified into model-based methods, data-driven methods, and hybrid methods [7].More detailed discussions about those three types of methods can be found in [8][9][10].As one category of data-driven methods, the stochastic model-based methods have been widely used attributable to their great potential in characterizing the stochastic dynamic degradation process and providing the probability distribution of RUL [11,12].Stochastic process models mainly include the Wiener process, the Gamma process, and the Inverse Gaussian process [13].Among them, the Wiener process model has attached increasing interest due to its explicit statistical interpretation and illustrious mathematical properties in describing the non-monotonic degradation process [14][15][16].
In practice, the degradation process of many engineering systems occurs in a stochastic way.Thus, the RUL is a random variable, which is difficult to estimate with certainty [6].
Therefore, in the Wiener process-based RUL prediction, the commonly used approach is to estimate the probability density function (PDF) of the system RUL by modeling the sensed degradation data.Due to some unobservable internal and external factors, the degradation process of a system is often influenced by multiple sources of variability, which leads to uncertainty in RUL prediction.Typically, the degradation process is affected by three main sources of variability: temporal variability, unit-to-unit variability, and measurement variability [17].To be specific, temporal variability depicts the inherent stochasticity of the degradation process over time, unit-to-unit variability indicates the heterogeneity of degradation trajectories among different units, and measurement variability describes the inevitable measurement error caused by noise in condition monitoring (CM) data [18].Hence, to improve the accuracy of RUL prediction, the multi-source variability should be incorporated into the degradation modeling process simultaneously [19].Over the past decade, various published works have described such variability in degradation modeling [20][21][22][23][24].However, most studies only focused on the RUL estimation with one or two sources of variability, which limits their adaptivity.
Recently, for RUL prediction with three-source variability, Si et al. [19] presented a linear Wiener process-based degradation model considering the three-source variability simultaneously and derived the analytical expressions of RUL distribution.In this work, the random drift coefficient and the underlying degradation state were updated jointly via a state-space model and the Kalman filtering (KF) technique.To deal with nonlinear patterns, Zheng et al. [25] extended the work in [19] and constructed a degradation model based on the nonlinear Wiener process, integrating both nonlinearity and three-source variability into the RUL prediction.For newly developed small sample systems, to address the lack of historical data and prior information, Wang et al. [26] proposed an adaptive RUL estimation method with three-source variability based on the expectation maximization (EM) algorithm.Aiming at the unbalanced historical degradation measurements, Yu et al. [27] established a nonlinear-drift-driven prognostic model considering three-source variability and formulated the approximate analytical expressions of the RUL for both online and offline estimation scenarios.On this basis, the MLE method in conjunction with a downsampling strategy was utilized for parameter estimation and the underflow issue averting in this study.It is noteworthy that the aforementioned linear or nonlinear Wiener processbased prognostic methods with three-source variability only focused on the single-phase degradation cases.However, in practical engineering, owing to changes in the external dynamic environment and internal degradation mechanisms, the degradation trajectories of products such as batteries [28], liquid coupling devices [29], and light emitting diodes [30] exhibit obvious two-phase characteristics with evident inflection points.Under these circumstances, the traditional single-phase model is insufficient to track the dynamics of such degradation processes.Therefore, formulating two-phase degradation models to guarantee the accuracy of RUL prediction is practical and necessary.
In recent years, researchers have reported various RUL prediction approaches considering multiple sources of variability for two-phase degradation processes.Zhang et al. [31] proposed a two-phase linear degradation model based on the Wiener process and derived the analytical forms of the RUL estimation.Chen et al. [32] extended the work in [31] and presented an extreme learning machine (ELM) algorithm to adaptively detect the random changing time of the two-phase degradation trajectory.For two-phase nonlinear degradation patterns, Lin et al. [33] formulated a nonlinear Wiener process-based degradation model and obtained an approximate analytical solution of the lifetime estimation.Hu et al. [34] constructed two RUL prediction models based on the nonlinear time-scale transformation model and the stochastic process model with purely time-dependent parameters.The commonality of the aforementioned works is that the unit-to-unit variability is considered.However, a common deficiency shared by these models is that the measurement variability is omitted.Taking into account the variability of measurement, Wang et al. [35] proposed a change-point Wiener process model with measurement errors to fit the twophase deteriorated OLEDs and adopt the hierarchical Bayesian method to implement parameters estimation, but the number of change-points was assumed to be known a priori.Guan et al. [36] established a two-phase degradation model with measurement errors and random drift terms for RUL estimation, however, the changing time of the model was fixed, and the nonlinearity was not considered in this work.To achieve RUL prediction in the case of imperfect prior degradation information, Chai et al. [37] proposed a linearnonlinear two-phase Wiener process-based degradation model considering measurement errors and the unknown degradation state at the change point simultaneously, which is a pioneering work.Nevertheless, the authors only described the unit-to-unit variability in parameter estimation, while ignoring the random effects of drift coefficients in the PDF expressions of RUL.In addition, the first phase of the proposed model is linear, which limits the method's adaptability.
The above-mentioned research achieved promising results in the two-phase linear or nonlinear degradation modeling with multiple sources of variability.However, most of the existing works about two-phase patterns only considered one or two sources of variability.In contrast, the research on RUL prediction using two-phase degradation models taking into account three-source variability simultaneously is very limited.Furthermore, degrading systems exhibiting two-phase nonlinear degradation patterns are extensively encountered in practice [38,39].Therefore, determining how to achieve two-phase nonlinear degradation modeling and RUL estimation considering three-source variability simultaneously is a compelling practical issue, which motivates our research in this paper.
To achieve this goal, several issues still need to be further investigated.First, the degradation state at the changing point of the two-phase degradation path is an unknown random variable until the changing point appears, and it is related to the changing time as well as the degradation rate of the first phase [31].Thus, it is natural to incorporate the uncertainty of the degradation state at the changing point into RUL prediction.Currently, the state transition probability function is often applied to solve this problem in most studies [33,34].However, these works did not consider the random degradation state at the changing point and the three-source variability simultaneously in the PDF of RUL.Second, it is noted that the changing points and the degradation rates exist difference for devices of the same batch owing to the unit-to-unit variability, and the degradation states are underlying due to the measurement noise.Therefore, determining how to realize parameter identification and real-time degradation state updating in the context of three-source variability poses an interesting challenge.
Motivated by these practical issues, the major contributions of this paper lie in the following aspects: (1) A two-phase nonlinear Wiener process-based degradation model is formulated, where the drift coefficient and the measurement error of each phase are assumed to be random variables to describe the three-source variability.(2) Taking into account the nonlinearity, the uncertainty of the degradation state at the changing point as well as the three-source variability simultaneously, the approximate analytical expressions of RUL estimation are derived under the concept of the first passage time (FPT).(3) Based on the historical degradation observations of multiple units from the same batch, the offline parameter estimation is conducted by the MLE method.Subsequently, with the newly obtained degradation data of the certain operating device, the random drift coefficients and the underlying degradation states are real-time updated by combining the Bayesian rule, state-space model, and KF technique.(4) Finally, a numerical simulation and a practical case study about the degradation data of the high-voltage pulse capacitors are implemented to verify the effectiveness and applicability of the proposed approach.
The remainder of this paper is organized as follows.In Section 2, the two-phase nonlinear Wiener process-based degradation model with three-source variability is formulated.Section 3 presents the approximate analytical solutions of the RUL estimation considering three-source variability.The model parameter estimation is conducted in Section 4.
Section 5 shows the implementation details and results analysis of the experiments.The conclusions are summarized in Section 6.

Degradation Modeling Description
To describe the degradation process with two-phase nonlinearity and three-source variability, the nonlinear Wiener process is employed in our work.Inspired by the single-phase nonlinear degradation model considering three-source variability discussed in Ref. [25] and the two-phase nonlinear degradation model with unit-to-unit variability presented by Refs.[33,34], a general nonlinear Wiener process-based degradation model can be described as, where X(t) denotes the actual degradation state at time t, x 0 is the initial value.For simplifying later calculation, we assume x 0 = 0. τ represents the changing time, thus, x τ is the degradation state at the changing time.λ 1 , λ 2 and σ 1 , σ 2 represent the drift and diffusion coefficients of each phase, respectively.Meanwhile, µ 1 (ρ; ϑ 1 ) and µ 2 (ρ − τ; ϑ 2 ) are nonlinear functions with time t and unknown parameters ϑ 1 , ϑ 2 , which are used to describe the nonlinear characteristics of X(t).B(t) denotes the standard Brown motion.For simplicity, we assume that the abovementioned parameters are independent of phase, meanwhile, the two phases are independent of each other [33].
Owing to the dynamics of {B(t), t ≥ 0}, the temporal variability of Equation (1) could be described.Such a modeling approach has been widely applied to characterize the stochastic degradation process of dynamic systems [27].Further, we assume that ) to describe the unit-to-unit variability, and they are statistically independent of {B(t), t ≥ 0}.Moreover, ϑ 1 , ϑ 2 , σ 1 , σ 2 are fixed parameters used to characterize the common degradation features of all systems from the same batch.In addition, due to the influence of the dynamic environment and the non-ideal instruments, the obtained observations inevitably contain noise.Therefore, to characterize the measurement variability between observed values and the true degradation states, the degradation measurement process {Y(t), t ≥ 0} is formulated as follows [25,37], where ε 1 , ε 2 are the measurement errors of each phase, and assumed to be independent and identically distributed (i.i.d.) with ε 1 ∼ N(0, σ 2 1ε ) and ε 2 ∼ N(0, σ 2 2ε ), respectively.It is further assumed that the measurement errors are independent of X(t).
Generally, the lifetime is usually defined as the FPT when the degradation process exceeds the preset failure threshold w [40].Based on the concept of FPT, the lifetime T of a degrading system can be defined as [19], Then, the expression of RUL at the current time t k can be defined as [25]: where l k represents the time from t k to the failure time, L k is the RUL with conditional PDF Initially, we only consider the temporal variability in lifetime estimation.In this case, the degradation process {X(t), t ≥ 0} described by Equation (1) could be directly observed.If the changing time and all parameters in Equation ( 1) are known fixed values, the PDF of the lifetime T, which only considers the temporal variability could be summarized as follows [33,34], (1) if 0 < t ≤ τ, (2 It is noticeable that the randomness of all parameters is neglected in Equations ( 5) and (6).For the two-phase degradation process with a random changing time, the degradation state at the changing point is usually unknown and could be obtained through a transition probability density from x 0 to x τ under an absorbable boundary w [31].Thus, taking into account the unit-to-unit variability and the randomness of the degradation state at the changing point jointly, the PDF of the lifetime T could be obtained via the law of total probability as follows, where h τ (x τ λ 1p , σ 1p ) is the transition probability density function considering the randomness of the first phase model, i.e., λ 1 ∼ N(λ 1p , σ 2 1p ).Then, according to the relationship between the lifetime and RUL [19,25], the PDF of RUL that considers both temporal variability and unit-to-unit variability could be derived as follows, Theorem 1.For the two-phase nonlinear degradation process in Equation ( 1) and the definition of RUL proposed in Equation ( 4), given the actual degradation state x k at the current time t k and ), the PDF of RUL with a certain changing time τ considering the temporal variability, unit-to-unit variability and the random degradation state at the changing point jointly can be formulated as, Case 1: The current time t k is smaller than the changing time τ (i.e., t k < τ) where S = S 1 − S 2 , T = T 1 − T 2 , and Case 2: The current time t k is larger than the changing time τ (i.e., t k ≥ τ).
Proof.See Appendix A.
It is worth noting that the aforementioned results assume that the degradation state x k could be directly and accurately observed.However, measurement errors are inevitable in practice.Therefore, to consider the impact of measurement errors, the measurement variability should be incorporated into the RUL estimation results of Theorem 1.

RUL Estimation Considering Three-Source Variability
Owing to the influence of measurement errors, the true degradation state x k at the current time t k is unknown.Therefore, we assume that x k ∼ N( xk|k , P k|k ) and x k could be updated based on the degradation observations Y 1:k = {y 1 , y 2 , • • • , y k } via the updating procedure proposed in Section 4.3.Then, the PDF of RUL with three-source variability could be derived based on x k ∼ N( xk|k , P k|k ) and Theorem 1.
Theorem 2. For the two-phase nonlinear degradation process in Equation ( 1) and the definition of RUL proposed in Equation ( 4), given the degradation observations Y 1:k up to the current time t k , the PDF of RUL with a certain changing time τ considering the random degradation state at the changing point and the three-source variability simultaneously can be formulated as follows, Case 1: where where and where In addition, it is worth mentioning that the changing time of the above results is a fixed value.In the case of three-source variability, considering the randomness of the changing point τ, the PDF of RUL could be derived as follows [31], where p(τ) represents the PDF of the changing time τ.Equation ( 13) could be solved by some numerical methods.

Model Parameter and Degradation State Estimation
Before conducting RUL prediction, the unknown model parameters should be identified first.Considering the common characteristics of devices from the same batch and the individual features of the specific operating device, the offline and online methods are adopted jointly to identify the model parameters.Firstly, based on the historical observations of devices within the same batch, the MLE method is used for offline parameter estimation.Then, the changing point detection method is constructed on this basis.Finally, according to the real-time monitoring data of one specific operating device, the online updating procedures of the model parameters are realized via the Bayesian rule and KF algorithm.

Offline Parameter Estimation
In this subsection, we mainly focus on the estimation of the unknown model parameters, which are denoted as ) are random variables that describe the unit-to-unit variability, and ϑ 1 , σ 1 , σ 1ε , ϑ 2 , σ 2 , σ 2ε are deterministic parameters that depict the common degradation features of devices from the same batch.It is noted that τ will be estimated via the changing point detection method in Section 4.2.
It is assumed that the historical observations of N devices within the same batch are known, i.e., Y 1: , and the corresponding actual degradation state data is denoted as m n } and m n denotes the available number of the observations.We define that the measured increment of the n-th device is For simplicity, the time interval is assumed to be a constant, i.e., ∆t = t n,j − t n,j−1 .Suppose τ n denotes the changing time of the n-th device, and it is assumed that τ n = τ n /∆t ∈ {0, 1, . . . ,m n } denotes the changing point location for simplifying the later calculation.Scilicet, the changing time τ n only appears at the measurement time {t n,0 , t n,1 , • • • , t n,m n }.As a result, ∆Y 1n = ∆y n,1 , ∆y n,2 , • • • , ∆y n, τ n represents the measured increments in the first phase, and ∆Y 2n = ∆y n, τ n +1 , ∆y n, τ n +2 , • • • , ∆y n,m n represents the measured increments in the second phase.
According to the definition in Equation ( 1) and the properties of the Wiener process, the measured increments ∆Y n = {∆y n,1 , ∆y n,2 , • • • , ∆y n,m n } follow the multivariate normal distribution with expectation and covariance matrix as follows, where Thus, the log-likelihood function of Θ could be expressed as follows, To facilitate calculation, let , and Then, the log-likelihood function in Equation ( 16) could be written as, Under the framework of the MLE algorithm, taking the first partial derivatives of ln(Θ|Y 1:N ) with respect to λ 1p and σ 2 1p , respectively.Then, the MLE of λ 1p , σ 2 1p could be obtained by setting these two derivatives to zero.

Changing Point Detection
For the two-phase nonlinear degradation process, changing point detection is crucial for parameter estimation and RUL prediction.It is assumed that the changing time is a random variable and follows the normal distribution, i.e., τ ∼ N(µ τ , σ 2 τ ).Then, for a specific operating device, the log-likelihood function of τ can be formulated as follows, where Similar to the offline parameter estimation process, the MLE of the deterministic parameters and the expressions of λ1p , σ2 1p , λ2p , σ2 2p could be derived.Then, by substituting these MLE values and expressions into Equation ( 22), a log-likelihood function that is only related to the changing point location τ n can be formulated.Thus, the optimal changing time τn of the n-th device could be expressed as follows, On this basis, maximizing ln L( τ n |X n ) by enumerating all possible values of τ n , the optimal changing time τn could be obtained.
In addition, for a specific operating device, if t k > τn at the current time t k , the changing time has appeared and is a certain value.However, if t k < τn , which means that the changing point has not appeared and is an unknown random variable.Consequently, we need to update its distribution.In this case, the above-estimated value τn can be treated as the observations of the changing time.Then, the mean and variance of τ could be formulated as [33],

Online Implicit State Updating
For the RUL estimation with three-source variability, the random drift coefficients λ 1 , λ 2 and the underlying degradation state x k need to be updated jointly based on the observations.In practice, the posterior estimates of the drift coefficients have little effect with x k [24].Without loss of generality, it is assumed that they are independent of each other.On this basis, the update mechanism of implicit states could be conducted through two steps: firstly, the random drift coefficients are updated by the Bayesian method, and then the KF algorithm is adopted to update x k .
As to a specific operating device, if the current time is t k , the degradation observation is denoted as Y 1:k = {y 1 , y 2 , • • • , y k }, and the corresponding actual degradation state is denoted as It is defined that the observation increment is ∆y j = y j − y j−1 , and the time interval is ∆t

dρ, and we have
To update the drift coefficients, let λ 1p,0 , σ 1p,0 and λ 2p,0 , σ 2p,0 attained in Section 4.1 represent the prior values of λ 1 and λ 2 , respectively.Then, according to the Bayesian rule [24], the following results could be obtained, where where ,0 + 1 To update the underlying degradation state, i.e., x k ∼ N( xk|k , P k|k ), the state and measurement formulas in Equations ( 1) and ( 2) should be converted into discrete time equations.In this case, a general state space model is constructed as follows, where s = 1, 2 denotes the phase.If ).In addition, it is assumed that {v k } k≥1 and {ε sk } k≥1 are independent of each other.
Iterating the above KF process step by step, the posterior distribution of x k could be solved analytically.Based on the implicit state updating results, the RUL estimation could be conducted.

Case Study
To verify the feasibility of the proposed method, a numerical simulation and a practical example of the high-voltage pulse capacitor are provided.For performance evaluation, the prediction results of the proposed method are compared with the two-phase nonlinear degradation model considering unit-to-unit variability (Lin's method) [33], the two-phase linear-nonlinear degradation model with measurement errors (Chai's method) [37], and the traditional single-phase nonlinear degradation model with three-source variability (Zheng's method) [25].The implementation details of the experiment are as follows.

Numerical Example
In this subsection, we focus on validating the effectiveness of our method for parameter identification and RUL estimation.For illustration purposes, the following power model is applied to define the nonlinear integral term in Equation (1), i.e., This kind of power function has been widely used in existing literature [22].Similar to the work of Feng et al. [20] we adopt the state space model to generate the simulation data and the parameters are given as: 4, σ 2 = 0.08, σ 2ε = 0.4 and the distribution parameters of the random changing time, i.e., µ τ = 50, σ τ = 2. On this basis, we generate 100 sets of degradation trajectories with random changing time, drift parameters, and measurement errors.Figure 1 shows several typical degradation paths of the sample, and it is observable that the degradation paths exhibit obvious two-phase nonlinear deteriorating patterns.Next, the unknown model parameters are identified based on the changing time detection and parameter estimation methods proposed in Section 4. The parameter estimation results with different sample sizes are detailed in Table 1.It can be found from Table 1 that the obtained results could gradually approach the true values as the size of the degradation paths is increased.In addition, the histogram of the detected changing time based on 100 sample paths is presented in Figure 2 and it manifests that the estimated changing points are concentrated around t = 50, which are very close to the true value of µ τ .Thus, the effectiveness of the offline parameter identification method is demonstrated.For the online parameter and the underlying degradation state updating, we chose one specific degradation path from the above-presented 100 degradation trajectories and the true drift parameters are λ 1 = 1.1133, λ 2 = 1.4690, and τ = 51.4 as shown in Figure 3.
Then the changing point detection procedure is applied for the online path, and the estimated changing time is 51.3, which is very close to the true value.Meanwhile, the parameter estimation results of n = 10 size in Table 1 are used as the prior information, when newly observed data are coming, the hyper-parameters of the drift coefficients and the underlying degradation state could be updated via the methods proposed in Section 4.3.Figure 4a shows the parameter updating process, it could be observed that the parameter updating curves could gradually approach the true values.Furthermore, it is clear from Figure 4a that as the observations accumulate, the values of σ 1p and σ 2p are gradually decreasing, which means the uncertainty of estimation is reduced.In addition, Figure 4b presents the comparison of the predicted underlying degradation path and the observed degradation path.It is clear from Figure 4b that our degradation state updating method could track the actual observed degradation state well.Based on the results of the parameter identification process, the PDFs of RUL estimation at different observed time points could be obtained, as shown in Figure 5a.It is worth noting that the preset failure threshold w = 570, thus the actual lifetime is T = 80.It can be observed from Figure 5a that the PDFs of RUL estimated by our model are closely distributed around the actual RUL values, verifying that our proposed method could effectively estimate the RUL of the two-phase nonlinear degradation process.For comparative purposes, we further obtained the mean RUL prediction results of our method, Lin's method [33], Chai's method [37], and Zheng's method [25], as shown in Figure 5b.It is observable from Figure 5b that the mean predicted RUL of our method is more accurate compared to the other three methods.Additionally, the mean RUL prediction curve of Lin's method [33] is relatively stable due to the degradation model they constructed considering two-phase nonlinearity.However, since they ignored the measurement variability, the deviations of their mean RUL prediction curve from the actual values are bigger than our method.Furthermore, the mean RUL curve of Chai's method [37] in the first phase has large bias, but it gradually becomes accurate in the second phase.Because Chai's method [37] only consider the nonlinearity of the second phase.In addition, the deviations between the mean RUL prediction curve of Zheng's method [25] and the actual RUL are the largest among the four methods in the first phase, especially near the changing point where the deviation is more pronounced.It is mainly because the degradation model they constructed is single-phase, which is insufficient to describe the two-phase nonlinear degradation process effectively.It is interesting to note that after the changing point has appeared, the deviations of Zheng's method [25] gradually decrease as the observations accumulate.The reason for this phenomenon is that the degradation trajectory is equivalent to a single-phase nonlinear degradation process after the changing point appears, thus, Zheng's method [25] that considers three-source variability simultaneously could effectively fit the degradation process.For performance evaluation, two metrics are employed to evaluate the performance of RUL prediction among different methods, including the mean square error (MSE) [22] and the absolute error (AE) [37].The MSE at each observation point could be defined as follows, where l k denotes the actual RUL at t k , and f L k |Y 1:k (l k Y 1:k ) is the corresponding conditional PDF of the RUL estimated by Theorem 2 in Section 3.2.
The second metric is the AE of actual RUL and the estimated RUL at each observation point, which could be represented by where lk denotes the estimated RUL at t k .
It is noted that in both criteria, the smallest value of MSE and AE corresponds to the best RUL prediction result.
The following Figure 6 presents the comparison results of the four methods.It can be seen from Figure 6a that our method provides smaller MSE values than the other three methods.Additionally, the AE values of our method in Figure 6b are also the smallest of the four models.The comparison results of these two criteria indicate that our RUL prediction method could effectively improve the prediction accuracy of RUL.Moreover, it can be observed that the MSE and AE values of Zheng's method [25] are the largest of the four methods near the changing point.Because the single-phase nonlinear degradation model constructed by Zheng's method [25] did not consider the impact of the random degradation state at the changing point on RUL estimation, and cannot estimate the RUL of the two-phase nonlinear degradation process well.In addition, Figure 6c compares the PDFs of the RUL distributions derived from our model with the other three methods at time 40.It is clear from Figure 6c that our method can cover the actual RUL well, which reflects the superiority of our method.Moreover, it can also be observed from Figure 6c that the PDF curve of Zheng's method [25] cannot cover the true value of RUL well, revealing that the bias of RUL predicted by Zheng's method [25] is very large in the first phase, which is consistent with the previous results of MSE and AE.It is noteworthy that the units are omitted in Figures 1-6 since the degradation data are generated via simulation.To quantitatively compare the RUL prediction performance of our method with the other methods, we further employ three evaluation metrics.The first metric is the total MSE (TMSE) [19], which is defined as the sum of the MSE at each observation point over the whole life cycle.Videlicet, if there are m observations, the TMSE could be formulated as, The second metric is the mean absolute error (MAE) [41], which measures the average absolute deviation between the predicted value and the actual value.The MAE could be represented as, where AE k is defined in Equation (30).The third metric is the cumulative relative accuracy (CRA) [27]., which evaluates the relative prediction accuracy of the RUL over time, and could be defined as, For these three quantitative metrics, the smallest value of TMSE and MAE corresponds to the best RUL prediction result, in contrast, a higher CRA value that approaches to 1 indicates higher RUL prediction accuracy.Then, the subsequent Table 2 shows the comparison results of the estimated RUL based on the three metrics mentioned above.It can be found from Table 2 that compared with the other three methods, our method has smaller TMSE and MAE values, whereas a bigger CRA value, which indicates that our method has higher accuracy in RUL prediction.

Practical Application
In this subsection, the data of high-voltage pulse capacitors are utilized to illustrate our approach [39].High-voltage-pulse capacitors can be utilized for electrical energy storing and releasing, which are often encountered in pulse lasers, radar, and particle accelerators [34].It is well known that capacitors and batteries can form hybrid energy storage systems, which are commonly used as energy sources for hybrid electric vehicles [42][43][44].However, the unit-to-unit variability mentioned in this work is defined as the heterogeneity between different units in the same batch of identical devices [33,34].Therefore, if degradation modeling with three-source variability is based on the degradation data of different types of devices, such as capacitors and batteries, the basic premise for the unit-to-unit variability is not satisfied.In addition, if the research object is a system composed of batteries and capacitors, it is another topic of RUL prediction, namely, the RUL estimation for multi-component systems.Currently, our work mainly focuses on the RUL prediction with three-source variability of different devices from the same batch, thus, only the capacitor degradation data is considered in this research.
According to the usage scenario, the high-voltage pulse capacitors have a short service time and a long storage time.Hence, the storage performance is a primary factor that influences the reliability of such capacitors.Capacitance can be used as a health indicator (HI) to evaluate the RUL of the capacitors; Figure 7 shows the degradation test data of five high-voltage pulse capacitors under the storage condition.The capacitance was observed every month and the relative capacitance variability was selected as the HI.Reference [34] revealed that the degradation paths of such capacitors could be modeled based on the two-phase nonlinear method.In addition, Reference [39] mentioned that the capacitance degradation path is affected by temperature and humidity, thus, measurement errors are inevitable.Therefore, it is appropriate to use the degradation data of these capacitors to verify our proposed two-phase nonlinear degradation model.Based on the degradation data of Capacitors 2-5, the offline parameter estimation results are obtained through the parameter identification method proposed in Section 4, and are listed in Table 3.Through the estimated values of b 1 , b 2 and σ 1ε , σ 2ε in Table 3, it can be found that nonlinearity and measurement errors exist in the degradation process of capacitors.It is noticeable that the results of Table 3 are treated as the prior information.For online implementation processes, the degradation data of Capacitors 1 is selected to illustrate the parameter updating and RUL prediction processes.Figure 8 presents the parameter updating results, and it can be found that the updated values of λ 1p , λ 2p are bigger than the offline values in Table 3.The reason is that the degradation path of Capacitors 1 is steeper than the mean paths of Capacitors 2-5, which can be seen in Figure 7.According to the literature [34,39], the failure threshold of capacitor is set as w = 5%.Then, based on the results of model updates, the PDFs of RUL could be obtained as shown in Figure 9a.It is observable from Figure 9a that the estimated values of RUL are almost consistent with the actual values, indicating that our method can effectively estimate the RUL of capacitor degradation data.For comparative purposes, we further presented the estimated RUL of Lin's method [33], Chai's method [37], and Zheng's method [25] as shown graphically in Figure 9b-d, respectively.It can be observed from Figure 9 that in both phases, the PDF curves of the proposed method are sharper compared to Lin's method [33] and Zheng's method [25].Additionally, although our PDF curves are similar in steepness to Chai's method [37], the PDF curves of our method are more compact around the actual RUL and the estimated RUL values are closer to the actual values compared to the other three methods.Thus, from the overall RUL prediction results, our method that considers three-source variability and two-phase nonlinearity has higher RUL prediction accuracy compared to the other three methods.To evaluate the performance of RUL prediction, three metrics, including MSE, AE, and relative error (RE) [19], are adopted for performance evaluation.Among them, the RE of the estimated RUL at t k can be defined as, where l k and lk denotes the estimated and the actual RUL at t k , respectively.Then, we calculated the MSE, AE, and RE of the predicted RUL, as shown in Figure 10.It can be found from Figure 10a that the MES values of our method maintain a relatively low level compared to the other three methods.Additionally, it can be observed from Figure 10b,c that the AE as well as the RE results of the proposed method are smaller than the other three methods, indicating the higher accuracy of our proposed method.Furthermore, it can be seen from Figure 10c that the RE curves of each method gradually increase in the later stage of the degradation process.The main reason is that as time accumulates, the actual value of RUL becomes smaller and smaller.Although the AE value is gradually decreasing over time, its proportion in the actual RUL value is larger than that in the early stage of the degradation process, thus, resulting in an increase trend in the later stage.In addition, it can be seen from Figure 10a that the MSE value of Chai's method [37] is slightly better than the MES value of our method.To investigate the reasons for this phenomenon, we further obtained the PDFs of the RUL estimate in the fifth, sixth, and seventh adjacent months, as shown in Figure 11.It can be found from Figure 11 that the PDF curves of our method at the three time points can cover the actual values of RUL well, and are more compact around the actual RUL, respectively.Moreover, the estimated RUL values of our method at different time points are closer to the actual values than the other three methods.
The key issue lies in the PDF curve of the estimated RUL at the sixth month in Figure 11b.As can be seen from Figure 11b, the deviation between our estimated RUL and the actual value is only slightly smaller than that of Chai's method [37].However, the PDF curve of Chai's method [37] is sharper than ours.According to the definition of MSE in Equation ( 29), the MSE value is jointly affected by the deviation between the estimated value of RUL and the actual value, as well as the corresponding PDF function.The MSE value will be smaller when the PDF of RUL is closely distributed around the actual RUL value [34].Therefore, the MSE value of Chai's method [37] is slightly smaller than ours in the sixth month.In fact, this is acceptable.Generally, as can be seen in Figure 10a, the MSE curve of our proposed method maintains a relatively low level compared to Chai's method [37].Therefore, the MSE result in the sixth month does not affect the conclusion that our RUL prediction accuracy is higher than Chai's method [37].
To quantitatively compare the predicted results based on the capacitor degradation data, we further calculated the TMSE, MAE, and CRA values of RUL estimation, as detailed in Table 4.It is observable from Table 4 that the TMSE of our method is the smallest.Thus, although the MSE value of Chai's method [37] in Figure 10a is smaller than ours in the sixth month, however, in general, the TMSE of our method is still better than Chai's method [37].In addition, among the four methods, our method has the smallest MAE value and the largest CRA value, because we simultaneously consider the two-phase nonlinearity and the three-source variability of the degradation process.The results of the above quantitative metrics indicate that our RUL prediction method has higher accuracy compared to the other three methods, which verifies the effectiveness and superiority of our method.

TMSE MAE CRA
Our method 4.9349 0.1091 0.9685 Lin's method [33] 7.0592 0.3637 0.8862 Chai's method [37] 7.5904 0.6182 0.8647 Zheng's method [25] 16.2556 1.0455 0.7873 In summary, the experimental results show that compared to the other three methods, the MSE, AE, and RE curves of our method maintain a relatively low level, and the RUL prediction results of our method have smaller TMSE, MAE values, whereas larger CRA values.Based on the above results, it can be found that our work can obtain competitive results, which indicates that the accuracy of RUL prediction is better than the other three methods, thus, verifying the feasibility and effectiveness of the proposed method in practical application.Furthermore, for two-phase nonlinear degradation devices, it is necessary to consider the nonlinearity and the three-source variability simultaneously to boost the performance of RUL prediction.

Conclusions
In this study, a two-phase nonlinear Wiener process-based degradation model subject to the three-source variability is formulated for RUL prediction.Based on the FPT concept, the approximate analytical form of the RUL is obtained considering the three-source variability and the random degradation state at the changing point simultaneously.To incorporate the historical observations of the devices within the same batch, the MLE method is adopted to estimate the unknown model parameters.By combining the Bayesian rule and KF algorithm, the drift coefficients and the underlying degradation state are adaptively updated in real-time with newly observed data.Finally, the effectiveness of the proposed method is verified via numerical and practical examples.The quantitative comparison results of MSE, AE, RE, TMSE, MAE, and CRA reveal that the two-phase nonlinear degradation model with three-source variability is more accurate than the existing methods that only partially consider the two-phase nonlinearity and three-source variability.Therefore, in RUL prediction for the two-phase nonlinear degradation devices, the impact of these uncertainties on degradation modeling should be considered simultaneously to improve the accuracy of RUL prediction.However, there are several directions that are worth further study.First, the measurement error is assumed as Gaussian, which may be inadequate when there are outliers in the degradation observations.In this case, the non-Gaussian measurement errors need to be considered, such as t-distribution.Second, the proposed model is formulated under the assumption of the progressive degradation process.However, shocks are often encountered in practice.Therefore, determining how to construct the two-phase nonlinear degradation model considering the interaction between the internal degradation mechanism and external shocks needs to be explored.Third, at present, our research only focuses on the RUL prediction for single-component systems, without considering the impact of correlation between different components on RUL prediction.Therefore, the degradation modeling and RUL prediction of multicomponent systems (such as hybrid energy systems consisting of capacitors and batteries) with three-source variability is also a valuable research direction in the future.

Nomenclature
The following abbreviations and notations are used in this manuscript: Proof.To obtain the RUL estimation result considering the unit-to-unit variability and the random degradation state at the changing point, the PDF of the lifetime should be derived first, which is given by Ref. [33].
For the two-phase nonlinear degradation model proposed in Equation ( 1), if the drift coefficient in each phase is a random variable and follows a normal distribution, i.e., ), the PDF of the lifetime T considering the unit-to-unit variability and the random degradation state at the changing point could be expressed as follows, where and Based on Equation (A1), according to the relationship between the lifetime T and RUL, the PDF of RUL that considers temporal variability, unit-to-unit variability, and the randomness of the degradation state at the changing point could be obtained as follows, Case 1: The current time t k is smaller than the changing time τ (i.e., t k < τ) where S = S 1 − S 2 , T = T 1 − T 2 , and Case 2: The current time t k is larger than the changing time τ (i.e., t k ≥ τ) The result of Equation (A3) is transformed from lifetime T to RUL based on the law of total probability as f T where the result of f T (t|λ 2 , x τ ) can be found in Equation ( 6), and the proof can be found in reference [33].
In this way, the proof has been completed.

Appendix B. Proof of Theorem 2
Proof.
(A) Based on f T (t|λ 1 ) in Equation ( 5) and λ 1 ∼ N(λ 1p , σ 2 1p ), the PDF of the lifetime considering the unit-to-unit variability could be formulated via Lemma A1 as, (B) Then based on the relationship between lifetime and RUL, we can further obtain the PDF of RUL considering the unit-to-unit variability as follows, (C) When considering the three-source variability simultaneously, it is necessary to further incorporate the measurement variability into the PDF of RUL that considers unit-to-unit variability, which could be presented as, (D) Then, based on x k ∼ N( xk|k , P k|k ) and Lemma A1, we can obtain the PDF of RUL with three-source variability when t k ≤ τ and 0 < l k + t k ≤ τ as follows, where, To solve the PDF of RUL considering three-source variability simultaneously, the PDF formulas in S − T of Theorem 1 can be transformed into the following form, where Then, based on Equation (A9) and x k ∼ N( xk|k , P k|k ), the PDF of RUL with threesource variability can be formulated as follows, separately.To facilitate calculation, the following Lemmas A2 and A3 are provided.

Lemma A2 ([45])
, the S formula in Equation (A9) is expanded to S a and S b as follows,   Then, the E x k |Y 1:k [T a ] in Equation (A29) could be formulated as follows,

Figure 2 .
Figure 2. The estimated changing time of 100 degradation paths.

Figure 4 .
Figure 4.The online updating process.(a) The parameter updating.(b) The underlying degradation state updating.

Figure 5 .
Figure 5.The RUL prediction results.(a) PDFs of the RUL.(b) The mean RUL of the four methods.

Figure 6 .
Figure 6.Performance evaluation of the RUL prediction.(a) MSE of the estimated RUL.(b) AE of the estimated RUL. (c) PDFs of the estimated RUL at time 40.

Figure 7 .
Figure 7. Relative capacitance variability of the capacitors with time.

Figure 10 .
Figure 10.Performance evaluation based on capacitance degradation data of Capacitor 1.(a) MSE of the predicted RUL.(b) AE of the predicted RUL.(c) RE of the predicted RUL.

Figure 11 .
Figure 11.Comparison of the estimated RUL at different months.(a) PDFs of the estimated RUL at the fifth month.(b) PDFs of the estimated RUL at the sixth month.(c) PDFs of the estimated RUL at the seventh month.
x 2 k − b T x k × exp − (c T −d T x k )whereD 1T = d 2 T − 2a T b S , D 2T = c T d T − b S b T .(C)Then, the common item before {•} in the formula T of Equation (A9) could be simplified as follows,I 2 2π(t k +l k −τ) 2 (σ 2 α2 +σ 2 β2) On this basis, the formula T in Equation (A9) could be split into T a and T b as follows,T a = exp − D 1T 2b S x 2 k + D 2T (e T x k ) × Φ g T +h T x Based on x k ∼ N( xk|k , P k|k ) and Lemma A3, E x k |Y 1:k [T a ] could be obtained as, E x k |Y 1:k [T a ] = S +D 1T P k|k ) exp − D 1T 2b S x 2 k + D 2T b S x k × (e T x k ) × exp − (x k − xk|k ) 2 2P k|k × Φ g T +h T x k f S dx k (A32)
Whereas if t k > τ, only the observations in the second phase could be applied, i.e., Y τ:k

Table 1 .
Parameter estimation values with different sample sizes.

Table 2 .
Comparison results of RUL prediction based on TMSE, MAE, and CRA.

Table 3 .
The parameter estimation results of the capacitance degradation data.

Table 4 .
Comparison results of RUL prediction based on TMSE, MAE, and CRA for Capacitor 1 degradation data.