[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Next Article in Journal
Outdoor Navigation Assistive System Based on Robust and Real-Time Visual–Auditory Substitution Approach
Previous Article in Journal
Survey on Energy Harvesting for Biomedical Devices: Applications, Challenges and Future Prospects for African Countries
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

School of Mechanical Engineering, Beijing Institute of Technology, Beijing 100081, China
*
Author to whom correspondence should be addressed.
Sensors 2024, 24(1), 165; https://doi.org/10.3390/s24010165
Submission received: 5 December 2023 / Revised: 24 December 2023 / Accepted: 26 December 2023 / Published: 27 December 2023
Figure 1
<p>The simulated degradation paths.</p> ">
Figure 2
<p>The estimated changing time of 100 degradation paths.</p> ">
Figure 3
<p>The online degradation path.</p> ">
Figure 4
<p>The online updating process. (<bold>a</bold>) The parameter updating. (<bold>b</bold>) The underlying degradation state updating.</p> ">
Figure 5
<p>The RUL prediction results. (<bold>a</bold>) PDFs of the RUL. (<bold>b</bold>) The mean RUL of the four methods.</p> ">
Figure 6
<p>Performance evaluation of the RUL prediction. (<bold>a</bold>) MSE of the estimated RUL. (<bold>b</bold>) AE of the estimated RUL. (<bold>c</bold>) PDFs of the estimated RUL at time 40.</p> ">
Figure 7
<p>Relative capacitance variability of the capacitors with time.</p> ">
Figure 8
<p>Online parameter updating processes of Capacitor 1.</p> ">
Figure 9
<p>PDFs of RUL prediction for Capacitor 1. (<bold>a</bold>) Our method. (<bold>b</bold>) Lin’s method. (<bold>c</bold>) Chai’ method. (<bold>d</bold>) Zheng’s method.</p> ">
Figure 10
<p>Performance evaluation based on capacitance degradation data of Capacitor 1. (<bold>a</bold>) MSE of the predicted RUL. (<bold>b</bold>) AE of the predicted RUL. (<bold>c</bold>) RE of the predicted RUL.</p> ">
Figure 11
<p>Comparison of the estimated RUL at different months. (<bold>a</bold>) PDFs of the estimated RUL at the fifth month. (<bold>b</bold>) PDFs of the estimated RUL at the sixth month. (<bold>c</bold>) PDFs of the estimated RUL at the seventh month.</p> ">
Versions Notes

Abstract

:
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.

1. 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 down-sampling strategy was utilized for parameter estimation and the underflow issue averting in this study. It is noteworthy that the aforementioned linear or nonlinear Wiener process-based 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 two-phase 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 linear-nonlinear 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.

2. 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,
X ( t ) = x 0 + λ 1 0 t μ 1 ( ρ ; ϑ 1 ) d ρ + σ 1 B ( t ) , 0 < t τ x τ + λ 2 τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ + σ 2 B ( t τ ) , t > τ
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 λ 1 ~ N ( μ 1 p , σ 1 p 2 ) and λ 2 ~ N ( μ 2 p , σ 2 p 2 ) 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],
Y ( t ) = X ( t ) + ε 1 , 0 < t τ X ( t ) + ε 2 , t > τ
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 , σ 1 ε 2 ) 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],
T   = inf t : X ( t ) w | X ( 0 ) < w
Then, the expression of RUL at the current time t k can be defined as [25]:
L k = inf l k > 0 : X ( t k + l k ) w
where l k represents the time from t k to the failure time, L k is the RUL with conditional PDF f L k | Y 1 : k ( l k | Y 1 : k ) , and Y 1 : k is the observations up to t k , i.e., Y 1 : k = y 1 , y 2 , , y k .

3. RUL Estimation with Three-Source Variability

3.1. RUL Estimation with Temporal Variability and Unit-to-Unit Variability

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 τ ,
f T ( t | λ 1 ) w x 0 λ 1 0 t μ 1 ( ρ ; ϑ 1 ) d ρ t μ 1 ( t ; ϑ 1 ) 2 π σ 1 2 t 3 exp w x 0 λ 1 0 t μ 1 ( ρ ; ϑ 1 ) d ρ 2 2 σ 1 2 t
(2) if τ < t ,
f T ( t | λ 2 , x τ ) w x τ λ 2 τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ ( t τ ) μ 2 ( t τ ; ϑ 2 ) 2 π σ 2 2 ( t τ ) 3 × exp w x τ λ 2 τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 2 σ 2 2 ( t τ )
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,
f T ( t ) = + f T ( t | λ 1 ) p ( λ 1 ) d λ 1 , 0 < t τ w + f T ( t | λ 2 , x τ ) p ( λ 2 ) h τ ( x τ | λ 1 p , σ 1 p ) d λ 2 d x τ , t > τ  
where h τ ( x τ | λ 1 p , σ 1 p ) is the transition probability density function considering the randomness of the first phase model, i.e., λ 1 ~ N ( λ 1 p , σ 1 p 2 ) .
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  λ 1 ~ N ( λ 1 p , σ 1 p 2 ) , λ 2 ~ N ( λ 2 p , σ 2 p 2 ) , 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 < τ )
f L ( l k ) 1 2 π l k 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k × w x k t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ l k μ 1 ( t k + l k ; ϑ 1 ) × ( w x k ) σ 1 p 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ + λ 1 p σ 1 2 l k t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k × exp w x k λ 1 p t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k , 0 < l k + t k τ S T , τ < l k + t k
where  S = S 1 S 2 ,  T = T 1 T 2 , and
S 1 = r a 2 2 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) exp ( λ α 2 λ β 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) × λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 + σ β 2 2 × Φ λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) + σ α 2 2 σ β 2 2 σ α 2 2 + σ β 2 2 × ϕ λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) S 2 = r b 2 2 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) exp ( λ α 2 λ β 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) × 1 Φ λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) T 1 = I 2 × r a 2 2 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) exp ( λ α 2 λ γ 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) × λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 + σ β 2 2 × Φ λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) + σ α 2 2 σ β 2 2 σ α 2 2 + σ β 2 2 × ϕ λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) T 2 = I 2 × r b 2 2 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) exp ( λ α 2 λ γ 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) × 1 Φ λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) λ α 2 = λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ , λ β 2 = w x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ , λ γ 2 = w + x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 ( w x k ) σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 , σ α 2 2 = σ 2 p 2 τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 + σ 2 2 ( t k + l k τ ) , σ β 2 2 = σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 + σ 1 2 ( τ t k ) , r a 2 = ( t k + l k τ ) σ 2 2 + σ 2 p 2 μ 2 ( t k + l k τ ; ϑ 2 ) τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ α 2 2 , r b 2 = ( t k + l k τ ) σ 2 2 λ a 2 λ 2 p ( t k + l k τ ) μ 2 ( t k + l k τ ; ϑ 2 ) σ α 2 2 , I 2 = exp 2 λ 1 p ( w x k ) σ 1 2 + 2 ( w x k ) 2 σ 1 p 4 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ + ( w x k ) 2 σ 1 p 2 σ 1 2 σ 1 2 + σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 4
Case 2: The current time  t k  is larger than the changing time  τ  (i.e.,  t k τ ).
f L ( l k ) 1 2 π l k 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 l k × w x k t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ l k μ 2 ( t k + l k τ ; ϑ 2 ) × σ 2 p 2 ( w x k ) t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ + λ 2 p σ 2 2 l k σ 2 p 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 + σ 2 2 l k × exp w x k λ 2 p t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 l 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.

