# Residual Lifetime Prediction with Multistage Stochastic Degradation for Equipment

Authors: Zhan Gao; Qi-guo Hu; Xiang-yang Xu
Journal: Complexity
 (2020)

---

## Abstract

Residual useful lifetime (RUL) prediction plays a key role of failure prediction and health management (PHM) in equipment. Aiming at the problems of residual life prediction without comprehensively considering multistage and individual differences in equipment performance degradation at present, we explore a prediction model that can fit the multistage random performance degradation. Degradation modeling is based on the random Wiener process. Moreover, according to the degradation monitoring data of the same batch of equipment, we apply the expectation maximization (EM) algorithm to estimate the prior distribution of the model. The real-time remaining life distribution of the equipment is acquired by merging prior information of real-time degradation data and historical degradation monitoring data. The accuracy of the proposed model is demonstrated by analyzing a practical case of metalized film capacitors, and the result shows that a better RUL estimation accuracy can be provided by our model compared with existing models.


---

## Body

## 1. Introduction

Prognostics and health management (PMH) has been a systematic method utilized to evaluate the reliability or residual life of actual life-cycle conditions in a system, predict failure degree, reduce the risk in operating, improve task completion rate, and make maintenance decisions. Prognostics is usually characterized by estimating the RUL of a system based on condition monitoring (CM) data [1, 2]. Once a prognosis is available, some health management tactics, such as maintenance, replacement, and logistic support, can be performed properly to obtain the operational objectives, which is required [3, 4].In general, RUL prediction is an important step in PMH. Over the last decade, there are a great deal of research studies as to RUL prediction have been developed for a large wide range of industrial products, for example, bearings [5, 6], gearboxes [7], lithium batteries [8], organic light emitting diodes [9], and laser generators [10]. The main idea of RUL is to realize the life prediction by obtaining the distribution or expectation of RUL according to the effective information such as equipment failure mechanism, CM data, and failure time data [11, 12]. For fields with high safety and reliability requirements, for instance, aerospace and electronic equipment, whose performance gradually decreases as service time increases, their final failure is manifested as degraded failure. At the same time, due to measurement errors, environmental interference, and other factors, the degradation process has a certain degree of randomness. At present, life prediction based on statistical data is a central issue in current research [13, 14]. Based on the probability density function (PDF) of the RUL with the randomness of the degradation process, the PDF of the RUL density can be obtained.In life prediction based on statistical data, stochastic processes are often used to fit the degradation process, such as Gamma process [15, 16], Wiener process [17, 18], and Inverse Gaussian process [19, 20]. The Wiener process is convenient to describe the degradation system that has a random nonmonotonic independent incremental process, which are influenced by external environment, internal state changes, and load conditions [21]. Li et al. in [22] proposed a novel RUL prediction models combined with degradation processes described by establishing the system age process and the system state process. They calculated the PDF of the RUL based on a closed-form solution by calculating the degradation trajectories values. Si et al. in [2] derived a random Wiener process model to predict the RUL of equipment by merging a recursive filtering algorithm and an EM algorithm for parameter estimation. With due consideration of the influence of measurement error, a nonlinear random Wiener model is proposed to predict remaining useful life [23].The above research has enriched the application of stochastic processes in life prediction; however, a large number of the current degradation models made assumptions that the degradation process is subject to a single degradation law. Due to the staged changes in internal failure mechanism and dynamic changes in external environmental condition, many degradation process of equipment is characterized by multistage degradation [24]. However, for the current research in the Wiener process of the performance degradation modeling, there are few of them who target on multistage performance degradation models. Therefore, Wang et al. in [25] further proposed a two-stage degradation model of both Gamma process and Wiener process. In reference to the attenuation process of the brightness of light emitting diodes, Wang et al. in [9] assumed that individual differences exist in the change points and utilized a two-stage Wiener process model and Bayesian estimation method to analyze the degradation data under logarithmic transformation. In addition, how to integrate the degradation data of similar equipment to estimate the hyperparameters of the degradation model is the key to the RUL prediction. A two-stage maximum likelihood estimation (TSMLE) method was utilized to calculate the hyperparameters of the degradation model to achieve life prediction in some literatures [26]. Zhang et al. in [27] proposed a real-time residual life prediction algorithm including the Cumulative Sum (CUSUM) method and the Brownian motion with drift. However, the CUSUM method is only suitable for change-point degradation analysis and lacks the use of prior knowledge; when degradation data is small, the estimation accuracy of the two-step method is not high. Additionally, the change points in the multistage stochastic degradation process are sometimes scattered, so the RUL prediction based on the random process of the change point analysis is often not accurate enough.Aiming at the above problems, in this paper, a multistage random degradation model based on the Wiener process is proposed to estimate the RUL of equipment aiming at above issues. Assuming that the parameters of the degradation model follow, respectively, certain random distributions to describe the differences between single equipment, an EM algorithm is used for iteratively estimating the prior distribution of the degradation model parameters based on the historical degradation data and historical failure time data. Finally, we apply the Bayesian algorithm to update the posterior distribution parameters of our model at each stage. By analyzing the degradation process of metallized film capacitors, the remaining life expectancy of individual capacitors can be predicted.