3.2. 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 ( x ^ k | 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 ( x ^ k | 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:  t k τ ,
(1) if  t k τ  and  0 < l k + t k τ ,
f L k | Y 1 : k ( l k | Y 1 : k ) M 1 N 1 × w λ 1 p t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ P k | k + υ 1 ( l k ) x ^ k | k P k | k + υ 1 ( l k ) 2 π l k 2 υ 1 2 ( l k ) P k | k + υ 1 ( l k ) × exp w x ^ k | k λ 1 p t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 2 P k | k + υ 1 ( l k )
where
υ 1 ( l k ) = t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k , M 1 = w υ 1 ( l k ) w σ 1 p 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ + λ 1 p σ 1 2 l k t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ l k μ 1 ( t k + l k ; ϑ 1 ) , N 1 = υ 1 ( l k ) σ 1 p 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ l k μ 1 ( t k + l k ; ϑ 1 )
(2) if   τ < l k + t k ,
f L k | Y 1 : k ( l k | Y 1 : k ) S n e w + T n e w
where
S n e w = 1 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) × d S exp a S 2 2 b S 1 + P k | k b S exp 2 a S x ^ k | k + a S 2 P k | k b 3 x ^ k | k 2 2 b S + P k | k Φ e S σ α 2 2 x ^ k | k a S σ α 2 2 P k | k b S + e S P k | k b S f S 2 1 + P k | k b S 2 + σ α 2 4 P k | k 1 + P k | k b S + 1 2 π P k | k + exp ( a S x k ) 2 2 b S × c S x k × exp ( x k x ^ k | k ) 2 2 P k | k × Φ e S σ α 2 2 x k f S d x k + r a 2 σ α 2 2 σ β 2 2 σ α 2 2 + σ β 2 2 × 1 2 π exp 1 2 L 3 S × 1 L 1 S 2 P k | k + 1 exp L 2 S L 1 S x ^ k | k 2 2 L 1 S 2 P k | k + 1 σ α 2 2 = σ 2 p 2 τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 + σ 2 2 ( t k + l k τ ) , σ β 2 2 = σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 + σ 1 2 ( τ t k ) , r a 2 = ( t k + l k τ ) σ 2 2 + σ 2 p 2 μ 2 ( t k + l k τ ; ϑ 2 ) τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ α 2 2 , r b 2 = ( t k + l k τ ) σ 2 2 λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ λ 2 p ( t k + l k τ ) μ 2 ( t k + l k τ ; ϑ 2 ) σ α 2 2 , a S = w λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ , b S = σ α 2 2 + σ β 2 2 , c S = r a 2 × σ α 2 2 σ α 2 2 + σ β 2 2 , d S = r a 2 × w λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ α 2 2 + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ β 2 2 σ α 2 2 + σ β 2 2 r b 2 , e S = w λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ α 2 2 + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ β 2 2 , f S = σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) , L 1 S = f S 2 + b S σ α 2 4 b S f S 2 , L 2 S = a S f S 2 + b S e S σ α 2 4 b S f S 2 f S 2 + b S σ α 2 4 , L 3 S = a S 2 f S 2 + b S e S 2 b S f S 2 a S f S 2 + b S e S σ α 2 2 2 b S f S 2 f S 2 + b S σ α 2 4
and
T n e w = exp 4 w b S λ 1 p σ 1 2 + σ 1 p 2 w c T 2 σ 1 4 2 b S σ 1 4 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) × f T 1 + D 1 T P k | k b S exp 2 D 2 T x ^ k | k + D 2 T 2 P k | k b S D 1 T x ^ k | k 2 2 b S + D 1 T P k | k Φ g T + h T x ^ k | k + D 2 T h T P k | k b S + D 1 T g T P k | k b S f S 2 1 + D 1 T P k | k b S 2 + h T 2 P k | k 1 + D 1 T P k | k b S + 1 2 π P k | k + exp D 1 T 2 b S x k 2 + D 2 T b S x k × e T x k × exp ( x k x ^ k | k ) 2 2 P k | k × Φ g T + h T x k f S d x k + r a 2 σ α 2 2 σ β 2 2 σ α 2 2 + σ β 2 2 × 1 2 π exp 1 2 R 3 T × 1 R 1 T 2 P k | k + 1 exp R 2 T R 1 T x ^ k | k 2 2 R 1 T 2 P k | k + 1
a T = 2 σ 1 p 2 σ 1 4 , b T = 2 λ 1 p σ 1 2 + 4 w σ 1 p 2 σ 1 4 , c T = w + λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ + 2 w σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 d T = σ 1 2 + 2 σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 , D 1 T = d T 2 2 a T b S , D 2 T = c T d T b S b T , e T = r a 2 σ α 2 2 σ 1 2 + 2 σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 σ α 2 2 + σ β 2 2 , f T = r a 2 × w λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 w σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 σ α 2 2 + λ 2 p σ β 2 2 τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ α 2 2 + σ β 2 2 r b 2 , g T = w λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 w σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 σ α 2 2 + λ 2 p σ β 2 2 τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ , h T = σ α 2 2 σ 1 2 + 2 σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 , R 1 T = D 1 T f S 2 + b S h T 2 b S f S 2 , R 2 T = D 2 T f S 2 b S g T h T b S f S 2 D 1 T f S 2 + b S h T 2 , R 3 T = g T 2 f S 2 D 2 T f S 2 b S g T h T 2 b S f S 2 D 1 T f S 2 + b S h T 2
Case 2: if  τ < t k ,
f L k | Y τ ˜ : k ( l k | Y τ ˜ : k ) = M 2 N 2 × w λ 2 p t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ P k | k + υ 2 ( l k ) x ^ k | k P k | k + υ 2 ( l k ) 2 π l k 2 υ 2 2 ( l k ) P k | k + υ 2 ( l k ) × exp w x ^ k | k λ 2 p t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 2 P k | k + υ 2 ( l k )
where
υ 2 ( l k ) = t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 l k , M 2 = w υ 2 ( l k ) w σ 2 p 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ + λ 2 p σ 2 2 l k t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ l k μ 2 ( t k + l k τ ; ϑ 2 ) N 2 = υ 2 ( l k ) σ 2 p 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ l k μ 2 ( t k + l k τ ; ϑ 2 ) .
Proof. 
See Appendix B. □
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],
  f L ( l k ) = t k + f L ( l k | τ ) p ( τ ) d τ
where p ( τ ) represents the PDF of the changing time τ . Equation (13) could be solved by some numerical methods.

4. 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.

4.1. Offline Parameter Estimation

In this subsection, we mainly focus on the estimation of the unknown model parameters, which are denoted as Θ = λ 1 p , σ 1 p , ϑ 1 , σ 1 , σ 1 ε , λ 2 p , σ 2 p , ϑ 2 , σ 2 , σ 2 ε . Among them, λ 1 ~ N ( λ 1 p , σ 1 p 2 ) , λ 2 ~ N ( λ 2 p , σ 2 p 2 ) 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 : N = { Y 1 , Y 2 , , Y N } , and the corresponding actual degradation state data is denoted as X 1 : N = { X 1 , X 2 , , X N } . The observations Y n = { y n , 0 , y n , 1 , , y n , m n } of the n - th device are measured at time { t n , 0 , t n , 1 , , t n , m n } and m n denotes the available number of the observations. We define that the measured increment of the n-th device is Δ y n , j = y n , j y n , j 1 , where j = 1 , 2 , , m n . Thus, Δ Y n = { Δ y n , 1 , Δ y n , 2 , , Δ y n , m n } . 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 1 n = { Δ y n , 1 , Δ y n , 2 , , Δ y n , τ ˜ n } represents the measured increments in the first phase, and Δ Y 2 n = { Δ 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,
  • if 0 < t n , j τ n ,
E Δ Y 1 n = λ 1 p Δ T 1 n , Σ 1 n = σ 1 p 2 Δ T 1 n Δ T 1 n + Ω 1 n
if t n , j > τ n ,
E Δ Y 2 n = λ 2 p Δ T 2 n , Σ 2 n = σ 2 p 2 Δ T 2 n Δ T 2 n + Ω 2 n
where Ω 1 n = σ 1 2 D 1 n + σ 1 ε 2 P 1 n , Ω 2 n = σ 2 2 D 2 n + σ 2 ε 2 P 2 n ,
P 1 n = 1 1 0 0 1 2 1 0 1 2 0 1 0 0 1 2 τ ˜ n × τ ˜ n , P 2 n = 1 1 0 0 1 2 1 0 1 2 0 1 0 0 1 2 ( m n τ ˜ n ) × ( m n τ ˜ n ) Δ T 1 n = 0 t n , 1 μ 1 ( ρ ; ϑ 1 ) d ρ , t n , 1 t n , 2 μ 1 ( ρ ; ϑ 1 ) d ρ , , t n , τ ˜ n 1 t n , τ ˜ n μ 1 ( ρ ; ϑ 1 ) d ρ Δ T 2 n = t n , τ ˜ n t n , τ ˜ n + 1 μ 2 ( ρ τ ; ϑ 2 ) d ρ , t n , τ ˜ n + 1 t n , τ ˜ n + 2 μ 2 ( ρ τ ; ϑ 2 ) d ρ , , t n , m n 1 t n , m n μ 2 ( ρ τ ; ϑ 2 ) d ρ D 1 n = d i a g ( Δ t n , 1 , Δ t n , 2 , , Δ t n , τ n ) , D 2 n = d i a g ( Δ t n , τ ˜ n + 1 , Δ t n , τ ˜ n + 2 , , Δ t n , m n )
Thus, the log-likelihood function of Θ could be expressed as follows,
ln ( Θ | Y 1 : N ) = ln ( 2 π ) 2 n = 1 N τ ˜ n 1 2 n = 1 N ln | Σ 1 n | 1 2 n = 1 N ( Δ Y 1 n λ 1 p Δ T 1 n ) Σ 1 n 1 ( Δ Y 1 n λ 1 p Δ T 1 n ) ln ( 2 π ) 2 n = 1 N ( m n τ ˜ n ) 1 2 n = 1 N ln | Σ 2 n | 1 2 n = 1 N ( Δ Y 2 n λ 2 p Δ T 2 n ) Σ 2 n 1 ( Δ Y 2 n λ 2 p Δ T 2 n )
To facilitate calculation, let σ ˜ 1 2 = σ 1 2 / σ 1 p 2 , σ ˜ 1 ε 2 = σ 1 ε 2 / σ 1 p 2 , Σ ˜ 1 n = Σ 1 n / σ 1 p 2 , and σ ˜ 2 2 = σ 2 2 / σ 2 p 2 , σ ˜ 2 ε 2 = σ 2 ε 2 / σ 2 p 2 , Σ ˜ 2 n = Σ 2 n / σ 2 p 2 . Then, the log-likelihood function in Equation (16) could be written as,
ln ( Θ | Y 1 : N ) = ln ( 2 π ) 2 n = 1 N τ ˜ n ln σ 1 p 2 2 n = 1 N τ ˜ n 1 2 n = 1 N ln | Σ ˜ 1 n | 1 2 σ 1 p 2 n = 1 N ( Δ Y 1 n λ 1 p Δ T 1 n ) Σ ˜ 1 n 1 ( Δ Y 1 n λ 1 p Δ T 1 n ) ln ( 2 π ) 2 n = 1 N ( m n τ ˜ n ) ln σ 2 p 2 2 n = 1 N ( m n τ ˜ n ) 1 2 n = 1 N ln | Σ ˜ 2 n | 1 2 σ 2 p 2 n = 1 N ( Δ Y 2 n λ 2 p Δ T 2 n ) Σ ˜ 2 n 1 ( Δ Y 2 n λ 2 p Δ T 2 n )
Under the framework of the MLE algorithm, taking the first partial derivatives of ln ( Θ | Y 1 : N ) with respect to λ 1 p and σ 1 p 2 , respectively. Then, the MLE of λ 1 p , σ 1 p 2 could be obtained by setting these two derivatives to zero.
λ ^ 1 p = n = 1 N Δ T 1 n Σ ˜ 1 n 1 Δ Y 1 n n = 1 N Δ T 1 n Σ ˜ 1 n 1 Δ T 1 n , σ ^ 1 p 2 = n = 1 N ( Δ Y 1 n λ ^ 1 p Δ T 1 n ) Σ ˜ 1 n 1 ( Δ Y 1 n λ ^ 1 p Δ T 1 n ) n = 1 N τ ˜ n
Similarly, the MLE of λ 2 p , σ 2 p 2 could be obtained as,
λ ^ 2 p = n = 1 N Δ T 2 n Σ ˜ 2 n 1 Δ Y 2 n n = 1 N Δ T 2 n Σ ˜ 2 n 1 Δ T 2 n , σ ^ 2 p 2 = n = 1 N ( Δ Y 2 n λ ^ 2 p Δ T 2 n ) Σ ˜ 2 n 1 ( Δ Y 2 n λ ^ 2 p Δ T 2 n ) n = 1 N ( m n τ ˜ n )
However, the results of λ ^ 1 p , σ ^ 1 p 2 and λ ^ 2 p , σ ^ 2 p 2 still contain the unknown parameters, i.e., ϑ 1 , σ ˜ 1 2 , σ ˜ 1 ε 2 and ϑ 2 , σ ˜ 2 2 , σ ˜ 2 ε 2 . Fortunately, the MLE of these parameters could be calculated based on λ ^ 1 p , σ ^ 1 p 2 and λ ^ 2 p , σ ^ 2 p 2 .
Since the two phases are independent of each other, thus, substituting Equation (18) into Equation (17), the profile likelihood function of ϑ 1 , σ ˜ 1 2 , σ ˜ 1 ε 2 according to λ ^ 1 p , σ ^ 1 p 2 could be expressed as follows,
ln L ( ϑ 1 , σ ˜ 1 2 , σ ˜ 1 ε 2 | Δ Y 1 n ) = ln ( 2 π ) 2 n = 1 N τ ˜ n ln σ ^ 1 p 2 2 n = 1 N τ ˜ n 1 2 n = 1 N τ ˜ n 1 2 n = 1 N ln | Σ ˜ 1 n |
Similarly, by substituting Equation (19) into Equation (17), the profile likelihood function of ϑ 2 , σ ˜ 2 2 , σ ˜ 2 ε 2 according to λ ^ 2 p , σ ^ 2 p 2 could be e formulated as follows.
ln L ( ϑ 2 , σ ˜ 2 2 , σ ˜ 2 ε 2 | Δ Y 2 n ) = ln ( 2 π ) 2 n = 1 N ( m n τ ˜ n ) ln σ ^ 2 p 2 2 n = 1 N ( m n τ ˜ n ) 1 2 n = 1 N ( m n τ ˜ n ) 1 2 n = 1 N ln | Σ ˜ 2 n |
The MLE of ϑ 1 , σ ˜ 1 2 , σ ˜ 1 ε 2 and ϑ 2 , σ ˜ 2 2 , σ ˜ 2 ε 2 could be obtained through a multi-dimensional search, which is implemented by the “fminsearch” function in MATLAB.
Then, substituting the MLE of ϑ 1 , σ ˜ 1 2 , σ ˜ 1 ε 2 and ϑ 2 , σ ˜ 2 2 , σ ˜ 2 ε 2 into Equations (18) and (19), respectively, the final estimates of λ 1 p , σ 1 p 2 and λ 2 p , σ 2 p 2 could be obtained. Then, the MLE of ϑ 1 , σ 1 , σ 1 ε and ϑ 2 , σ 2 , σ 2 ε could be calculated via inverting the relations σ ˜ 1 2 = σ 1 2 / σ 1 p 2 , σ ˜ 1 ε 2 = σ 1 ε 2 / σ 1 p 2 , Σ ˜ 1 n = Σ 1 n / σ 1 p 2 and σ ˜ 2 2 = σ 2 2 / σ 2 p 2 , σ ˜ 2 ε 2 = σ 2 ε 2 / σ 2 p 2 , Σ ˜ 2 n = Σ 2 n / σ 2 p 2 , accordingly.

4.2. 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,
ln L ( Θ | Δ Y n ) = τ ˜ n 2 ln ( 2 π ) 1 2 ln | Σ 1 n | 1 2 ( Δ Y 1 n λ 1 p Δ T 1 n ) Σ 1 n 1 ( Δ Y 1 n λ 1 p Δ T 1 n ) m n τ ˜ n 2 ln ( 2 π ) 1 2 ln | Σ 2 n | 1 2 ( Δ Y 2 n λ 2 p Δ T 2 n ) Σ 2 n 1 ( Δ Y 2 n λ 2 p Δ T 2 n )
where Δ Y n = ( Δ y n , 1 , Δ y n , 2 , , Δ y n , k ) .
Similar to the offline parameter estimation process, the MLE of the deterministic parameters and the expressions of λ ^ 1 p , σ ^ 1 p 2 , λ ^ 2 p , σ ^ 2 p 2 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,
τ ^ n = Δ t × arg max τ ˜ n ln L ( τ ˜ n | Δ Y n )
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],
μ τ = 1 N n = 1 N τ ^ n , σ τ = 1 N n = 1 N ( τ ^ n μ τ ) 2

4.3. 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 X 1 : k = x 1 , x 2 , , x k . It is defined that the observation increment is Δ y j = y j y j 1 , and the time interval is Δ t j = t j t j 1 , j = 1 , , k . If t k τ , all the observations Y 1 : k = y 1 , y 2 , , y k can be used for model updating, thus, Δ Y k = Δ y 1 , Δ y 2 , , Δ y k . Meanwhile, let Δ T k = t k 1 t k μ 1 ( ρ ; ϑ 1 ) d ρ , and we have Δ T k = Δ T 1 , Δ T 2 , , Δ T k . Whereas if t k > τ , only the observations in the second phase could be applied, i.e., Y τ ˜ : k = y τ ˜ , y τ ˜ + 1 , , y k , thus, Δ Y k = Δ y τ ˜ + 1 , Δ y τ ˜ + 2 , , Δ y k . Let Δ T k = t k 1 t k μ 2 ( ρ τ ; ϑ 2 ) d ρ , and consequently, Δ T k = Δ T τ ˜ + 1 , Δ T τ ˜ + 2 , , Δ T k .
To update the drift coefficients, let λ 1 p , 0 , σ 1 p , 0 and λ 2 p , 0 , σ 2 p , 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,
if t k τ ,
λ 1 | Y 1 : k N ( λ 1 p , k , σ 1 p , k 2 )
where
λ 1 p , k = Δ T k Ω k 1 Δ Y k σ 1 p , 0 2 + λ 1 p , 0 Δ T k Ω k 1 Δ T k σ 1 p , 0 2 + 1 , σ 1 p , k 2 = σ 1 p , 0 2 Δ T k Ω k 1 Δ T k σ 1 p , 0 2 + 1
if t k > τ ,
λ 2 | Y τ ˜ : k N ( λ 2 p , k , σ 2 p , k 2 )
where
λ 2 p , k = Δ T k Ω k 1 Δ Y k σ 2 p , 0 2 + λ 2 p , 0 Δ T k Ω k 1 Δ T k σ 2 p , 0 2 + 1 , σ 2 p , k 2 = σ 2 p , 0 2 Δ T k Ω k 1 Δ T k σ 2 p , 0 2 + 1
To update the underlying degradation state, i.e., x k ~ N ( x ^ k | 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,
x k = x k 1 + λ k 1 Δ T k + v k y k = x k + ε s k
where s = 1 , 2 denotes the phase. If t k τ , thus,  s = 1 v k = σ 1 B ( t k ) B ( t k 1 ) v k ~ N ( 0 , σ 1 2 Δ t k ) , and ε k ~ N ( 0 , σ 1 ε 2 ) . Whereas if  t k > τ , we have s = 2 v k = σ 2 B ( t k τ ) B ( t k 1 τ ) v k ~ N ( 0 , σ 2 2 Δ t k ) , and ε k ~ N ( 0 , σ 2 ε 2 ) . In addition, it is assumed that v k k 1 and ε s k k 1 are independent of each other.
The KF algorithm is adopted to solve Equation (27). For x k ~ N ( x ^ k | k , P k | k ) , if t k τ , we define x ^ k | k = E x k | Y 1 : k , Θ and P k | k = var x k | Y 1 : k , Θ , if t k > τ , we define x ^ k | k = E x k | Y τ ˜ + 1 : k , Θ , P k | k = var x k | Y τ ˜ + 1 : k , Θ . Then, the KF process could be summarized as follows,
State estimation:
x ^ k | k 1 = x ^ k 1 | k 1 + λ s p , k 1 Δ T k x ^ k | k = x ^ k | k 1 + K ( k ) ( y k x ^ k | k 1 ) K ( k ) = P k | k 1 ( P k | k 1 + σ s ε 2 ) 1 P k | k 1 = P k 1 | k 1 + σ s 2 ( t k t k 1 )
Variance update:
P k | k = 1 K ( k ) P k | k 1
where s = 1 , 2 denotes the phase, when s = 1 , λ 1 p , k 1 could be obtained based on Equation (25), when s = 2 , λ 2 p , k 1 could be obtained based on Equation (26). In addition, the initial values of x ^ k | k , P k | k are set as x ^ 0 | 0 = 0 , P 0 | 0 = 0 based on x 0 = 0 in Equation (1).
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.

5. 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.

5.1. 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., 0 t μ 1 ( ρ ; ϑ 1 ) d ρ = t b 1 , τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ = ( t τ ) b 2 . 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: λ 1 p = 1.2 ,   σ 1 p = 0.2 ,   b 1 = 1.5 ,   σ 1 = 0.05 ,   σ 1 ε = 0.3 ,   λ 2 p = 1.5 ,   σ 2 p = 0.25 ,   b 2 = 1.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 σ 1 p and σ 2 p 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,
M S E k = 0 + ( l ˜ k l k ) 2 f L k | Y 1 : k ( l k | Y 1 : k ) d l k  
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
A E k = l ˜ k l ^ k  
where l ^ k 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 Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 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,
T M S E = k = 1 m M S E k
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,
M A E = 1 m k = 1 m A E k
where A E 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,
C R A = 1 m k = 1 m 1 A E k l ˜ k
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.
Overall, the numerical simulation results demonstrate that the proposed prediction RUL method considering three-source variability and two-phase nonlinearity is effective, and could improve the estimation accuracy.

5.2. 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 λ 1 p , λ 2 p 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,
R E k = l ˜ k l ^ k l ˜  
where l ˜ k and l ^ k 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.
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.

6. 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 multi-component systems (such as hybrid energy systems consisting of capacitors and batteries) with three-source variability is also a valuable research direction in the future.

Author Contributions

Conceptualization, X.C.; methodology, X.C.; validation, X.C., Y.H. and J.L.; formal analysis, X.C.; investigation, J.L. and X.C.; resources, Y.H. and J.L.; data curation, X.C.; writing—original draft preparation, X.C.; writing—review and editing, Y.H. and J.L.; visualization, X.C.; supervision, Y.H. and J.L.; project administration, Y.H. and J.L. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported in part by the Key R&D Programs of the Ministry of Science and Technology of China under Grant No. 2020YFB1712602.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are openly available in Ref. [39].

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

The following abbreviations and notations are used in this manuscript:
Abbreviations
PHMPrognostics and health management
RULRemaining useful life
PDFProbability density function
CDFCumulative distribution function
CMCondition monitoring
KFKalman filtering
EMExpectation maximization
ELMExtreme learning machine
FPTFirst passage time
MSEMean square error
AEAbsolute error
HIhealth indicator
Notation
X ( t ) Actual degradation state at time t
x 0 Initial value
τ Time of the changing point
x τ Degradation state at the changing time
λ 1 Drift coefficient of the first phase
λ 2 Drift coefficient of the second phase
σ 1 Diffusion coefficients of the first phase
σ 2 Diffusion coefficients of the second phase
μ 1 ( ρ ; ϑ 1 ) Nonlinear function of the first phase
μ 2 ( ρ τ ; ϑ 2 ) Nonlinear function of the second phase
ϑ 1 Nonlinear   parameter   in   μ 1 ( ρ ; ϑ 1 )
ϑ 2 Nonlinear   parameter   in   μ 2 ( ρ τ ; ϑ 2 )
B ( t ) Standard Brown motion
λ 1 p , σ 1 p Mean   and   standard   deviation   of   λ 1
λ 2 p , σ 2 p Mean   and   standard   deviation   of   λ 2
Y ( t ) Measurement process
ε 1 , ε 2 Measurement errors of each phase
σ 1 ε , σ 2 ε Standard   deviation   of   ε 1   and   ε 2
w Failure threshold
T Lifetime
L k RUL   at   the   current   time   t k
Y 1 : k Observations   up   to   t k
f T ( ) PDF of the lifetime
f L ( ) PDF of the RUL
Φ CDF of the standard normal distribution
ϕ PDF of the standard normal distribution
h τ ( x τ ) Transition probability density function
x k Actual   degradation   state   at   the   current   time   t k
x ^ k | k , P k | k Mean   and   variance   of   x k
Θ Unknown model parameter vector
N Number of the tested devices from the same batch
Y 1 : N Historical data of N devices from the same batch
X 1 : N Actual degradation state data of N devices from the same batch
Y n Observations   of   the  n-th device
Δ y n , j Measured   increment   at   time   t j of   the  n-th device
Δ Y n Measured increment vector
m n Available   number   of   the   observations   for   the  n-th device
τ n Changing   time   of   the  n-th device
τ ˜ n Changing   point   location   of   the  n-th device
μ τ , σ τ Mean and standard deviation of τ
Y 1 : k Degradation   observations   at   the   current   time   t k
X 1 : k actual   degradation   states   at   the   current   time   t k

Appendix A. Proof of Theorem 1

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., λ 1 ~ N ( λ 1 p , σ 1 p 2 ) , λ 2 ~ N ( λ 2 p , σ 2 p 2 ) , 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,
f T ( t ) w 0 t μ 1 ( ρ ; ϑ 1 ) d ρ t μ 1 ( t ; ϑ 1 ) × σ 1 p 2 w 0 t μ 1 ( ρ ; ϑ 1 ) d ρ + λ 1 p σ 1 2 t 0 t μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 t 2 π t 2 0 t μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 t exp w λ 1 p 0 t μ 1 ( ρ ; ϑ 1 ) d ρ 2 2 0 t μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 t , 0 < t τ Q R , t > τ  
where Q = Q 1 Q 2 , R = R 1 R 2 , and
Q 1 = r a 1 2 2 π ( t τ ) 2 ( σ α 1 2 + σ β 1 2 ) exp ( λ α 1 λ β 1 ) 2 2 ( σ α 1 2 + σ β 1 2 ) × λ β 1 σ α 1 2 + λ α 1 σ β 1 2 σ α 1 2 + σ β 1 2 × Φ λ β 1 σ α 1 2 + λ α 1 σ β 1 2 σ α 1 2 σ β 1 2 ( σ α 1 2 + σ β 1 2 ) + σ α 1 2 σ β 1 2 σ α 1 2 + σ β 1 2 × ϕ λ β 1 σ α 1 2 + λ α 1 σ β 1 2 σ α 1 2 σ β 1 2 ( σ α 1 2 + σ β 1 2 ) Q 2 = r b 1 2 2 π ( t τ ) 2 ( σ α 1 2 + σ β 1 2 ) exp ( λ α 1 λ β 1 ) 2 2 ( σ α 1 2 + σ β 1 2 ) × 1 Φ λ β 1 σ α 1 2 + λ α 1 σ β 1 2 σ α 1 2 σ β 1 2 ( σ α 1 2 + σ β 1 2 ) R 1 = I 1 × r a 1 2 2 π ( t τ ) 2 ( σ α 1 2 + σ β 1 2 ) exp ( λ α 1 λ γ 1 ) 2 2 ( σ α 1 2 + σ β 1 2 ) × λ γ 1 σ α 1 2 + λ α 1 σ β 1 2 σ α 1 2 + σ β 1 2 × Φ λ γ 1 σ α 1 2 + λ α 1 σ β 1 2 σ α 1 2 σ β 1 2 ( σ α 1 2 + σ β 1 2 ) + σ α 1 2 σ β 1 2 σ α 1 2 + σ β 1 2 × ϕ λ γ 1 σ α 1 2 + λ α 1 σ β 1 2 σ α 1 2 σ β 1 2 ( σ α 1 2 + σ β 1 2 ) R 2 = I 1 × r b 1 2 2 π ( t τ ) 2 ( σ α 1 2 + σ β 1 2 ) exp ( λ α 1 λ γ 1 ) 2 2 ( σ α 1 2 + σ β 1 2 ) × 1 Φ λ γ 1 σ α 1 2 + λ α 1 σ β 1 2 σ α 1 2 σ β 1 2 ( σ α 1 2 + σ β 1 2 ) λ α 1 = λ 2 p τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ , λ β 1 = w λ 1 p 0 τ μ 1 ( ρ ; ϑ 1 ) d ρ , λ γ 1 = w λ 1 p 0 τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 w σ 1 p 2 0 τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 , σ α 1 2 = σ 2 p 2 τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 + σ 2 2 ( t τ ) , σ β 1 2 = σ 1 p 2 0 τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 + σ 1 2 τ r a 1 = ( t τ ) σ 2 2 + σ 2 p 2 μ 2 ( t τ ; ϑ 2 ) τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ σ α 1 2 , r b 1 = ( t τ ) σ 2 2 λ α 1 λ 2 p ( t τ ) μ 2 ( t τ ; ϑ 2 ) σ α 1 2 I 1 = exp 2 λ 1 p w σ 1 2 + 2 ( w 2 σ 1 p 4 0 τ μ 1 ( ρ ; ϑ 1 ) d ρ + w 2 σ 1 p 2 σ 1 2 ) ( σ 1 2 + σ 1 p 2 0 τ μ 1 ( ρ ; ϑ 1 ) d ρ ) σ 1 4
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 < τ )
f L ( l k ) 1 2 π l k 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k × w x k t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ l k μ 1 ( t k + l k ; ϑ 1 ) × ( w x k ) σ 1 p 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ + λ 1 p σ 1 2 l k t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k × exp w x k λ 1 p t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k , 0 < l k + t k τ S T , τ < l k + t k
where S = S 1 S 2 , T = T 1 T 2 , and
S 1 = r a 2 2 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) exp ( λ α 2 λ β 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) × λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 + σ β 2 2 × Φ λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) + σ α 2 2 σ β 2 2 σ α 2 2 + σ β 2 2 × ϕ λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) S 2 = r b 2 2 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) exp ( λ α 2 λ β 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) × 1 Φ λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) T 1 = I 2 × r a 2 2 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) exp ( λ α 2 λ γ 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) × λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 + σ β 2 2 × Φ λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) + σ α 2 2 σ β 2 2 σ α 2 2 + σ β 2 2 × ϕ λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) T 2 = I 2 × r b 2 2 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) exp ( λ α 2 λ γ 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) × 1 Φ λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) λ α 2 = λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ , λ β 2 = w x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ , λ γ 2 = w + x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 ( w x k ) σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 , σ α 2 2 = σ 2 p 2 τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 + σ 2 2 ( t k + l k τ ) , σ β 2 2 = σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 + σ 1 2 ( τ t k ) , r a 2 = ( t k + l k τ ) σ 2 2 + σ 2 p 2 μ 2 ( t k + l k τ ; ϑ 2 ) τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ α 2 2 , r b 2 = ( t k + l k τ ) σ 2 2 λ a 2 λ 2 p ( t k + l k τ ) μ 2 ( t k + l k τ ; ϑ 2 ) σ α 2 2 , I 2 = exp 2 λ 1 p ( w x k ) σ 1 2 + 2 ( w x k ) 2 σ 1 p 4 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ + ( w x k ) 2 σ 1 p 2 σ 1 2 σ 1 2 + σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 4
Case 2: The current time t k is larger than the changing time τ (i.e., t k τ )
f L ( l k ) 1 2 π l k 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 l k × w x k t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ l k μ 2 ( t k + l k τ ; ϑ 2 ) × σ 2 p 2 ( w x k ) t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ + λ 2 p σ 2 2 l k σ 2 p 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 + σ 2 2 l k × exp w x k λ 2 p t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 l k
The result of Equation (A3) is transformed from lifetime T to RUL based on the law of total probability as f T ( t ) = + f T ( t | λ 2 , x τ ) p ( λ 2 ) d λ 2 , 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. 
 