## 2. Residual Life Prediction Model of Equipment

### 2.1. Model Analysis and Assumptions

Generally, a degradation model of the Wiener process can be described as(1)Xt=ut+σBt,where u is the drift coefficient, σ is the diffusion coefficient, and Bt is the standard Brownian motion.When we predict the remaining life of the equipment, the following assumptions should be made:(1)
If the performance degradation of the equipment exceeds the failure threshold, degradation failure will occur.(2)
The performance degradation of equipment has the characteristics of multiple stages, each performance degradation stage obeys the Wiener process, and it is divided by a certain fixed degradation amount.(3)
The drift coefficientμ and diffusion coefficient σ of the Wiener process are random variables. We assume that ω=1/σ2 and f=u,ω as joint prior distribution of u and ω, which obeys the Gaussian-Gamma distribution. f (u,ω) can be expressed as(2)fμ,ω=ω1/22πcexp−ωμ−c22d⋅baΓaωa−1exp−bω,where a, b, c, and d represent the hyperparameters in prior distribution of u and ω.


### 2.2. Performance Degradation Modeling

In fact, the degradation process of the equipment is characterized by multistage and randomness. A multistage degradation model of equipment can be expressed as(3)Xt=X0+X1tI0,t1t+Xt1+X2t−t1It1,t2t+⋯+Xtn+Xnt−tnItn,∞t,where It is the indicative function, X0 is the initial value of the Wiener degradation model, and ti represents the time that the performance degradation of each stage reaching its boundary value. Xkk=1,2,…,n is the performance degradation process of each stage.The degradation process of the equipment in each stage can be described by using the Wiener process without loss of generality, so the performance degradation amountXt of the equipment can be expressed as(4)Xt=X0+μ1t+σ1BtI0,t1t+D1+μ2t−t1+σ2Bt−t1It1,t2t+⋯+Dn+μnt−tn+σnBt−tnItn,∞t,where Dk is the boundary value of kth performance degradation stage, Dn represents the degradation failure threshold, and μk and ωk are the drift coefficient and diffusion coefficient of the kth stage of Wiener degradation, respectively.The remaining lifeLτ of the equipment represents the time from the performance degradation value Xτ at current time τ to the failure threshold Dn According to the definition of performance failure, the remaining life is expressed as(5)Lτ=inft:Xt+τ≥Dn,t>0.DenoteT as the life of the equipment, according to the relationship between the life and the remaining life, and it can be expressed as(6)T=inft:Xt≥Dn,t>0=Lτ+τ.Denoteξk as the time when the device reaches the threshold value Dk of the kth stage of performance degradation from the initial moment, and ξk can be expressed as(7)ξk=inft:Xt≥Dk,t>0,=ξk−1+inft:X′t≥Dk−Dk−1,t>0,where X′t=μkt+σkBt. Because the failure time of the Wiener process obeys the Inverse Gaussian distribution, ξk obeys the Inverse Gaussian distribution fΔξΔξ|μk,ωk, and it is formulated as(8)fΔξΔξ|μk,ωk=ωk1/2ΔDk2πΔξk3expωkΔDk−μkΔξk22Δξk.Therefore, the probability density distribution function in thekth stage degradation process of the equipment can be formulated as(9)fΔξk|μk,ωk=∫0∞fΔξξk|ξk−1dξk=∫0∞fΔξξk−ξk−1dξk.When the performance degradation valueXτ at τ of the equipment meets the condition Xξk−1≤Xτ≤Xξk, the remaining life of the degradation stage of equipment can be expressed as(10)Lτ=ξk−τ+Δξk+1+Δξk+2+⋯+Δξn.The probability density distribution function of the remaining lifeLτ of equipment degradation is expressed as(11)fLτLτ|μk,ωk,μk+1,ωk+1,…,μn,ωn=∫τ∞∫τ∞⋯∫τ∞fΔξlτ+τ−ξn−1fΔξξn−1−ξn−2,…,fΔξξk−τdξkdξk+1,…,dξn−1.Because ofEΔξk=ΔDk/μk, the expectation of remaining life Lτ can be expressed as(12)ELτ=Wk−Xτμk+ΔDk+1μk+1+⋯ΔDnμn.


## 2.1. Model Analysis and Assumptions

Generally, a degradation model of the Wiener process can be described as(1)Xt=ut+σBt,where u is the drift coefficient, σ is the diffusion coefficient, and Bt is the standard Brownian motion.When we predict the remaining life of the equipment, the following assumptions should be made:(1)
If the performance degradation of the equipment exceeds the failure threshold, degradation failure will occur.(2)
The performance degradation of equipment has the characteristics of multiple stages, each performance degradation stage obeys the Wiener process, and it is divided by a certain fixed degradation amount.(3)
The drift coefficientμ and diffusion coefficient σ of the Wiener process are random variables. We assume that ω=1/σ2 and f=u,ω as joint prior distribution of u and ω, which obeys the Gaussian-Gamma distribution. f (u,ω) can be expressed as(2)fμ,ω=ω1/22πcexp−ωμ−c22d⋅baΓaωa−1exp−bω,where a, b, c, and d represent the hyperparameters in prior distribution of u and ω.


## 2.2. Performance Degradation Modeling

In fact, the degradation process of the equipment is characterized by multistage and randomness. A multistage degradation model of equipment can be expressed as(3)Xt=X0+X1tI0,t1t+Xt1+X2t−t1It1,t2t+⋯+Xtn+Xnt−tnItn,∞t,where It is the indicative function, X0 is the initial value of the Wiener degradation model, and ti represents the time that the performance degradation of each stage reaching its boundary value. Xkk=1,2,…,n is the performance degradation process of each stage.The degradation process of the equipment in each stage can be described by using the Wiener process without loss of generality, so the performance degradation amountXt of the equipment can be expressed as(4)Xt=X0+μ1t+σ1BtI0,t1t+D1+μ2t−t1+σ2Bt−t1It1,t2t+⋯+Dn+μnt−tn+σnBt−tnItn,∞t,where Dk is the boundary value of kth performance degradation stage, Dn represents the degradation failure threshold, and μk and ωk are the drift coefficient and diffusion coefficient of the kth stage of Wiener degradation, respectively.The remaining lifeLτ of the equipment represents the time from the performance degradation value Xτ at current time τ to the failure threshold Dn According to the definition of performance failure, the remaining life is expressed as(5)Lτ=inft:Xt+τ≥Dn,t>0.DenoteT as the life of the equipment, according to the relationship between the life and the remaining life, and it can be expressed as(6)T=inft:Xt≥Dn,t>0=Lτ+τ.Denoteξk as the time when the device reaches the threshold value Dk of the kth stage of performance degradation from the initial moment, and ξk can be expressed as(7)ξk=inft:Xt≥Dk,t>0,=ξk−1+inft:X′t≥Dk−Dk−1,t>0,where X′t=μkt+σkBt. Because the failure time of the Wiener process obeys the Inverse Gaussian distribution, ξk obeys the Inverse Gaussian distribution fΔξΔξ|μk,ωk, and it is formulated as(8)fΔξΔξ|μk,ωk=ωk1/2ΔDk2πΔξk3expωkΔDk−μkΔξk22Δξk.Therefore, the probability density distribution function in thekth stage degradation process of the equipment can be formulated as(9)fΔξk|μk,ωk=∫0∞fΔξξk|ξk−1dξk=∫0∞fΔξξk−ξk−1dξk.When the performance degradation valueXτ at τ of the equipment meets the condition Xξk−1≤Xτ≤Xξk, the remaining life of the degradation stage of equipment can be expressed as(10)Lτ=ξk−τ+Δξk+1+Δξk+2+⋯+Δξn.The probability density distribution function of the remaining lifeLτ of equipment degradation is expressed as(11)fLτLτ|μk,ωk,μk+1,ωk+1,…,μn,ωn=∫τ∞∫τ∞⋯∫τ∞fΔξlτ+τ−ξn−1fΔξξn−1−ξn−2,…,fΔξξk−τdξkdξk+1,…,dξn−1.Because ofEΔξk=ΔDk/μk, the expectation of remaining life Lτ can be expressed as(12)ELτ=Wk−Xτμk+ΔDk+1μk+1+⋯ΔDnμn.


## 3. Bayesian Deduction of Hyperparameters

When the boundary valueDk of each performance degradation stage is determined, μ and σ of the Wiener process jointly determining the probability density distribution of the remaining life, the parameters μ and ω need to be solved. Because the degradation process of individual equipment is different, in order to describe it conveniently, we assume that μ and ω obey the following conjugate prior distribution:(13)ω=σ−2∼Gammaa,b,μ|ω∼Nc,dω.When the real-time degradation dataX=X1,X2,…Xn is obtained, we can utilize the Bayesian formula to derive the posterior distribution πμ,ω|X, which is formulated as(14)πμ,ω|X=LX|μ,ω⋅πμ,ω∫0∞∫−∞∞LX|μ,ω⋅πμ,ωdμdω,where LX|μ,ω presents log-likelihood function for parameters μ and ω, which is expressed as(15)LX|μ,ω=∏i=1nω1/22πΔΛtiexp−ωΔXi−μΔΛti22ΔΛti.πμ,ω is the joint prior density function of μ and ω, which is formulated as(16)πμ,ω=ω1/22πdexp−ωμ−c22dbaωa−1Γaexp−bω.In view of the derivation above, we combine equations (14)–(16) to obtain the following form:(17)πμ,ω|X∝LX|μ,ω⋅πμ,ω∝ωn+1/2+a−1⋅exp−ω2μ2Λtn−2μXn+∑i=1nΔXi2ΔΛti−ω2μ−c2d−bω∝,ωn+1/2+a−1exp−ωb+c22d−Xnd+c22Λtnd2+d+∑i=1nΔXi22ΔΛti⋅ω1/2expω2μ−Xnd+c/Λtnd+12d/Λtnd+1.Because the conjugate prior distribution of random parameters and its posterior distribution have the same distribution function, the difference is only the change of the hyperparameter value in different degradation stages. According to equation (17), we can obtain the posterior estimation values ak′,bk′,ck′,anddk′ of the hyperparameters in the kth performance degradation stage deduced as(18)ak′=nk2+ak,bk′=bk+ck22dk−Xn,kdk+ck22Λtn,kdk2+dk+∑i=1nΔXi,k22ΔΛti,k,ck′=Xn,kdk+ckΛtn,kdk+1,dk′=dkΛtn,kdk+1.The posterior meanμ′ and ω′ of μ and ω can be obtained from a′,b′,c′,andd′, and from equation (17), the kth performance degradation stage is expressed as(19)πμk,ωk|ΔXk=bk′ak′ωkak′−1Γak′exp−ωkbk′⋅ωk′2πdk′exp−ωkμk−ck′22πdk′.After getting the result ofπμk,ωk|ΔXk, we can obtain the posterior mean values μk′andωk′ in kth stage to realize the prediction of the remaining life of the equipment.


## 4. Parameter Estimation of Degradation Model

### 4.1. Determination of the Complete Maximum Likelihood Function of the Hyperparameters

When the boundary value of each stage of performance degradation is determined, the probability density distribution of the remaining lifeLτ can be determined by the Wiener process parameters μ and ω in each performance degradation stage. Therefore, parameters μ and ω is required to estimate. Because we assumed that the unknown parameters μ and ω are random variables, the unknown hyperparameters in the prior distribution of parameters can be estimated based on the historical performance degradation data and historical failure time data [28].For the same kind of equipment, we can collect the degradation data ofN equipments as priori information in CM time. Assuming that the performance degradation data of the ith (i = 1,2, …, N) equipment at the jth (j = 1,2, …, n) measurement is denoted as Xij, the corresponding measurement time is tij. The kΔxkand Δtk of the historical data sample of the kth performance degradation stage are, respectively, expressed as(20)Δxk=Δx11,kΔx12,k⋯Δx1n,kΔx21,kΔx22,k⋯Δx2n,k⋮⋮⋮ΔxN1,kΔxN2,k⋯ΔxNn,k,Δtk=Δt11,kΔt12,k⋯Δt1n,kΔt21,kΔt22,k⋯Δt2n,k⋮⋮⋮ΔtN1,kΔtN2,k⋯ΔtNn,k.The probability density function of the degradation incrementΔXkt, in the kth stage of degradation based on the characteristics of the Wiener process, can be expressed as(21)fΔXk=12πΔtkσk2exp−ΔXk−μkΔtk2σk2Δtk.Therefore, the complete likelihood functions of hyperparametersak,bk,ck,anddkof the historical degradation data tij,k,Xij,k in the kth stage are expressed as(22)LΔXk|ak,bk,ck,dk=∏i=1N∏j=1nhΔXij,k⋅∏i=1Nπμi,k,ωi,k=∏i=1N∏j=1nωi,k1/22πΔtij,kexp−ωi,kΔXij,k−μi,kΔtij,k22Δtij,k⋅∏i=1Nbkakωi,kak−1Γakexp−ωi,kbk⋅ωi,k1/22πdkexp−ωi,kμi,k−ck22dk.In theory, we can obtain the priori estimates of hyperparameters for thekth stage Wiener process a^k, b^k, c^k, and d^k through the traditional maximum likelihood estimation (MLE), and it can be expressed as(23)ψa^k−lna^k=∑i=1Nlnωi,kN+lnN−ln∑i−1Nωi,k,b^k=N⋅a^k∑i=1Nωi,k,c^k=∑i=1Nμi,kωi,k,d^k=∑i=1Nωi,kμi,k2−2c^kμi,kωi,k+c^k2ωi,kN,where ψ⋅ is the digamma distribution function, μi,kandωi,k are implicit data, and the unknown items of the implicit data must be output first, otherwise, we cannot obtain the results of a^k, b^k, c^k, and d^kdirectly.


### 4.2. Estimation of Hyperparameter

The EM algorithm is an effective method for maximum likelihood estimation [29]. It is mainly used for parameter estimation of missing data, for instance, the EM algorithm can be applied to estimate the parameters of likelihood function when observed data is incomplete or the expression of the likelihood function is so complicated that the traditional estimation method of the maximum likelihood function could be failure.Each recursive iteration process of the EM algorithm containsE steps and M steps. The first step of the EM algorithm (E-step): given the observed data ΔXk and current hyperparameter ak,bk,ck,anddk estimated value, we can get the conditional expected value of unknown latent variable Z=μ,ω in complete data log-likelihood function LΔXk|ak,bk,ck,dk. The second step of the EM algorithm (M-step): solve the ak,bk,ck,anddk and make them as parameter estimate values for the next iteration.InE-step, if iteration process performs r times in the kth stage of the degradation process, al,k,bl,k,cl,k, and dl,k are the estimated value in lth (1 ≤ l ≤ r). The expectation that each unknown item contains implicit data μi,k,ωi,k in the l + 1 step iteration process is(24)Eωi,k=a^l,k+Ni/2b^l,k+c^l,k2/2d^l,k−d^l,kXiNi,k+c^l,k2/2d^ltiNi,kd^l,k+1+∑j=1NiΔXij,k2/2ΔΛtij,k,Elnωi,k=ψa^l,k+Ni/2−lnb^l,k+c^l,k2/2dl.k−dl,kXiNi,k+c^l,k2/2dl,kdl,ktiNi,k+1+∑j=1NiΔXij,k2/2tij,k,Eωi,k,μi,k=Eωi,k⋅d^l,kXiNi.k+c^l,ktiNi,kd^l,k+1,Eωi,k,μi,k2=Eωi,k⋅d^l,kXiNi,k+c^l,ktiNi,kd^l,k+12+d^l,ktiNi,kd^l,k+1.InM-step, we substitute the expectations of μi,k and ωi,k that are obtained in E-step into equation (23), and the estimated values of hyperparameters a^l+1, b^l+1, c^l+1, and d^l+1after l + 1 iterations under the kth stage degradation can be obtained as(25)ψa^l+1,k−lna^l+1,k=∑i=1NElnωi,kN+lnN−ln∑i=1NEωi,k,b^l+1,k=N⋅a^l+1,k∑i=1NEωi,k,c^l+1,k=∑i=1NEωi,kμi,k∑i=1NEωi,k,d^l+1,k=∑i=1NEωi,kμi,k2−2c^l+1,kEωi,kμi,k+c^l+1,k2Eωi,kN.If the error betweena^r−1, b^r−1, c^r−1, and d^r−1 and a^r, b^r, c^r, and d^r reaches the set condition after r iterations, then the a priori estimate of the hyperparameter is obtained.


### 4.3. Real-Time Remaining Life Distribution of Equipment

If joint posterior distribution of the parametersμ and ω are obtained, the marginal probability distributions can be obtained, respectively. Because the joint posterior distribution of μ and ω obeys Gaussian-gamma distribution, the posterior distribution of μ obeys μ∼t2a′d′ and the posterior distribution of ω obeys 1/ω′∼IGa′,1/b′. IG (·) is the Inverse Gamma distribution, so the posterior distribution of μ and ω can be expressed as(26)μk=Eμk|ΔXk=dk′,ωk=Eωk|ΔXk=ak′bk′.Moreover, the real-time remaining life density distribution of the equipment can be further expressed as(27)fLτLτ=fLτLτ|μ^1,ω^1,μ^2,ω^2,…,μ^k,ω^k.


### 4.4. Prediction Steps for Remaining Life of Single Equipment

When the real-time performance degradation data of a single equipment is obtained, the remaining life of a single aero engine can be predicted in real time, and the specific steps are as follows:Step 1
: we estimate the unknown hyperparameters in the prior distribution of the parametersμ and ω of the degradation model at each stage, getting estimated values a^, b^, c^, and d^by the EM algorithm by historical degradation data X and failure time data ξ.Step 2
: when the real-time performance degradation data ΔX of the equipment is obtained, the model parameters in each performance degradation stage are updated in real time using equation (18) according to the Bayesian method. Moreover, we get the estimated hyperparameter values a^′, b^′, c^′, and d^′ in posterior distribution in our model.Step 3
: according to equation (26) the posterior estimated values μ^ and ω^ of our model in each performance degradation stage are calculated, respectively, and substituting them into equation (27), the remaining life density distribution of equipment after Bayesian updating can be obtained. Finally, we can estimate the expectancy value of equipment remaining life though calculating equation (12).


## 4.1. Determination of the Complete Maximum Likelihood Function of the Hyperparameters

When the boundary value of each stage of performance degradation is determined, the probability density distribution of the remaining lifeLτ can be determined by the Wiener process parameters μ and ω in each performance degradation stage. Therefore, parameters μ and ω is required to estimate. Because we assumed that the unknown parameters μ and ω are random variables, the unknown hyperparameters in the prior distribution of parameters can be estimated based on the historical performance degradation data and historical failure time data [28].For the same kind of equipment, we can collect the degradation data ofN equipments as priori information in CM time. Assuming that the performance degradation data of the ith (i = 1,2, …, N) equipment at the jth (j = 1,2, …, n) measurement is denoted as Xij, the corresponding measurement time is tij. The kΔxkand Δtk of the historical data sample of the kth performance degradation stage are, respectively, expressed as(20)Δxk=Δx11,kΔx12,k⋯Δx1n,kΔx21,kΔx22,k⋯Δx2n,k⋮⋮⋮ΔxN1,kΔxN2,k⋯ΔxNn,k,Δtk=Δt11,kΔt12,k⋯Δt1n,kΔt21,kΔt22,k⋯Δt2n,k⋮⋮⋮ΔtN1,kΔtN2,k⋯ΔtNn,k.The probability density function of the degradation incrementΔXkt, in the kth stage of degradation based on the characteristics of the Wiener process, can be expressed as(21)fΔXk=12πΔtkσk2exp−ΔXk−μkΔtk2σk2Δtk.Therefore, the complete likelihood functions of hyperparametersak,bk,ck,anddkof the historical degradation data tij,k,Xij,k in the kth stage are expressed as(22)LΔXk|ak,bk,ck,dk=∏i=1N∏j=1nhΔXij,k⋅∏i=1Nπμi,k,ωi,k=∏i=1N∏j=1nωi,k1/22πΔtij,kexp−ωi,kΔXij,k−μi,kΔtij,k22Δtij,k⋅∏i=1Nbkakωi,kak−1Γakexp−ωi,kbk⋅ωi,k1/22πdkexp−ωi,kμi,k−ck22dk.In theory, we can obtain the priori estimates of hyperparameters for thekth stage Wiener process a^k, b^k, c^k, and d^k through the traditional maximum likelihood estimation (MLE), and it can be expressed as(23)ψa^k−lna^k=∑i=1Nlnωi,kN+lnN−ln∑i−1Nωi,k,b^k=N⋅a^k∑i=1Nωi,k,c^k=∑i=1Nμi,kωi,k,d^k=∑i=1Nωi,kμi,k2−2c^kμi,kωi,k+c^k2ωi,kN,where ψ⋅ is the digamma distribution function, μi,kandωi,k are implicit data, and the unknown items of the implicit data must be output first, otherwise, we cannot obtain the results of a^k, b^k, c^k, and d^kdirectly.


## 4.2. Estimation of Hyperparameter

The EM algorithm is an effective method for maximum likelihood estimation [29]. It is mainly used for parameter estimation of missing data, for instance, the EM algorithm can be applied to estimate the parameters of likelihood function when observed data is incomplete or the expression of the likelihood function is so complicated that the traditional estimation method of the maximum likelihood function could be failure.Each recursive iteration process of the EM algorithm containsE steps and M steps. The first step of the EM algorithm (E-step): given the observed data ΔXk and current hyperparameter ak,bk,ck,anddk estimated value, we can get the conditional expected value of unknown latent variable Z=μ,ω in complete data log-likelihood function LΔXk|ak,bk,ck,dk. The second step of the EM algorithm (M-step): solve the ak,bk,ck,anddk and make them as parameter estimate values for the next iteration.InE-step, if iteration process performs r times in the kth stage of the degradation process, al,k,bl,k,cl,k, and dl,k are the estimated value in lth (1 ≤ l ≤ r). The expectation that each unknown item contains implicit data μi,k,ωi,k in the l + 1 step iteration process is(24)Eωi,k=a^l,k+Ni/2b^l,k+c^l,k2/2d^l,k−d^l,kXiNi,k+c^l,k2/2d^ltiNi,kd^l,k+1+∑j=1NiΔXij,k2/2ΔΛtij,k,Elnωi,k=ψa^l,k+Ni/2−lnb^l,k+c^l,k2/2dl.k−dl,kXiNi,k+c^l,k2/2dl,kdl,ktiNi,k+1+∑j=1NiΔXij,k2/2tij,k,Eωi,k,μi,k=Eωi,k⋅d^l,kXiNi.k+c^l,ktiNi,kd^l,k+1,Eωi,k,μi,k2=Eωi,k⋅d^l,kXiNi,k+c^l,ktiNi,kd^l,k+12+d^l,ktiNi,kd^l,k+1.InM-step, we substitute the expectations of μi,k and ωi,k that are obtained in E-step into equation (23), and the estimated values of hyperparameters a^l+1, b^l+1, c^l+1, and d^l+1after l + 1 iterations under the kth stage degradation can be obtained as(25)ψa^l+1,k−lna^l+1,k=∑i=1NElnωi,kN+lnN−ln∑i=1NEωi,k,b^l+1,k=N⋅a^l+1,k∑i=1NEωi,k,c^l+1,k=∑i=1NEωi,kμi,k∑i=1NEωi,k,d^l+1,k=∑i=1NEωi,kμi,k2−2c^l+1,kEωi,kμi,k+c^l+1,k2Eωi,kN.If the error betweena^r−1, b^r−1, c^r−1, and d^r−1 and a^r, b^r, c^r, and d^r reaches the set condition after r iterations, then the a priori estimate of the hyperparameter is obtained.


## 4.3. Real-Time Remaining Life Distribution of Equipment

If joint posterior distribution of the parametersμ and ω are obtained, the marginal probability distributions can be obtained, respectively. Because the joint posterior distribution of μ and ω obeys Gaussian-gamma distribution, the posterior distribution of μ obeys μ∼t2a′d′ and the posterior distribution of ω obeys 1/ω′∼IGa′,1/b′. IG (·) is the Inverse Gamma distribution, so the posterior distribution of μ and ω can be expressed as(26)μk=Eμk|ΔXk=dk′,ωk=Eωk|ΔXk=ak′bk′.Moreover, the real-time remaining life density distribution of the equipment can be further expressed as(27)fLτLτ=fLτLτ|μ^1,ω^1,μ^2,ω^2,…,μ^k,ω^k.


## 4.4. Prediction Steps for Remaining Life of Single Equipment

When the real-time performance degradation data of a single equipment is obtained, the remaining life of a single aero engine can be predicted in real time, and the specific steps are as follows:Step 1
: we estimate the unknown hyperparameters in the prior distribution of the parametersμ and ω of the degradation model at each stage, getting estimated values a^, b^, c^, and d^by the EM algorithm by historical degradation data X and failure time data ξ.Step 2
: when the real-time performance degradation data ΔX of the equipment is obtained, the model parameters in each performance degradation stage are updated in real time using equation (18) according to the Bayesian method. Moreover, we get the estimated hyperparameter values a^′, b^′, c^′, and d^′ in posterior distribution in our model.Step 3
: according to equation (26) the posterior estimated values μ^ and ω^ of our model in each performance degradation stage are calculated, respectively, and substituting them into equation (27), the remaining life density distribution of equipment after Bayesian updating can be obtained. Finally, we can estimate the expectancy value of equipment remaining life though calculating equation (12).


## 5. A Practical Case Study

Metalized film capacitor, as one of an important component in Inertial Confinement Fusion (ICF) laser equipment, is the key of life prediction to maintain the laser equipment and to develop a spare parts strategy [30]. However, some problems are unavoidable such as impurity or defect areas in organic film of a metallized film capacitor when it is in the production process. It is prone to form a discharge channel broken down by the external voltage during multiple charge and discharge processes, and because of weak electrical resistance in areas, capacitor result in “self-healing” phenomenon. A single self-healing causes a very small loss of capacitance; however, as the cycle of charge and discharge process increases, the capacitance value of the capacitor decreases continuously. When the loss of capacitance exceeds a certain threshold, the degradation rate of capacitance increases. Usually, we define 5% of the capacity loss as the end-of-life index [31]. The variation of the capacitance degradation under normal working conditions of the eight capacitors is shown in Figure 1, and t is the number of charge and discharge.Figure 1
Capacitance degradation curves.It can be observed from Figure1 that the capacitance of most capacitors shows a decreasing trend at the increase of the number of charge and discharge. The degradation rate of metallized film capacitors is relatively slow in the initial stage. When the capacitance loss reaches about 0.23 μF, the capacitance degradation speed accelerates. We assume that the performance degradation threshold D1 = 0.23. Degradation of another capacitor during working of charge and discharge is shown in Table 1, and we record this capacitor as No. 9 capacitor. Degradation failure occurs when the capacitance change ΔC of the capacitor exceeds 5% of the initial capacitance value C0. We record the failure threshold Dn = 5% C0.Table 1
Capacitance of No. 9 capacitor.
t/shots0100020003000C(t)/μF56.7156.6856.6456.56t/shots4000500060007000C(t)/μF56.3856.1056.0655.82Take degradation data in first five of the 8 capacitors as historical monitoring data for priori parameter estimation, and take the degradation data in last three of them as real-time data for posterior parameter updating accordingly. We assume that the initial values of the hyperparameter is (1, 1, 1, 1), and when the error between the adjacent estimated values after iteration is not greater than 10−6 to end EM iteration process, the priori estimated hyperparameter values a^, b^, c^, and d^ obtained. Moreover, the iterative convergence process of hyperparameter estimation is shown in Figure 2, and the a priori estimated values of hyperparameters at each stage are shown in Table 2.Figure 2
Iterative convergence process of hyperparameter estimation.
(a)(b)(c)(d)Table 2
Estimation of hyperparameter in prior distribution.
StageEstimation value of posteriorab (10−5)c (10−5)d (10−4)11.77101.13653.16180.761220.91532.06415.72051.6162After determining the prior distribution of the model parameters in each performance degradation stage, the real-time performance degradation data of a single capacitor is obtained. Moreover, according to the Bayesian method, the model parameters in each performance degradation stage can be updated in real time. Take No. 9 capacitor as an example for life prediction, and we rely on real-time performance degradation data of last three of capacitors to estimate the hyperparameters of the posterior distribution in the first stage and second stage. Tables3 and 4 show the estimated values after Bayesian updating in different stages. Finally, we substitute the estimated values of Bayesian updating in the first stage and second stage into equation (27), and the probability density distribution of the remaining life can be obtained. Figures 3 and 4 show the real-time remaining life probability density distribution diagrams of No. 9 capacitor in the first and second stages, respectively.Table 3
Updated hyperparameters based on the Bayesian method (1st stage).
Estimation value of posteriort/shots1000200030004000a′1 (10–5)1.77102.27102.77100.9071b′1 (10–5)1.30201.00070.97852.0572c′1 (10–5)3.16183.11283.08625.6825d′1 (10–4)0.76120.72740.66061.6197RL (104/shot)2.8322.7242.6032.470Table 4
Updated hyperparameters based on the Bayesian method (2nd stage).
Estimation value of posteriort/shots500060007000a′20.93531.41531.9153b′22.06412.05872.0529c′25.71785.64215.5129d′21.61621.60981.6103RL (104/shot)2.4022.2242.178Figure 3
Probability density distribution of residual life (first stage).Figure 4
Probability density distribution of residual life (second stage).The EM algorithm has certain advantage in estimation of hyperparameter. Comparing the EM algorithm with the Two-stage Maximum Likelihood Estimation (TSMLE) [26] for the error of the hyperparameter estimation results, we utilize the sum of relative errors Err of 4 hyperparameter estimations to show the precision of the estimation results, and Err is expressed as(28)Err=a^−aa+b^−bb+c^−cc+d^−dd.The estimation error between the two methods is shown in Table5. The result of Err in the EM algorithm and the TSMLE method can be seen. Obviously, the EM algorithm has certain accuracy advantage in the estimation of hyperparameters.Table 5
Comparison of error for the estimation value of hyperparameters.
StageErrEMTSMLE10.08910.200520.12780.4503Table6 shows the Relative Error (RE) of RUL prediction for No.9 Capacitor in our model, model of single-stage degradation [31] and the model of Brownian motion with drift [32]. It can be seen that the relative error of considering the multistage degradation process is smaller than the other models, which proves that the proposed model can provide better accuracy for predicting the remaining life of metallized film capacitors. In addition, due to environmental influence, defects, or impurities of organic film in the same type capacitors, individual capacitor has different degradation. Therefore, the Wiener process distribution parameter estimation based on Bayesian estimation of prior degradation information is more in line with engineering reality.Table 6
RE of RUL prediction for No. 9 Capacitor.
CM time (t/shot)REModel of single-stage degradationOur modelModel of Brownian motion with drift10000.01200.00970.008620000.01950.01350.014230000.02190.02030.022440000.04540.03120.049750000.02680.01270.026460000.04050.03400.043570000.03040.02640.0312


## 6. Conclusion Remarks

In this paper, we propose a degradation model based on the Wiener process combined with the EM algorithm and Bayesian method to estimate the RUL of equipment. A RUL estimation method is constructed with due considering the multistage degradation and individual differences in equipment. By assuming that parameters of the Wiener process are random variables which obey a certain distribution, we utilize EM algorithm to estimate prior parameter of the model based on historical degradation monitoring data and updating posterior parameters based on CM data by Bayesian. Finally, we can obtain different CDF of RUL and expectation value of RUL.By comparing the proposed model with existing models, we utilize RE to measure prediction accuracy, and the result shows that the accuracy of proposed model is better than existing models. The Wiener process that considers individual differences can better describe the degradation process of the equipment, and multistage degradation is more suitable in engineering practical. In addition, we can see that the estimation accuracy of EM is better than the TSMLE method in hyperparameter estimation.



---

*Source: 8847703-2020-11-17.xml*