Case 1: (1) If t k τ and 0 < l k + t k τ ,
To facilitate calculation, the following Lemma A1 is provided.
Lemma A1  
([19]). If  Z ~ N ( μ , σ 2 )  and  w 1 , w 2 , A , B R + , then
E Z ( w 1 A Z ) exp ( w 2 B Z ) 2 2 C = C B 2 σ 2 + C w 1 A B w 2 σ 2 + μ C B 2 σ 2 + C exp ( w 2 B μ ) 2 2 ( B 2 σ 2 + C )
(A) Based on f T ( t | λ 1 ) in Equation (5) and λ 1 ~ N ( λ 1 p , σ 1 p 2 ) , the PDF of the lifetime considering the unit-to-unit variability could be formulated via Lemma A1 as,
f T ( t ) = + f T ( t | λ 1 ) p ( λ 1 ) d λ 1 + w x 0 λ 1 0 t μ 1 ( ρ ; ϑ 1 ) d ρ t μ 1 ( t ; ϑ 1 ) 2 π σ 1 2 t 3 exp w x 0 λ 1 0 t μ 1 ( ρ ; ϑ 1 ) d ρ 2 2 σ 1 2 t p ( λ 1 ) d λ 1 w 0 t μ 1 ( ρ ; ϑ 1 ) d ρ t μ 1 ( t ; ϑ 1 ) × w σ 1 p 2 0 t μ 1 ( ρ ; ϑ 1 ) d ρ + λ 1 p σ 1 2 t 0 t μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 t 2 π t 2 0 t μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 t × exp w λ 1 p 0 t μ 1 ( ρ ; ϑ 1 ) d ρ 2 2 0 t μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 t
(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,
f L ( l k ) w x k t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ l k μ 1 ( t k + l k ; ϑ 1 ) × ( w x k ) σ 1 p 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ + λ 1 p σ 1 2 l k t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k 2 π l k 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k × exp w x k λ 1 p t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k
(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,
f L k | Y 1 : k ( l k | Y 1 : k ) E x k | Y 1 : k w x k t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ l k μ 1 ( t k + l k ; ϑ 1 ) × ( w x k ) σ 1 p 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ + λ 1 p σ 1 2 l k t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k 2 π l k 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k × exp w x k λ 1 p t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k
(D) Then, based on x k ~ N ( x ^ k | 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,
f L k | Y 1 : k ( l k | Y 1 : k ) M 1 N 1 × w λ 1 p t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ P k | k + υ 1 ( l k ) x ^ k | k P k | k + υ 1 ( l k ) 2 π l k 2 υ 1 2 ( l k ) P k | k + υ 1 ( l k ) exp w x ^ k | k λ 1 p t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 2 P k | k + υ 1 ( l k )
where,
υ 1 ( l k ) = t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ 2 σ 1 p 2 + σ 1 2 l k , M 1 = w υ 1 ( l k ) w σ 1 p 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ + λ 1 p σ 1 2 l k t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ l k μ 1 ( t k + l k ; ϑ 1 ) , N 1 = υ 1 ( l k ) σ 1 p 2 t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ t k t k + l k μ 1 ( ρ ; ϑ 1 ) d ρ l k μ 1 ( t k + l k ; ϑ 1 )
(2) If t k τ and τ < l k + t k ,
(A) 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,
S = S 1 S 2 = 1 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) exp ( λ α 2 λ β 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) × r a 2 × λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 + σ β 2 2 r b 2 × Φ λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) + r a 2 σ α 2 2 σ β 2 2 σ α 2 2 + σ β 2 2 × ϕ λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) T = T 1 T 2 = I 2 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) exp ( λ α 2 λ γ 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) × r a 2 × λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 + σ β 2 2 r b 2 × Φ λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) + r a 2 σ α 2 2 σ β 2 2 σ α 2 2 + σ β 2 2 × ϕ λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 )
where
λ α 2 = λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ , λ β 2 = w x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ , λ γ 2 = w + x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 ( w x k ) σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 , σ α 2 2 = σ 2 p 2 τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 + σ 2 2 ( t k + l k τ ) , σ β 2 2 = σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 + σ 1 2 ( τ t k ) , r a 2 = ( t k + l k τ ) σ 2 2 + σ 2 p 2 μ 2 ( t k + l k τ ; ϑ 2 ) τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ α 2 2 , r b 2 = ( t k + l k τ ) σ 2 2 λ a 2 λ 2 p ( t k + l k τ ) μ 2 ( t k + l k τ ; ϑ 2 ) σ α 2 2 , I 2 = exp 2 λ 1 p ( w x k ) σ 1 2 + 2 ( w x k ) 2 σ 1 p 4 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ + ( w x k ) 2 σ 1 p 2 σ 1 2 σ 1 2 + σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 4
(B) Then, based on Equation (A9) and x k ~ N ( x ^ k | k , P k | k ) , the PDF of RUL with three-source variability can be formulated as follows,
f L k | Y 1 : k ( l k | Y 1 : k ) = E x k | Y 1 : k S T = E x k | Y 1 : k S E x k | Y 1 : k T
(C) Equation (A10) can be obtained through solving E x k | Y 1 : k S and E x k | Y 1 : k T separately. To facilitate calculation, the following Lemmas A2 and A3 are provided.
Lemma A2 
([45]). If  Z ~ N ( μ , σ 2 ) w , A R  and   B R + , 
E Z exp ( w A Z ) 2 2 B = B A 2 σ 2 + B exp ( w A μ ) 2 2 ( A 2 σ 2 + B )
Lemma A3 
([19]). If  Z ~ N ( μ , σ 2 ) ,  w , A , B , C , D R  and  1 - 2 B σ 2 > 0 , then,
E Z exp ( AZ + BZ 2 ) Φ ( C + DZ ) = 1 1 2 B σ 2 exp 2 A μ + A 2 σ 2 + 2 B μ 2 2 1 2 B σ 2 Φ C + D μ + A D σ 2 2 B C σ 2 1 - 2 B σ 2 2 + D 2 σ 2 1 - 2 B σ 2
(D) Then, to calculate E x k | Y 1 : k S , the S formula in Equation (A9) is expanded to S a and S b as follows,
S a = exp ( λ α 2 λ β 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) × r a 2 × λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 + σ β 2 2 r b 2 × Φ λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) = exp λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ w x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 2 ( σ α 2 2 + σ β 2 2 ) × r a 2 × w x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ α 2 2 + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ β 2 2 σ α 2 2 + σ β 2 2 r b 2 × Φ w x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ α 2 2 + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 )
S b = exp ( λ α 2 λ β 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) × ϕ λ β 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) = exp λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ w x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 2 ( σ α 2 2 + σ β 2 2 ) × ϕ w x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ α 2 2 + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 )
where
σ α 2 2 = σ 2 p 2 τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 + σ 2 2 ( t k + l k τ ) , σ β 2 2 = σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 + σ 1 2 ( τ t k ) , r a 2 = ( t k + l k τ ) σ 2 2 + σ 2 p 2 μ 2 ( t k + l k τ ; ϑ 2 ) τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ α 2 2 , r b 2 = ( t k + l k τ ) σ 2 2 λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ λ 2 p ( t k + l k τ ) μ 2 ( t k + l k τ ; ϑ 2 ) σ α 2 2 ,
(E) Based on Equation (A13), E x k | Y 1 : k S could be expressed as follows,
E x k | Y 1 : k S = 1 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) × E x k | Y 1 : k S a + r a 2 σ α 2 2 σ β 2 2 σ α 2 2 + σ β 2 2 E x k | Y 1 : k S b
(F) To solve E x k | Y 1 : k S a , we define the following parameter simplification formulas,
a S = w λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ , b S = σ α 2 2 + σ β 2 2 , c S = r a 2 × σ α 2 2 σ α 2 2 + σ β 2 2 , d S = r a 2 × w λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ α 2 2 + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ β 2 2 σ α 2 2 + σ β 2 2 r b 2 , e S = w λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ α 2 2 + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ β 2 2 , f S = σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 )
Then, the formula S a in Equation (A13) could be simplified as,
S a = exp ( a S x k ) 2 2 b S × d S + c S x k × Φ e S σ α 2 2 x k f S
(G) Thus, based on x k ~ N ( x ^ k | k , P k | k ) and Lemma A3, E x k | Y 1 : k S a could be obtained as follows,
E x k | Y 1 : k S a = E x k | Y 1 : k exp ( a S x k ) 2 2 b S × d S × Φ e S σ α 2 2 x k f S + E x k | Y 1 : k exp ( a S x k ) 2 2 b S × c S x k × Φ e S σ α 2 2 x k f S = d S exp a S 2 2 b S 1 + P k | k b S exp 2 a S x ^ k | k + a S 2 P k | k b 3 x ^ k | k 2 2 b S + P k | k Φ e S σ α 2 2 x ^ k | k a S σ α 2 2 P k | k b S + e S P k | k b S f S 2 1 + P k | k b S 2 + σ α 2 4 P k | k 1 + P k | k b S + 1 2 π P k | k + exp ( a S x k ) 2 2 b S × c S x k × exp ( x k x ^ k | k ) 2 2 P k | k × Φ e S σ α 2 2 x k f S d x k
(H) To solve E x k | Y 1 : k S b , based on Equation (15), the S b in Equation (A13) could be simplified as,
S b = exp ( a S x k ) 2 2 b S × 1 2 π exp ( e S σ α 2 2 x k ) 2 2 f S 2
Thus, E x k | Y 1 : k S b could be formulated as follows,
E x k | Y 1 : k S b = E x k | Y 1 : k exp ( a S x k ) 2 2 b S × 1 2 π exp ( e S σ α 2 2 x k ) 2 2 f S 2 = 1 2 π E x k | Y 1 : k exp 1 2 f S 2 + b S σ α 2 4 b S f S 2 x k 2 2 a S f S 2 + b S e S σ α 2 2 b S f S 2 x k + a S 2 f S 2 + b S e S 2 b S f S 2 = 1 2 π E x k | Y 1 : k exp 1 2 ( L 1 S x k L 2 S ) 2 + L 3 S
where
L 1 S = f S 2 + b S σ α 2 4 b S f S 2 , L 2 S = a S f S 2 + b S e S σ α 2 4 b S f S 2 f S 2 + b S σ α 2 4 , L 3 S = a S 2 f S 2 + b S e S 2 b S f S 2 a S f S 2 + b S e S σ α 2 2 2 b S f S 2 f S 2 + b S σ α 2 4
Then, according to Lemma A2, E x k | Y 1 : k S b could be obtained as,
E x k | Y 1 : k S b = 1 2 π E x k | Y 1 : k exp 1 2 ( L 1 S x k L 2 S ) 2 + L 3 S = 1 2 π exp 1 2 L 3 S × E x k | Y 1 : k exp 1 2 ( L 1 S x k L 2 S ) 2 = 1 2 π exp 1 2 L 3 S × 1 L 1 S 2 P k | k + 1 exp L 2 S L 1 S x ^ k | k 2 2 L 1 S 2 P k | k + 1
(I) Finally, bringing the results of E x k | Y 1 : k S a and E x k | Y 1 : k S b into Equation (A14), the final result of E x k | Y 1 : k S , i.e., the formula S n e w in Theorem 2, could be obtained. Due to space limitations, the result of S n e w is omitted here.
Similarly, we can solve E x k | Y 1 : k T in Equation (A10) as follows.
(A) Firstly, I 2 in the formula T of Equation (A9) should be simplified as,
I 2 = exp 2 λ 1 p ( w x k ) σ 1 2 + 2 ( w x k ) 2 σ 1 p 4 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ + ( w x k ) 2 σ 1 p 2 σ 1 2 σ 1 2 + σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 4 = exp 2 w λ 1 p σ 1 2 + σ 1 p 2 w σ 1 4 exp 2 σ 1 p 2 σ 1 4 x k 2 2 λ 1 p σ 1 2 + 4 w σ 1 p 2 σ 1 4 x k = exp 2 w λ 1 p σ 1 2 + σ 1 p 2 w σ 1 4 exp a T x k 2 b T x k
where a T = 2 σ 1 p 2 σ 1 4 , b T = 2 λ 1 p σ 1 2 + 4 w σ 1 p 2 σ 1 4 .
(B) Further, as to
exp ( λ α 2 λ γ 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) = exp λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ + w x k + λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ + 2 ( w x k ) σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 2 2 ( σ α 2 2 + σ β 2 2 )
we define
c T = w + λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ + 2 w σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 d T = 1 + 2 σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 = σ 1 2 + 2 σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2
Thus, the Equation (A22) could be simplified as,
exp c T d T x k 2 2 b S
Based on Equations (A21) and (A24), we have
I 2 exp ( λ α 2 λ γ 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) = exp 2 w λ 1 p σ 1 2 + σ 1 p 2 w σ 1 4 exp a T x k 2 b T x k × exp c T d T x k 2 2 b S = exp 2 w λ 1 p σ 1 2 + σ 1 p 2 w σ 1 4 exp d T 2 2 a T b S 2 b S x k 2 + c T d T b S b T b S x k c T 2 2 b S = exp 2 w λ 1 p σ 1 2 + σ 1 p 2 w σ 1 4 exp c T 2 2 b S exp D 1 T 2 b S x k 2 + D 2 T b S x k
where D 1 T = d T 2 2 a T b S , D 2 T = 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 ) exp ( λ α 2 λ γ 2 ) 2 2 ( σ α 2 2 + σ β 2 2 ) = exp 2 w λ 1 p σ 1 2 + σ 1 p 2 w σ 1 4 exp c T 2 2 b S exp D 1 T 2 b S x k 2 + D 2 T b S x k 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) = exp 4 w b S λ 1 p σ 1 2 + σ 1 p 2 w c T 2 σ 1 4 2 b S σ 1 4 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) × exp D 1 T 2 b S x k 2 + D 2 T b S x k
(D) On this basis, the formula T in Equation (A9) could be split into T a and T b as follows,
T a = exp D 1 T 2 b S x k 2 + D 2 T b S x k × r a 2 × λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 + σ β 2 2 r b 2 × Φ λ γ 2 σ α 2 2 + λ α 2 σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 ) = exp D 1 T 2 b S x k 2 + D 2 T b S x k × r a 2 × w + x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 ( w x k ) σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 σ α 2 2 + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ β 2 2 σ α 2 2 + σ β 2 2 r b 2 × Φ w + x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 ( w x k ) σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 σ α 2 2 + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 )
T b = exp D 1 T 2 b S x k 2 + D 2 T b S x k × ϕ w + x k λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 ( w x k ) σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 σ α 2 2 + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ β 2 2 σ α 2 2 σ β 2 2 ( σ α 2 2 + σ β 2 2 )
(E) Thus, the simplified form of E x k | Y 1 : k T could be obtained based on T a , T b ,
E x k | Y 1 : k T = exp 4 w b S λ 1 p σ 1 2 + σ 1 p 2 w c T 2 σ 1 4 2 b S σ 1 4 2 π ( t k + l k τ ) 2 ( σ α 2 2 + σ β 2 2 ) × E x k | Y 1 : k T a + r a 2 × σ α 2 2 σ β 2 2 σ α 2 2 + σ β 2 2 E x k | Y 1 : k T b
(F) To solve E x k | Y 1 : k T a , we define
e T = r a 2 σ α 2 2 σ 1 2 + 2 σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 σ α 2 2 + σ β 2 2 , f T = r a 2 × w λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 w σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 σ α 2 2 + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ β 2 2 σ α 2 2 + σ β 2 2 r b 2 , g T = w λ 1 p t k τ μ 1 ( ρ ; ϑ 1 ) d ρ 2 w σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2 σ α 2 2 + λ 2 p τ t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ σ β 2 2 , h T = σ α 2 2 σ 1 2 + 2 σ 1 p 2 t k τ μ 1 ( ρ ; ϑ 1 ) d ρ σ 1 2
Then, the E x k | Y 1 : k T a in Equation (A29) could be formulated as follows,
E x k | Y 1 : k T a = E x k | Y 1 : k exp D 1 T 2 b S x k 2 + D 2 T b S x k × f T + e T x k × Φ g T + h T x k f S = E x k | Y 1 : k exp D 1 T 2 b S x k 2 + D 2 T b S x k × f T × Φ g T + h T x k f S + E x k | Y 1 : k exp D 1 T 2 b S x k 2 + D 2 T b S x k × e T x k × Φ g T + h T x k f S
(G) Based on x k ~ N ( x ^ k | 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 = f T 1 + D 1 T P k | k b S exp 2 D 2 T x ^ k | k + D 2 T 2 P k | k b 3 D 1 T x ^ k | k 2 2 b S + D 1 T P k | k Φ g T + h T x ^ k | k + D 2 T h T P k | k b S + D 1 T g T P k | k b S f S 2 1 + D 1 T P k | k b S 2 + h T 2 P k | k 1 + D 1 T P k | k b S + 1 2 π P k | k + exp D 1 T 2 b S x k 2 + D 2 T b S x k × e T x k × exp ( x k x ^ k | k ) 2 2 P k | k × Φ g T + h T x k f S d x k
(H) Similarly, according to Equation (A30), E x k | Y 1 : k T b could be formulated as,
E x k | Y 1 : k T b = E x k | Y 1 : k exp D 1 T 2 b S x k 2 + D 2 T b S x k × 1 2 π exp ( g T + h T x k ) 2 2 f S 2 = 1 2 π E x k | Y 1 : k exp 1 2 D 1 T f S 2 + b S h T 2 b S f T 2 x k 2 2 D 2 T f S 2 b S g T h T b S f S 2 x k + g T 2 f T 2 = 1 2 π E x k | Y 1 : k exp 1 2 R 1 T x k R 2 T 2 R 3 T = 1 2 π exp 1 2 R 3 T E x k | Y 1 : k exp 1 2 R 1 T x k R 2 T 2
where
R 1 T = D 1 T f S 2 + b S h T 2 b S f S 2 , R 2 T = D 2 T f S 2 b S g T h T b S f S 2 D 1 T f S 2 + b S h T 2 , R 3 T = g T 2 f S 2 D 2 T f S 2 b S g T h T 2 b S f S 2 D 1 T f S 2 + b S h T 2
(I) Based on x k ~ N ( x ^ k | k , P k | k ) and Lemma A1, E x k | Y 1 : k T b could be obtained as follows,
E x k | Y 1 : k T b = 1 2 π exp 1 2 R 3 T × 1 R 1 T 2 P k | k + 1 exp R 2 T R 1 T x ^ k | k 2 2 R 1 T 2 P k | k + 1
(J) Finally, bringing the final results of E x k | Y 1 : k T a and E x k | Y 1 : k T b into Equation (A29), the final results of E x k | Y 1 : k T , i.e., the formula T n e w in Theorem 2, could be obtained. Due to space limitations, the result of T n e w is omitted here.
Case 2: If τ < t k ,
(A) When τ < t k , the degradation process is equivalent to a single-phase nonlinear process, and the randomness of the degradation state at the changing point can be omitted. Therefore, based on λ 2 ~ N ( λ 2 p , σ 2 p 2 ) and the f T ( t | λ 2 , x τ ) in Equation (6), the PDF of the lifetime T could be obtained through Lemma A1 as follows,
f T ( t ) = + f T ( t | λ 2 , x τ ) p ( λ 2 ) d λ 2 + w x τ λ 2 τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ ( t τ ) μ 2 ( t τ ; ϑ 2 ) 2 π σ 2 2 ( t τ ) 3 exp w x τ λ 2 τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 2 σ 2 2 ( t τ ) p ( λ 2 ) d λ 2 1 2 π ( t τ ) 2 τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 ( t τ ) × w x τ τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ ( t τ ) μ 2 ( t τ ; ϑ 2 ) × ( w x τ ) σ 2 p 2 τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ + λ 2 p σ 2 2 ( t τ ) τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 ( t τ ) × exp w x τ λ 2 p τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 2 τ t μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 ( t τ )
(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,
f L ( l k ) w x k t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ l k μ 2 ( t k + l k τ ; ϑ 2 ) × ( w x k ) σ 2 p 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ + λ 2 p σ 2 2 l k t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 l k 2 π l k 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 l k × exp w x k λ 2 p t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 l k
(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 is presented as,
f L k | Y τ ˜ : k ( l k | Y τ ˜ : k ) = E x k | Y τ ˜ : k 1 2 π l k 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 l k × w x k t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ l k μ 2 ( t k + l k τ ; ϑ 2 ) × ( w x k ) σ 2 p 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ + λ 2 p σ 2 2 l k t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 l k × exp w x k λ 2 p t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 l k
where τ ˜ denote the index of the changing point location.
(D) Then, based on x k ~ N ( x ^ k | k , P k | k ) and Lemma A1, we could obtain the PDF of RUL with three-source variability simultaneously when τ < t k as follows,
f L k | Y τ ˜ : k ( l k | Y τ ˜ : k ) = M 2 N 2 × w λ 2 p t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ P k | k + υ 2 ( l k ) x ^ k | k P k | k + υ 2 ( l k ) 2 π l k 2 υ 2 2 ( l k ) P k | k + υ 2 ( l k ) × exp w x ^ k | k λ 2 p t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 2 P k | k + υ 2 ( l k )
where
υ 2 ( l k ) = t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ 2 σ 2 p 2 + σ 2 2 l k , M 2 = w υ 2 ( l k ) w σ 2 p 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ + λ 2 p σ 2 2 l k t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ l k μ 2 ( t k + l k τ ; ϑ 2 ) N 2 = υ 2 ( l k ) σ 2 p 2 t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ t k t k + l k μ 2 ( ρ τ ; ϑ 2 ) d ρ l k μ 2 ( t k + l k τ ; ϑ 2 )
In this way, the proof has been completed. □

References

  1. Li, R.; Verhagen, W.J.C.; Curran, R. A systematic methodology for Prognostic and Health Management system architecture definition. Reliab. Eng. Syst. Safe 2020, 193, 106598. [Google Scholar] [CrossRef]
  2. Zhang, Y.; Fang, L.; Qi, Z.; Deng, H. A Review of Remaining Useful Life Prediction Approaches for Mechanical Equipment. IEEE Sens. J. 2023, 23, 29991–30006. [Google Scholar] [CrossRef]
  3. Vogl, G.W.; Weiss, B.A.; Helu, M. A review of diagnostic and prognostic capabilities and best practices for manufacturing. J. Intell. Manuf. 2019, 30, 79–95. [Google Scholar] [CrossRef]
  4. Hu, Y.; Miao, X.; Si, Y.; Pan, E.; Zio, E. Prognostics and health management: A review from the perspectives of design, development and decision. Reliab. Eng. Syst. Safe 2022, 217, 108063. [Google Scholar] [CrossRef]
  5. Kim, S.; Choi, J.; Kim, N.H. Challenges and Opportunities of System-Level Prognostics. Sensors 2021, 21, 7655. [Google Scholar] [CrossRef]
  6. Si, X.; Wang, W.; Hu, C.; Zhou, D. Remaining useful life estimation—A review on the statistical data driven approaches. Eur. J. Oper. Res. 2011, 213, 1–14. [Google Scholar] [CrossRef]
  7. Lui, Y.H.; Li, M.; Downey, A.; Shen, S.; Nemani, V.P.; Ye, H.; VanElzen, C.; Jain, G.; Hu, S.; Laflamme, S.; et al. Physics-based prognostics of implantable-grade lithium-ion battery for remaining useful life prediction. J. Power Sources 2021, 485, 229327. [Google Scholar] [CrossRef]
  8. Guo, J.; Li, Z.; Li, M. A Review on Prognostics Methods for Engineering Systems. IEEE Trans. Reliab. 2020, 69, 1110–1129. [Google Scholar] [CrossRef]
  9. Lei, Y.; Li, N.; Guo, L.; Li, N.; Yan, T.; Lin, J. Machinery health prognostics: A systematic review from data acquisition to RUL prediction. Mech. Syst. Signal Process. 2018, 104, 799–834. [Google Scholar] [CrossRef]
  10. Sayyad, S.; Kumar, S.; Bongale, A.; Kotecha, K.; Abraham, A. Remaining Useful-Life Prediction of the Milling Cutting Tool Using Time–Frequency-Based Features and Deep Learning Models. Sensors 2023, 23, 5659. [Google Scholar] [CrossRef]
  11. Lei, Y.; Li, N.; Lin, J. A New Method Based on Stochastic Process Models for Machine Remaining Useful Life Prediction. IEEE Trans. Instrum. Meas. 2016, 65, 2671–2684. [Google Scholar] [CrossRef]
  12. Zhang, J.X.; Du, D.B.; Si, X.S.; Liu, Y.; Hu, C.H. Prognostics Based on Stochastic Degradation Process: The Last Exit Time Perspective. IEEE Trans. Reliab. 2021, 70, 1158–1176. [Google Scholar] [CrossRef]
  13. Zhang, S.; Zhai, Q.; Shi, X.; Liu, X. A Wiener Process Model With Dynamic Covariate for Degradation Modeling and Remaining Useful Life Prediction. IEEE Trans. Reliab. 2023, 72, 214–223. [Google Scholar] [CrossRef]
  14. Zhang, Z.; Si, X.; Hu, C.; Lei, Y. Degradation data analysis and remaining useful life estimation: A review on Wiener-process-based methods. Eur. J. Oper. Res. 2018, 271, 775–796. [Google Scholar] [CrossRef]
  15. Cai, B.; Fan, H.; Shao, X.; Liu, Y.; Liu, G.; Liu, Z.; Ji, R. Remaining useful life re-prediction methodology based on Wiener process: Subsea Christmas tree system as a case study. Comput. Ind. Eng. 2021, 151, 106983. [Google Scholar] [CrossRef]
  16. Si, X.S.; Wang, W.; Hu, C.H.; Zhou, D.H.; Pecht, M.G. Remaining Useful Life Estimation Based on a Nonlinear Diffusion Degradation Process. IEEE Trans. Reliab. 2012, 61, 50–67. [Google Scholar] [CrossRef]
  17. Zhang, Z.X.; Si, X.S.; Hu, C.H.; Zhang, Q.; Li, T.M.; Xu, C.Q. Planning Repeated Degradation Testing for Products with Three-Source Variability. IEEE Trans. Reliab. 2016, 65, 640–647. [Google Scholar] [CrossRef]
  18. Wang, Y.; Liu, Q.; Lu, W.; Peng, Y. A general time-varying Wiener process for degradation modeling and RUL estimation under three-source variability. Reliab. Eng. Syst. Safe 2023, 232, 109041. [Google Scholar] [CrossRef]
  19. Si, X.S.; Wang, W.; Hu, C.H.; Zhou, D.H. Estimating Remaining Useful Life With Three-Source Variability in Degradation Modeling. IEEE Trans. Reliab. 2014, 63, 167–190. [Google Scholar] [CrossRef]
  20. Feng, L.; Wang, H.; Si, X.; Zou, H. A State-Space-Based Prognostic Model for Hidden and Age-Dependent Nonlinear Degradation Process. IEEE Trans. Autom. Sci. Eng. 2013, 10, 1072–1086. [Google Scholar] [CrossRef]
  21. Ye, Z.S.; Wang, Y.; Tsui, K.L.; Pecht, M. Degradation Data Analysis Using Wiener Processes With Measurement Errors. IEEE Trans. Reliab. 2013, 62, 772–780. [Google Scholar] [CrossRef]
  22. Si, X.S. An Adaptive Prognostic Approach via Nonlinear Degradation Modeling: Application to Battery Data. IEEE Trans. Ind. Electron. 2015, 62, 5082–5096. [Google Scholar] [CrossRef]
  23. Wang, D.; Tsui, K. Brownian motion with adaptive drift for remaining useful life prediction: Revisited. Mech. Syst. Signal Process. 2018, 99, 691–701. [Google Scholar] [CrossRef]
  24. Cai, Z.; Wang, Z.; Chen, Y.; Guo, J.; Xiang, H. Remaining useful lifetime prediction for equipment based on nonlinear implicit degradation modeling. J. Syst. Eng. Electron. 2020, 31, 194–205. [Google Scholar] [CrossRef]
  25. Zheng, J.F.; Si, X.S.; Hu, C.H.; Zhang, Z.X.; Jiang, W. A Nonlinear Prognostic Model for Degrading Systems with Three-Source Variability. IEEE Trans. Reliab. 2016, 65, 736–750. [Google Scholar] [CrossRef]
  26. Wang, X.; Hu, C.; Si, X.; Zheng, J.; Pei, H.; Yu, Y.; Ren, Z. An Adaptive Prognostic Approach for Newly Developed System With Three-Source Variability. IEEE Access 2019, 7, 53091–53102. [Google Scholar] [CrossRef]
  27. Yu, W.; Tu, W.; Kim, I.Y.; Mechefske, C. A nonlinear-drift-driven Wiener process model for remaining useful life estimation considering three sources of variability. Reliab. Eng. Syst. Safe 2021, 212, 107631. [Google Scholar] [CrossRef]
  28. Burgess, W.L. Valve Regulated Lead Acid battery float service life estimation using a Kalman filter. J. Power Sources 2009, 191, 16–21. [Google Scholar] [CrossRef]
  29. Yan, W.A.; Song, B.W.; Duan, G.L.; Shi, Y.M. Real-time reliability evaluation of two-phase Wiener degradation process. Commun. Stat. Theory Methods 2017, 46, 176–188. [Google Scholar] [CrossRef]
  30. Wang, P.; Tang, Y.; Joo Bae, S.; He, Y. Bayesian analysis of two-phase degradation data based on change-point Wiener process. Reliab. Eng. Syst. Safe 2018, 170, 244–256. [Google Scholar] [CrossRef]
  31. Zhang, J.X.; Hu, C.H.; He, X.; Si, X.S.; Liu, Y.; Zhou, D.H. A Novel Lifetime Estimation Method for Two-Phase Degrading Systems. IEEE Trans. Reliab. 2019, 68, 689–709. [Google Scholar] [CrossRef]
  32. Chen, X.; Liu, Z.; Wang, J.; Yang, C.; Long, B.; Zhou, X. An Adaptive Prediction Model for the Remaining Life of an Li-Ion Battery Based on the Fusion of the Two-Phase Wiener Process and an Extreme Learning Machine. Electronics 2021, 10, 540. [Google Scholar] [CrossRef]
  33. Lin, J.; Liao, G.; Chen, M.; Yin, H. Two-phase degradation modeling and remaining useful life prediction using nonlinear wiener process. Comput. Ind. Eng. 2021, 160, 107533. [Google Scholar] [CrossRef]
  34. Hu, C.; Xing, Y.; Du, D.; Si, X.; Zhang, J. Remaining useful life estimation for two-phase nonlinear degradation processes. Reliab. Eng. Syst. Safe 2023, 230, 108945. [Google Scholar] [CrossRef]
  35. Wang, P.; Tang, Y.; Bae, S.J.; Xu, A. Bayesian Approach for Two-Phase Degradation Data Based on Change-Point Wiener Process With Measurement Errors. IEEE Trans. Reliab. 2018, 67, 688–700. [Google Scholar] [CrossRef]
  36. Guan, Q.; Wei, X.; Bai, W.; Jia, L. Two-stage degradation modeling for remaining useful life prediction based on the Wiener process with measurement errors. Qual. Reliab. Eng. Int. 2022, 38, 3485–3512. [Google Scholar] [CrossRef]
  37. Lin, W.; Chai, Y. Remaining useful life prediction for nonlinear two-phase degradation process with measurement errors and imperfect prior information. Meas. Sci. Technol. 2023, 34, 055018. [Google Scholar] [CrossRef]
  38. Lin, J.; Lin, Z.; Liao, G.; Yin, H. A Novel Product Remaining Useful Life Prediction Approach Considering Fault Effects. IEEE/CAA J. Autom. Sin. 2021, 8, 1762–1773. [Google Scholar] [CrossRef]
  39. Feng, J.; Sun, Q.; Jin, T. Storage Life Prediction for a High-Performance Capacitor Using Multi-Phase Wiener Degradation Model. Commun. Stat. Simul. Comput. 2012, 41, 1317–1335. [Google Scholar] [CrossRef]
  40. Xu, X.; Tang, S.; Yu, C.; Xie, J.; Han, X.; Ouyang, M. Remaining Useful Life Prediction of Lithium-ion Batteries Based on Wiener Process Under Time-Varying Temperature Condition. Reliab. Eng. Syst. Safe 2021, 214, 107675. [Google Scholar] [CrossRef]
  41. Yan, J.; Liu, Y.; Li, L.; Ren, X. Wind Turbine Condition Monitoring Using the SSA-Optimized Self-Attention BiLSTM Network and Changepoint Detection Algorithm. Sensors 2023, 23, 5873. [Google Scholar] [CrossRef] [PubMed]
  42. Ali, N.; Liu, Z.; Armghan, H.; Armghan, A. Double integral sliding mode controller for wirelessly charging of fuel cell-battery-super capacitor based hybrid electric vehicle. J. Energy Storage 2022, 51, 104288. [Google Scholar] [CrossRef]
  43. Lemian, D.; Bode, F. Battery-Supercapacitor Energy Storage Systems for Electrical Vehicles: A Review. Energies 2022, 15, 5683. [Google Scholar] [CrossRef]
  44. Mohammed, A.S.; Atnaw, S.M.; Salau, A.O.; Eneh, J.N. Review of optimal sizing and power management strategies for fuel cell/battery/super capacitor hybrid electric vehicles. Energy Rep. 2023, 9, 2213–2228. [Google Scholar] [CrossRef]
  45. Tang, S.; Guo, X.; Yu, C.; Zhou, Z.; Zhou, Z.; Zhang, B. Real time remaining useful life prediction based on nonlinear Wiener based degradation processes with measurement errors. J. Cent. South Univ. 2014, 21, 4509–4517. [Google Scholar] [CrossRef]
Figure 1. The simulated degradation paths.
Figure 1. The simulated degradation paths.
Sensors 24 00165 g001
Figure 2. The estimated changing time of 100 degradation paths.
Figure 2. The estimated changing time of 100 degradation paths.
Sensors 24 00165 g002
Figure 3. The online degradation path.
Figure 3. The online degradation path.
Sensors 24 00165 g003
Figure 4. The online updating process. (a) The parameter updating. (b) The underlying degradation state updating.
Figure 4. The online updating process. (a) The parameter updating. (b) The underlying degradation state updating.
Sensors 24 00165 g004
Figure 5. The RUL prediction results. (a) PDFs of the RUL. (b) The mean RUL of the four methods.
Figure 5. The RUL prediction results. (a) PDFs of the RUL. (b) The mean RUL of the four methods.
Sensors 24 00165 g005
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 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.
Sensors 24 00165 g006
Figure 7. Relative capacitance variability of the capacitors with time.
Figure 7. Relative capacitance variability of the capacitors with time.
Sensors 24 00165 g007
Figure 8. Online parameter updating processes of Capacitor 1.
Figure 8. Online parameter updating processes of Capacitor 1.
Sensors 24 00165 g008
Figure 9. PDFs of RUL prediction for Capacitor 1. (a) Our method. (b) Lin’s method. (c) Chai’ method. (d) Zheng’s method.
Figure 9. PDFs of RUL prediction for Capacitor 1. (a) Our method. (b) Lin’s method. (c) Chai’ method. (d) Zheng’s method.
Sensors 24 00165 g009
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 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.
Sensors 24 00165 g010
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.
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.
Sensors 24 00165 g011
Table 1. Parameter estimation values with different sample sizes.
Table 1. Parameter estimation values with different sample sizes.
Size λ 1 p σ 1 p σ 1 σ 1 ε b 1 λ 2 p σ 2 p σ 2 σ 2 ε b 2 μ τ σ τ
n = 10 1.32790.33700.04500.29631.49941.65600.28510.09720.39671.402750.58001.6302
n = 50 1.25590.25000.05180.29831.50021.47940.25570.09730.40001.401649.61802.1199
n = 100 1.22430.23110.05050.29951.50011.48650.24820.10410.40241.401749.74001.9923
True value1.20.20.050.31.51.50.250.080.41.4502
Table 2. Comparison results of RUL prediction based on TMSE, MAE, and CRA.
Table 2. Comparison results of RUL prediction based on TMSE, MAE, and CRA.
MetricTMSEMAECRA
Our method1.2030 × 1040.86340.9653
Lin’s method [33]2.3221 × 1043.74150.8795
Chai’s method [37]3.1273 × 1043.95290.8976
Zheng’s method [25]6.9231 × 1048.37260.8037
Table 3. The parameter estimation results of the capacitance degradation data.
Table 3. The parameter estimation results of the capacitance degradation data.
Variable λ 1 p σ 1 p σ 1 σ 1 ε b 1 λ 2 p σ 2 p σ 2 σ 2 ε b 2 μ τ σ τ
value0.03810.00870.16140.04972.08000.22320.04070.22740.03441.52176.25000.8292
Table 4. Comparison results of RUL prediction based on TMSE, MAE, and CRA for Capacitor 1 degradation data.
Table 4. Comparison results of RUL prediction based on TMSE, MAE, and CRA for Capacitor 1 degradation data.
MetricTMSEMAECRA
Our method4.93490.10910.9685
Lin’s method [33]7.05920.36370.8862
Chai’s method [37]7.59040.61820.8647
Zheng’s method [25]16.25561.04550.7873
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Cui, X.; Lu, J.; Han, Y. Remaining Useful Life Prediction for Two-Phase Nonlinear Degrading Systems with Three-Source Variability. Sensors 2024, 24, 165. https://doi.org/10.3390/s24010165

AMA Style

Cui X, Lu J, Han Y. Remaining Useful Life Prediction for Two-Phase Nonlinear Degrading Systems with Three-Source Variability. Sensors. 2024; 24(1):165. https://doi.org/10.3390/s24010165

Chicago/Turabian Style

Cui, Xuemiao, Jiping Lu, and Yafeng Han. 2024. "Remaining Useful Life Prediction for Two-Phase Nonlinear Degrading Systems with Three-Source Variability" Sensors 24, no. 1: 165. https://doi.org/10.3390/s24010165

APA Style

Cui, X., Lu, J., & Han, Y. (2024). Remaining Useful Life Prediction for Two-Phase Nonlinear Degrading Systems with Three-Source Variability. Sensors, 24(1), 165. https://doi.org/10.3390/s24010165

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop