# Degradation Data-Driven Remaining Useful Life Estimation in the Absence of Prior Degradation Knowledge

Authors: Yong Yu; Changhua Hu; Xiaosheng Si; Jianxun Zhang
Journal: Journal of Control Science and Engineering
 (2017)

---

## Abstract

Recent developments in prognostic and health management have been targeted at utilizing the observed degradation signals to estimate residual life distributions. Current degradation models mainly focus on a population of “identical” devices or an individual device with population information, not a single component in the absence of prior degradation knowledge. However, the fast development of science and technology provides us with many kinds of new systems, and we just have the real-time monitoring information to analyze the reliability for them. The fusion algorithm presented herein addresses this challenge by combining the excellent modeling ability of Bayesian updating method for the multilevel data and the prominent estimation ability of ECM algorithm for incomplete data. Residual life distributions and posterior distributions are first calculated through the Bayesian updating method based on random initial a priori distributions. Then the a priori distributions are revised and improved for future predictions by the ECM algorithm. Once a new signal is observed, we can reuse the fusion algorithm to improve the accuracy of residual life distributions. The applicability of this fusion algorithm is validated by a set of simulation experiments.


---

## Body

## 1. Introduction

Modern engineering systems are overwhelmingly complex because of increasing requirements on their functionalities and qualities. These systems often have a high standard of system reliability because a single failure can lead to catastrophic consequences with profound impacts, extreme costs, and potential safety hazards. It will take an exceedingly long time for a system to fail, so prognostics for systems have become extremely difficult, even if the actual operating conditions are severe and rigorous. Therefore, effective methods that can predict failure progression and evaluate the reliability of the system have long been sought. When there are sufficient monitoring data and efficient computational capability, prognostics for components based on observed degradation data is promising and effective [1–3].The general path model is a typical method utilizing the observed degradation signals of these degraded systems. Lu and Meeker [4] introduced the model to the degradation literature in 1993, for the first time. In their model, the fixed-effects parameters affect the populations’ characteristics and the random-effects parameters describe an individual unit’s characteristics. Once the parameters are known, the residual life is deterministic. Therefore, the core work is to estimate the unknown fixed and random parameters. Lu and Meeker used a two-stage method to estimate the unknown parameters. Lu et al. [5] extended the degradation model and suggested likelihood-based estimation methods. However, these are not suitable for all types of degradation data. Su et al. [6] considered random sample sizes and random repeated measurement times for each product unit. They discussed the advantages and disadvantages of two-stage least-squares (LS) estimation, maximum modified likelihood (MML) estimation, and maximum likelihood estimation (MLE). They showed that the LS estimators are not consistent in the case of random sample sizes. However, MLE can provide a consistent estimator and has smaller biases and variances compared to the LS and MML estimates. In the further study, Weaver et al. [7] also used MLE to estimate the unknown parameters. They extended the research and examined effects of sample size on the estimation precision. Under the mixed-effect path model, Wu and Shao [8] built the asymptotic properties of the (weighted) least-squares estimators. They used these properties to calculate approximate confidence intervals and point estimates for percentiles of the failure-time distribution. They used the weighted least-squares estimators to predict the resistor of metal film and the metal fatigue crack length.However, the above papers only focus on the estimation of the unknown parameters about population devices and need a fair amount of samples in the test. In order to solve this problem, Robinson and Crowder [9] described a fully Bayesian approach which allows a small sample size. They used a variety of simple prior distributions and observed that this aspect has little effect on the posterior distributions of these data, showing that the information in the degradation data dominates. So the Bayesian approach is more suitable for the parameters estimation of an individual device with population information compared with two-stage LS estimation and MLE. Gebraeel et al. [10] developed two different exponential degradation signal models. One model assumes that the error fluctuations follow an iid random error process; however, the second model considers that the error terms follow a Brownian motion process. In their paper, they used the Bayesian updating methods to combine the distribution information of the parameters across the population and the monitored degradation data from the individual device. The Bayesian updating methods can update the stochastic parameters of degradation models, every time a new degradation signal comes. Gebraeel [11] extended the Bayesian updating procedure by assuming that the stochastic parameters in the exponential degradation models follow a bivariate normal distribution. Chakraborty et al. [12] further extended the updating procedure and investigated the difference of the life time distributions when the stochastic parameters do not follow the normal distribution. They also built methods for calculating Remaining Useful Life when the stochastic parameters of the exponential model follow more general distributions. Chen and Tsui [13] adopted a piecewise log-linear degradation model and assumed the change time of the two different phases is random. This assumption explicitly accounts for the characteristics of different degradation phases. When new observations were available, they updated the a posteriori information of the model stochastic parameters including regression coefficients and the variance of the error term by using Bayesian methods. They also suggested a new method which took the correlations into consideration, among degradation predictions, to compute the RUL distribution with better accuracy. Their approach can be naturally extended to more general degradation models.The above Bayesian approaches can be applied to predicting the RUL of an individual device with population information. However, the fast development of science and technology provides us with kinds of newly made systems, and we just have the observed monitoring degradation information to analyze the reliability of them. Traditional ways to predict system failures often fail because the domain knowledge and expert experience are limited and historical data is nonexistent.Therefore, effective methods that can predict failures of these newly made systems in the absence of prior degradation knowledge have long been sought. Considering the excellent modeling ability of Bayesian updating method for the multilevel data and the prominent estimation ability of ECM algorithm for incomplete data, the goal of our paper is to develop a collaboration method between Bayesian updating method and ECM algorithm to estimate the Remaining Useful Life (RUL) of the newly made system just with real-time sensing data. To verify the applicability of this fusion algorithm, a set of simulation experiments are conducted.The remainder of this paper is organized as follows. Section2 develops a separate Bayesian updating method and calculates the RUL distributions for the exponential degradation model with a Brownian error term, under the assumption that the prior distributions are known. Section 3 explains the procedure of the fusion process in detail and estimates the prior information for the Bayesian updating method by our fusion algorithm. Section 4 illustrates the validity of the collaboration method by a set of simulation experiments. The paper concludes with some discussions and guidance for the estimation of RUL of the single component without prior degradation knowledge in Section 5.


## 2. Bayesian Updating Method and Residual Life Distribution

In order to develop our collaboration method, we first introduce the Bayesian updating method and its estimation result of the Remaining Useful Life in this section. We will adopt the exponential degradation signal model developed by Gebraeel et al. [10, 11], Kaiser and Gebraeel [14], and Elwany and Gebraeel [15] and assume that the error fluctuation of the degradation model follows a Brownian motion process. As our objective is to compute the RUL of a single system without prior distribution knowledge, we believe that the adopted exponential degradation model and the error fluctuation are adequate for the given degradation path in Figure 1. For a further discussion of model selection and evaluation, see Li et al. [16].Figure 1
Vibration-based degradation signal.Under the above assumption, we could use the Bayesian updating procedure to compute the unknown random parameters of the exponential degradation model. Once we have got the calculated posterior distributions of these random parameters in the exponential degradation model, we can derive the residual life distribution of the component. However, in our paper we only have the real-time monitoring degradation signal for the newly made component, and prior distributions are nonexistent. In order to utilize the Bayesian updating procedure, we first assume we have got the accurate and informative a priori information. The method of estimating these prior distributions will be illustrated and detailed in Section3.


### 2.1. The Degradation Signal Model

First, we review the general definition of the Brownian motion process.Definition 1.
A standard Brownian motion process,W(t),t≥0, possesses the following properties:1.
Ift0<t1<⋯<tn, then Wt0,W(t1)-W(t0),…,W(tn)-W(tn-1) are mutually independent.2.
Ifs,t≥0, then(1)PWs+t-Ws∈A=∫A2πt-1/2e-x2/2tdx.3.
With probability onet→W(t) is continuous.The first part of the definition describes the fact that the processW(t) has independent increments. The second part means that the increment W(s+t)-W(s) follows a normal distribution with mean zero and variance t. The third part describes the fact that W(t),t≥0, almost certainly has continuous paths. For an in-depth discussion about the Brownian motion process and its properties, see Durrett [17].Then, we letX(t) denote the real-time monitoring signal as a continuous stochastic process, with respect to time t. We define the functional form of X(t) as(2)Xt=φ+θexp⁡βt+εt-σ2t2=φ+θexp⁡βtexp⁡εt-σ2t2,where φ is the fixed intercept and represents the initial degradation, and β is a normal random variable such that the mean of β is μ1 and the variance of β is σ12. θ is a lognormal random variable with mean μ0 and variance σ02, and ε(t)=σW(t) is a Brownian motion with mean zero and variance σ2t. Under the assumption that θ, β, and ε(t) are mutually independent, it is obvious to obtain that Eexp⁡εt-σ2t/2=1, and thus E(S(t)∣θ,β)=φ+θexp⁡(βt).Furthermore, we find that it is easy to calculate with the logarithmic degradation data. Thus, we defineS(t) as follows:(3)St=ln⁡Xt-φ=ln⁡θ+βt+εt-σ2t2.By defining θ′=ln⁡θ, β′=β-σ2/2, we can further simplify S(t) as follows:(4)St=θ′+β′t+εt.


### 2.2. Bayesian Updating Process of Stochastic Parameters

We letSk=Stk-Stk-1,k=2,3,… denote the difference value of the monitored logged degradation signal at times tk and tk-1, with S1=S(t1).Then, suppose we have obtained the logged difference valueS1:k={S1,S2,…,Sk} at times t1,…,tk. And the error increments, ε(ti)-ε(ti-1),i=2,…,k, are independent normal random variables. If the stochastic parameters, θ′ and β′, are given, we can define the conditional joint density function of S1:k as(5)pS1:k∣θ′,β′=12πσ2kexp⁡-S1-θ′-β′t122σ2t1-∑i=2kSi-β′ti-ti-122σ2ti-ti-1.Generally, however,θ′ and β′ will be unknown. Based on the former assumption, we suppose we have got accurate and informative priors. And we let π0(θ′) and π1(β′) denote the prior distributions on θ′ and β′, respectively, where π0(θ′)=N(μ0,σ02) and π1(β′)=N(μ1′,σ12), μ1′=μ1-σ2/2. Then, given the logged difference data, S1:k, obtained at times t1,…,tk, the posterior joint distribution of (θ′,β′) can be expressed as follows.Theorem 2.
Given the obtained logged degradation signal,S1:k, the conditional posterior joint distribution of (θ′,β′) is θ′,β′∣S1:k~N(μθ′,k,σθ′,k2,μβ′,k,σβ′,k2,ρk), where(6)μθ′,k=S1σ02+μ0σ2t1σ12tk+σ2-σ02t1σ12∑i=1kSi+μ1′σ2σ02+σ2t1σ12tk+σ2-σ02σ12t1,μβ′,k=σ12∑i=1kSi+μ1′σ2σ2t1+σ02-σ12S1σ02+μ0σ2t1σ02+σ2t1σ12tk+σ2-σ02σ12t1,σθ′,k2=σ2σ02t1σ12tk+σ2σ02+σ2t1σ12tk+σ2-σ02σ12t1,σβ′,k2=σ2σ12σ02+σ2t1σ02+σ2t1σ12tk+σ2-σ02σ12t1,ρk=-σ0σ1t1σ02+σ2t1σ12tk+σ2.Proof.
Given the prior distributions ofθ′ and β′, π0(θ′)=N(μ0,σ02) and π1(β′)=N(μ1′,σ12), we can express the conditional posterior joint distribution of (θ′,β′), p(θ′,β′∣S1:k), as follows:(7)pθ′,β′∣S1:k∝pS1:k∣θ′,β′π0θ′π1β′∝12πσθ′,kσβ′,k1-ρk2·exp⁡-121-ρk2·θ′-μθ′,k2σθ′,k2-2ρkθ′-μθ′,kβ′-μβ′,kσθ′,kσβ′,k+β′-μβ′,k2σβ′,k2.As our focus is on the fusion algorithm, for a further detailed calculation of the posterior distribution of (θ′,β′), see Gebraeel et al. [10].


### 2.3. Estimation of the Residual Life Distribution

Every time a new degradation signal comes, we can compute a new posterior distribution of(θ′,β′). As the objective of our paper is to estimate the distribution of the RUL of the monitored system, we suppose that the system’s failure occurs when the observed degradation signal reaches the failure threshold, w, and thus we need to estimate the time until the degradation signal reaches w. In our paper, we assume that the threshold value is a constant value.The objective of prognostics is to compute the distribution of the failure time until the degradation signal reaches the thresholdw. To achieve this goal, we first calculate the posterior distribution of (θ′,β′). Then, we let the random variable S(t+tk) denote the logged degradation signal value obtained at time t+tk,t>0, given S1:k obtained at times t1,…,tk. Under the above assumption, the distribution of S(t+tk) given S1:k can be expressed as follows.Theorem 3.
Given the observed difference value of the logged degradation data,S1:k, the distribution of S(t+tk) is S(t+tk)∣S1:k~N(μ~(t+tk),σ~2(t+tk)), where(8)μ~t+tk≜∑i=1kSi+μβ′t=Stk+μβ′,kt,σ~2t+tk≜σβ′,k2t2+σ2t.Proof.
First, note thatSt=θ′+β′t+εt and we can express St+tk=Stk+β′t+εt+tk–εtk, where Stk=∑i=1kSi. Therefore, given S1:k, St+tk follows a normal distribution with mean μ~t+tk=Stk+Eβ′t=Stk+μβ′,kt and variance σ~2t+tk=t2Vβ′+Vεt+tk-εtk=σβ′,k2t2+σ2t.Then, we letT denote the RUL of the monitored system, and we know that T meets ST+tk=ln⁡w. Given S1:k, the conditional cumulative distribution function (CDF) of T, FT∣S1:k(t)=pT≤t∣S1:k, can be expressed as(9)PT≤t∣S1:k=1-PSt+tk≤ln⁡w∣S1:k=1-PZ<ln⁡w-μ~t+tkσ~2t+tk=PZ≥ln⁡w-μ~t+tkσ~2t+tk=Φgt,where Z represents a standard normal variable, Φ(∙) denotes the CDF of the standard normal variable, and g(t)=(μ~(t+tk)-ln⁡w)/σ~2(t+tk).We compute the residual life distribution of the system at timetk, under the condition that the observed degradation signal does not reach the threshold w; that is, ∑i=1kSi=Stk<ln⁡w. Thus, we get(10)limt→0⁡gt=limt→0⁡∑i=1kSi+μβ′,kt-ln⁡wσβ′,k2t2+σ2t=-∞.Therefore, we get limt→0⁡FT∣S1:k(t)=0, which means that the domain of the RUL, T, is (0,∞). We can express the conditional probability distribution function (PDF) of T, given S1:k, as(11)fT∣S1:k,T≥0t=ϕgtg′t,where ϕ(·) denotes the PDF of the standard normal variable.Given the conditional PDF ofT, we can write the expectation of RUL, at time tk, as(12)ET∣S1:k,T≥0t=∫0∞fT∣S1:k,T≥0tdt.In order to simplify the integral of expectation of RUL, we use the failure equationS(T+tk)=ln⁡w to compute the RUL, at time tk, approximately. And according to Theorem 3, given the difference value of logged degradation signal, S1:k, the distribution of S(t+tk) follows S(t+tk)∣S1:k~N(μ~(t+tk),σ~2(t+tk)), and the mean of this distribution, μ~(t+tk), closely approximates S(t+tk). Thus, we can write μ~(T+tk)=ln⁡w and express the RUL, at time tk, as(13)RULk=ln⁡w-Stkμβ′,k,where RULk is closely approximated to the expectation of RUL.In other literature, the prior distribution parameters,σ2, μ0, σ02, μ1′, σ12, are computed from historical monitoring data or derived from domain knowledge and expert experience. However, in this section we only have the real-time condition degradation information and do not know the accurate prior distribution parameters. What is more, the priors are fixed in the whole Bayesian updating procedure in the previous articles, and if the prior distribution parameters are inaccurate, the posterior distribution for these unknown parameters would have great errors. So our paper will solve these problems in the next section.


## 2.1. The Degradation Signal Model

First, we review the general definition of the Brownian motion process.Definition 1.
A standard Brownian motion process,W(t),t≥0, possesses the following properties:1.
Ift0<t1<⋯<tn, then Wt0,W(t1)-W(t0),…,W(tn)-W(tn-1) are mutually independent.2.
Ifs,t≥0, then(1)PWs+t-Ws∈A=∫A2πt-1/2e-x2/2tdx.3.
With probability onet→W(t) is continuous.The first part of the definition describes the fact that the processW(t) has independent increments. The second part means that the increment W(s+t)-W(s) follows a normal distribution with mean zero and variance t. The third part describes the fact that W(t),t≥0, almost certainly has continuous paths. For an in-depth discussion about the Brownian motion process and its properties, see Durrett [17].Then, we letX(t) denote the real-time monitoring signal as a continuous stochastic process, with respect to time t. We define the functional form of X(t) as(2)Xt=φ+θexp⁡βt+εt-σ2t2=φ+θexp⁡βtexp⁡εt-σ2t2,where φ is the fixed intercept and represents the initial degradation, and β is a normal random variable such that the mean of β is μ1 and the variance of β is σ12. θ is a lognormal random variable with mean μ0 and variance σ02, and ε(t)=σW(t) is a Brownian motion with mean zero and variance σ2t. Under the assumption that θ, β, and ε(t) are mutually independent, it is obvious to obtain that Eexp⁡εt-σ2t/2=1, and thus E(S(t)∣θ,β)=φ+θexp⁡(βt).Furthermore, we find that it is easy to calculate with the logarithmic degradation data. Thus, we defineS(t) as follows:(3)St=ln⁡Xt-φ=ln⁡θ+βt+εt-σ2t2.By defining θ′=ln⁡θ, β′=β-σ2/2, we can further simplify S(t) as follows:(4)St=θ′+β′t+εt.


## 2.2. Bayesian Updating Process of Stochastic Parameters

We letSk=Stk-Stk-1,k=2,3,… denote the difference value of the monitored logged degradation signal at times tk and tk-1, with S1=S(t1).Then, suppose we have obtained the logged difference valueS1:k={S1,S2,…,Sk} at times t1,…,tk. And the error increments, ε(ti)-ε(ti-1),i=2,…,k, are independent normal random variables. If the stochastic parameters, θ′ and β′, are given, we can define the conditional joint density function of S1:k as(5)pS1:k∣θ′,β′=12πσ2kexp⁡-S1-θ′-β′t122σ2t1-∑i=2kSi-β′ti-ti-122σ2ti-ti-1.Generally, however,θ′ and β′ will be unknown. Based on the former assumption, we suppose we have got accurate and informative priors. And we let π0(θ′) and π1(β′) denote the prior distributions on θ′ and β′, respectively, where π0(θ′)=N(μ0,σ02) and π1(β′)=N(μ1′,σ12), μ1′=μ1-σ2/2. Then, given the logged difference data, S1:k, obtained at times t1,…,tk, the posterior joint distribution of (θ′,β′) can be expressed as follows.Theorem 2.
Given the obtained logged degradation signal,S1:k, the conditional posterior joint distribution of (θ′,β′) is θ′,β′∣S1:k~N(μθ′,k,σθ′,k2,μβ′,k,σβ′,k2,ρk), where(6)μθ′,k=S1σ02+μ0σ2t1σ12tk+σ2-σ02t1σ12∑i=1kSi+μ1′σ2σ02+σ2t1σ12tk+σ2-σ02σ12t1,μβ′,k=σ12∑i=1kSi+μ1′σ2σ2t1+σ02-σ12S1σ02+μ0σ2t1σ02+σ2t1σ12tk+σ2-σ02σ12t1,σθ′,k2=σ2σ02t1σ12tk+σ2σ02+σ2t1σ12tk+σ2-σ02σ12t1,σβ′,k2=σ2σ12σ02+σ2t1σ02+σ2t1σ12tk+σ2-σ02σ12t1,ρk=-σ0σ1t1σ02+σ2t1σ12tk+σ2.Proof.
Given the prior distributions ofθ′ and β′, π0(θ′)=N(μ0,σ02) and π1(β′)=N(μ1′,σ12), we can express the conditional posterior joint distribution of (θ′,β′), p(θ′,β′∣S1:k), as follows:(7)pθ′,β′∣S1:k∝pS1:k∣θ′,β′π0θ′π1β′∝12πσθ′,kσβ′,k1-ρk2·exp⁡-121-ρk2·θ′-μθ′,k2σθ′,k2-2ρkθ′-μθ′,kβ′-μβ′,kσθ′,kσβ′,k+β′-μβ′,k2σβ′,k2.As our focus is on the fusion algorithm, for a further detailed calculation of the posterior distribution of (θ′,β′), see Gebraeel et al. [10].


## 2.3. Estimation of the Residual Life Distribution

Every time a new degradation signal comes, we can compute a new posterior distribution of(θ′,β′). As the objective of our paper is to estimate the distribution of the RUL of the monitored system, we suppose that the system’s failure occurs when the observed degradation signal reaches the failure threshold, w, and thus we need to estimate the time until the degradation signal reaches w. In our paper, we assume that the threshold value is a constant value.The objective of prognostics is to compute the distribution of the failure time until the degradation signal reaches the thresholdw. To achieve this goal, we first calculate the posterior distribution of (θ′,β′). Then, we let the random variable S(t+tk) denote the logged degradation signal value obtained at time t+tk,t>0, given S1:k obtained at times t1,…,tk. Under the above assumption, the distribution of S(t+tk) given S1:k can be expressed as follows.Theorem 3.
Given the observed difference value of the logged degradation data,S1:k, the distribution of S(t+tk) is S(t+tk)∣S1:k~N(μ~(t+tk),σ~2(t+tk)), where(8)μ~t+tk≜∑i=1kSi+μβ′t=Stk+μβ′,kt,σ~2t+tk≜σβ′,k2t2+σ2t.Proof.
First, note thatSt=θ′+β′t+εt and we can express St+tk=Stk+β′t+εt+tk–εtk, where Stk=∑i=1kSi. Therefore, given S1:k, St+tk follows a normal distribution with mean μ~t+tk=Stk+Eβ′t=Stk+μβ′,kt and variance σ~2t+tk=t2Vβ′+Vεt+tk-εtk=σβ′,k2t2+σ2t.Then, we letT denote the RUL of the monitored system, and we know that T meets ST+tk=ln⁡w. Given S1:k, the conditional cumulative distribution function (CDF) of T, FT∣S1:k(t)=pT≤t∣S1:k, can be expressed as(9)PT≤t∣S1:k=1-PSt+tk≤ln⁡w∣S1:k=1-PZ<ln⁡w-μ~t+tkσ~2t+tk=PZ≥ln⁡w-μ~t+tkσ~2t+tk=Φgt,where Z represents a standard normal variable, Φ(∙) denotes the CDF of the standard normal variable, and g(t)=(μ~(t+tk)-ln⁡w)/σ~2(t+tk).We compute the residual life distribution of the system at timetk, under the condition that the observed degradation signal does not reach the threshold w; that is, ∑i=1kSi=Stk<ln⁡w. Thus, we get(10)limt→0⁡gt=limt→0⁡∑i=1kSi+μβ′,kt-ln⁡wσβ′,k2t2+σ2t=-∞.Therefore, we get limt→0⁡FT∣S1:k(t)=0, which means that the domain of the RUL, T, is (0,∞). We can express the conditional probability distribution function (PDF) of T, given S1:k, as(11)fT∣S1:k,T≥0t=ϕgtg′t,where ϕ(·) denotes the PDF of the standard normal variable.Given the conditional PDF ofT, we can write the expectation of RUL, at time tk, as(12)ET∣S1:k,T≥0t=∫0∞fT∣S1:k,T≥0tdt.In order to simplify the integral of expectation of RUL, we use the failure equationS(T+tk)=ln⁡w to compute the RUL, at time tk, approximately. And according to Theorem 3, given the difference value of logged degradation signal, S1:k, the distribution of S(t+tk) follows S(t+tk)∣S1:k~N(μ~(t+tk),σ~2(t+tk)), and the mean of this distribution, μ~(t+tk), closely approximates S(t+tk). Thus, we can write μ~(T+tk)=ln⁡w and express the RUL, at time tk, as(13)RULk=ln⁡w-Stkμβ′,k,where RULk is closely approximated to the expectation of RUL.In other literature, the prior distribution parameters,σ2, μ0, σ02, μ1′, σ12, are computed from historical monitoring data or derived from domain knowledge and expert experience. However, in this section we only have the real-time condition degradation information and do not know the accurate prior distribution parameters. What is more, the priors are fixed in the whole Bayesian updating procedure in the previous articles, and if the prior distribution parameters are inaccurate, the posterior distribution for these unknown parameters would have great errors. So our paper will solve these problems in the next section.


## 3. Estimating the Prior Information by Fusion Algorithm

In Section2, we have estimated the RUL of the single component by Bayesian updating procedure under the assumption that a priori distributions were known. However, we only have real-time observations and the prior distributions needed in Bayesian updating process are nonexistent. So our paper will develop a collaboration algorithm between Bayesian updating and ECM algorithm to estimate these a priori parameters.


### 3.1. The ECM Algorithm

We letΘ=[σ2,μ0,σ02,μ1′,σ12] denote the unknown prior distribution parameters. Given the difference value of logged degradation signal S1:k, we can express the log-likelihood function of S1:k as(14)lkΘ=ln⁡pS1:k∣Θ=∑j=2kln⁡pSj∣S1:j-1,Θ,where p(S1:k∣Θ) is the joint PDF of S1:k. Thus, we can write the MLE of Θ at time tk as(15)Θ^=arg⁡maxΘ⁡lkΘ,where Θ^ is the variable value of Θ corresponding to maximum of lk(Θ). Our goal is to estimate an appropriate prior distribution parameter Θ^. However, in formula (4) θ′ and β′ are stochastic and unobserved, so the calculation of formula (15) is hard to complete.In order to avoid the above problem, we propose a collaboration algorithm between Bayesian updating and ECM algorithm and use the fusion algorithm to estimate prior distributions based on degradation signals. We letΘk=[σk2,μ0,k,μ1,k′,σ0,k2,σ1,k2] denote the estimated parameters based on degradation signals, S1:k, and Θ^k(i)=[σ^k2(i),μ^0,k(i),μ^1,k′(i),σ^0,k2(i),σ^1,k2(i)] denote the estimated variable value of Θk of the ith iteration result in the ECM algorithm. We treat θ′ and β′ as the hidden variables, because θ′ and β′ are stochastic and unobserved in the whole degradation process. Our objective is to calculate the variable value of Θk, so the E-step of ECM algorithm can be written as (16)lΘk∣Θ^ki=Eθ′,β′∣S1:k,Θ^kiln⁡pS1:k,θ′,β′∣Θk,where Eθ′,β′∣S1:k,Θ^k(i) is the conditional expectation of hidden variables, θ′ and β′, given Θ^ki and S1:k.Next we express the CM-step of ECM algorithm as(17)Θ^ki+1/S=arg⁡maxΘ⁡lΘk∣Θ^ki=σ^k2i,μ^0,ki,μ^1,k′i,σ^0,k2i,σ^1,k2i,Θ^ki+2/S=arg⁡maxΘ⁡lΘk∣Θ^ki+1/S=σ^k2i+1,μ^0,ki,μ^1,k′i,σ^0,k2i,σ^1,k2i,Θ^ki+3/S=arg⁡maxΘ⁡lΘk∣Θ^ki+2/S=σ^k2i+1,μ^0,ki+1,μ^1,k′i,σ^0,k2i,σ^1,k2i,Θ^ki+4/S=arg⁡maxΘ⁡lΘk∣Θ^ki+3/S=σ^k2i+1,μ^0,ki+1,μ^1,k′i+1,σ^0,k2i,σ^1,k2i⋮Θ^ki+1=arg⁡maxΘ⁡lΘk∣Θ^ki+S-1/S=σ^k2i+1,μ^0,ki+1,μ^1,k′i+1,σ^0,k2i+1,σ^1,k2i,where S(S=5) is the dimension of vector Θ, and Θ^k(i+1) is the maximum value of Θ based on Θ^k(i+(S-1)/S) from the derivative equation ∂l(Θk∣Θ^k(i+(S-1)/S))/∂σ^1,k2(i)=0. As our focus is on the algorithm fusion process, for the theory of ECM algorithm and its convergence analysis, see Meng and Rubin [18], Van et al. [19], and Liu and Rubin [20].


### 3.2. The Collaboration between Bayesian Updating and ECM Algorithm

After we have reviewed the ECM algorithm, we begin our fusion process. Figure2 shows the procedure of the fusion process.Figure 2
The collaboration algorithm between Bayesian updating and ECM algorithm.First, we use random initial a priori distributions to start the fusion algorithm, when we collect the degradation signalS(t1) at time t1. Then, we get the posterior distribution of (θ′,β′), p(θ′,β′∣S1:k), and a residual life distribution by Bayesian updating method. Of course, the results are inaccurate.Then, we use the posterior distribution of(θ′,β′), p(θ′,β′∣S1:k), to substitute the distribution, given Θ^k(i), of hidden variables in the E-step. So, the E-step of the ECM algorithm can be rewritten as(18)lΘk∣pθ′,β′∣S1:k=Eθ′,β′∣pθ′,β′∣S1:kln⁡pS1:k,θ′,β′∣Θk.Because of the different E-step, we get a rewritten CM-step which is different from the one in the ECM algorithm. In order to get a more accurate estimated value, Θ^k, in ECM algorithm we need multiple iterations. However, we can get the optimal solution of Θk only through S(S=5) steps calculation, in our fusion algorithm. This will be proved in Theorem 4. In our fusion algorithm the rewritten CM-step can be expressed as(19)Θ^k1/5=arg⁡maxΘ⁡lΘk∣pθ′,β′∣S1:k=σk2,μ0,k,μ1,k′,σ0,k2,σ1,k2,Θ^k2/5=arg⁡maxΘ⁡lΘk∣pθ′,β′∣S1:k1/5=σ^k2,μ0,k,μ1,k′,σ0,k2,σ1,k2,Θ^k3/5=arg⁡maxΘ⁡lΘk∣pθ′,β′∣S1:k2/5=σ^k2,μ^0,k,μ1,k′,σ0,k2,σ1,k2,Θ^k4/5=arg⁡maxΘ⁡lΘk∣pθ′,β′∣S1:k3/5=σ^k2,μ^0,k,μ^1,k′,σ0,k2,σ1,k2,Θ^k=arg⁡maxΘ⁡lΘk∣pθ′,β′∣S1:k4/5=σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ1,k2.In this five-step calculation of formula (19), we can get the optimal estimations σk2=σ^k2, μ0,k=μ^0,k, μ1,k′=μ^1,k′, σ0,k2=σ^0,k2, σ1,k2=σ^1,k2, respectively. The results of the rewritten CM-step are the a priori distributions of the next Bayesian updating procedure. Each time we collect a new degradation signal S(tk) at time tk, we can recalculate the a priori distributions and the residual life distribution.


### 3.3. Estimation of Prior Information

After we have finished the fusion process, we begin to calculate the optimal variable value ofΘk by our fusion algorithm. The log-likelihood function of complete data can be expressed as(20)ln⁡pS1:k,θ′,β′∣Θk=ln⁡pS1:k∣θ′,β′,Θk+ln⁡pθ′,β′∣Θk=-k+22ln⁡2π-k2ln⁡σk2-S1-θ′-β′t122σk2t1-∑i=2kSi-β′ti-ti-122σk2ti-ti-1-12ln⁡σ0,k2-12ln⁡σ1,k2-θ′-μ0,k22σ0,k2-β′-μ1,k′22σ1,k2.According to formulas (18) and (19), we can compute the rewritten E-step as follows:(21)lΘk∣pθ′,β′∣S1:k=Eθ′,β′∣pθ′,β′∣S1:kln⁡pS1:k,θ′,β′∣Θk=-k+22ln⁡2π-k2ln⁡σk2-S12-2S1μθ′,k+μβ′,kt1+μθ′,k2+σθ′,k2+2t1ρkσθ′,kσβ′,k+μθ′,kμβ′,k+t12μβ′,k2+σβ′,k22σk2t1-∑i=2kSi2-2Siμβ′,kti-ti-1+μβ′,k2+σβ′,k2ti-ti-122σk2ti-ti-1-12ln⁡σ0,k2-12ln⁡σ1,k2-μθ′,k2+σθ′,k2-2μθ′,kμ0,k+μ0,k22σ0,k2-μβ′,k2+σβ′,k2-2μβ′,kμ1,k′+μ1,k′22σ1,k2.Then, based on the rewritten CM-step, we can get the optimal σ^k2, μ^0,k, μ^1,k′, σ^0,k2, σ^1,k2 by the following steps, respectively.Step 1.
By calculating∂l(Θk∣p(θ′,β′∣S1:k)/∂σk2=0, we can get the optimal(22)σ^k2=1kS12-2S1μθ′,k+μβ′,kt1+μθ′,k2+σθ′,k2+2t1ρkσθ′,kσβ′,k+μθ′,kμβ′,k+t12μβ′,k2+σβ′,k2t1+∑i=2kSi2-2Siμβ′,kti-ti-1+μβ′,k2+σβ′,k2ti-ti-12ti-ti-1.Step 2.
By calculating∂l(Θk∣p(θ′,β′∣S1:k)(1/5)/∂μ0,k=0, we can get the optimal μ^0,k=μθ′,k.Step 3.
By calculating∂l(Θk∣p(θ′,β′∣S1:k)(2/5)/∂μ1,k′=0, we can get the optimal μ^1,k′=μβ′,k.Step 4.
By calculating∂l(Θk∣p(θ′,β′∣S1:k)(3/5)/∂σ0,k2=0, we can get the optimal σ^0,k2=σθ′,k2.Step 5.
By calculating∂l(Θk∣p(θ′,β′∣S1:k)(4/5)/∂σ1,k2=0, we can get the optimal σ^1,k2=σβ′,k2.According to the above five steps,Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2] can be expressed as(23)σ^k2=1kS12-2S1μθ′,k+μβ′,kt1+μθ′,k2+σθ′,k2+2t1ρkσθ′,kσβ′,k+μθ′,kμβ′,k+t12μβ′,k2+σβ′,k2t1+∑i=2kSi2-2Siμβ′,kti-ti-1+μβ′,k2+σβ′,k2ti-ti-12ti-ti-1,μ^0,k=μθ′,k,σ^0,k2=σθ′,k2,μ^1,k′=μβ′,k,σ^1,k2=σβ′,k2.We get the optimal σ^k2, μ^0,k, μ^1,k′, σ^0,k2, σ^1,k2 in each step, respectively. However, it does not mean that Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2] is the only maximum point of ∂l(Θk∣p(θ′,β′∣S1:k)/∂Θk=0. So, we need a further discussion.Theorem 4.
According to the result of formula (23), the estimated Θ^k is the only maximum point of ∂l(Θk∣p(θ′,β′∣S1:k)/∂Θk=0.Proof.
First, from formula (21) we can get Θ^k which is the only solution of ∂l(Θk∣p(θ′,β′∣S1:k)/∂Θk=0. Next, we prove that Θ^k is corresponding to maximum l(Θk∣p(θ′,β′∣S1:k). We can perform the second derivative of log-likelihood function on vector Θk as (24)∂2lΘk∣pθ′,β′∣S1:k∂Θk∂ΘkT=k2σk4-τσk600000-1σ0,k20μ0,k-μθ′,kσ0,k4000-1σ1,k20μ1,k′-μβ′,kσ1,k40μ0,k-μθ′,kσ0,k4012σ0,k4-ψ1σ0,k6000μ1,k′-μβ′,kσ1,k4012σ1,k4-ψ2σ1,k6,where(25)τ=S12-2S1μθ′,k+μβ′,kt1+μθ′,k2+σθ′,k2+2t1ρkσθ′,kσβ′,k+μθ′,kμβ′,k+t12μβ′,k2+σβ′,k2t1+∑i=2kSi2-2Siμβ′,kti-ti-1+μβ′,k2+σβ′,k2ti-ti-12ti-ti-1,ψ1=μθ′,k2+σθ′,k2-2μθ′,kμ0,k+μ0,k2,ψ2=μβ′,k2+σβ′,k2-2μβ′,kμ1,k′+μ1,k′2.Then we use the estimated variable value Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2] in (23), and we can find the order principal minor of the matrix in (24) as follows:(26)Δ1Θk=Θ^k<0,Δ2Θk=Θ^k>0,Δ3Θk=Θ^k<0,Δ4Θk=Θ^k>0,Δ5Θk=Θ^k<0.Thus, we can conclude that (24) is a negative-definite matrix when Θk=Θ^k, which means that Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2] is the only maximum point.This fusion procedure between Bayesian updating and ECM algorithm can be performed each time a new degradation signal is observed. That is to say, each time the degradation signalS(tk) is observed, we can recalculate Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2] and obtain new estimates of residual life for the newly made system. What is more, the initial values of priors in the Bayesian updating for the first time are unrestricted, and once we get a new Θ^k, it would be used in the next Bayesian updating as the priors. In Section 4, we will use simulation method to further evaluate the performance of our collaboration algorithm and illustrate that the unrestricted initial values of priors have little influence on the estimated vector Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2].


## 3.1. The ECM Algorithm

We letΘ=[σ2,μ0,σ02,μ1′,σ12] denote the unknown prior distribution parameters. Given the difference value of logged degradation signal S1:k, we can express the log-likelihood function of S1:k as(14)lkΘ=ln⁡pS1:k∣Θ=∑j=2kln⁡pSj∣S1:j-1,Θ,where p(S1:k∣Θ) is the joint PDF of S1:k. Thus, we can write the MLE of Θ at time tk as(15)Θ^=arg⁡maxΘ⁡lkΘ,where Θ^ is the variable value of Θ corresponding to maximum of lk(Θ). Our goal is to estimate an appropriate prior distribution parameter Θ^. However, in formula (4) θ′ and β′ are stochastic and unobserved, so the calculation of formula (15) is hard to complete.In order to avoid the above problem, we propose a collaboration algorithm between Bayesian updating and ECM algorithm and use the fusion algorithm to estimate prior distributions based on degradation signals. We letΘk=[σk2,μ0,k,μ1,k′,σ0,k2,σ1,k2] denote the estimated parameters based on degradation signals, S1:k, and Θ^k(i)=[σ^k2(i),μ^0,k(i),μ^1,k′(i),σ^0,k2(i),σ^1,k2(i)] denote the estimated variable value of Θk of the ith iteration result in the ECM algorithm. We treat θ′ and β′ as the hidden variables, because θ′ and β′ are stochastic and unobserved in the whole degradation process. Our objective is to calculate the variable value of Θk, so the E-step of ECM algorithm can be written as (16)lΘk∣Θ^ki=Eθ′,β′∣S1:k,Θ^kiln⁡pS1:k,θ′,β′∣Θk,where Eθ′,β′∣S1:k,Θ^k(i) is the conditional expectation of hidden variables, θ′ and β′, given Θ^ki and S1:k.Next we express the CM-step of ECM algorithm as(17)Θ^ki+1/S=arg⁡maxΘ⁡lΘk∣Θ^ki=σ^k2i,μ^0,ki,μ^1,k′i,σ^0,k2i,σ^1,k2i,Θ^ki+2/S=arg⁡maxΘ⁡lΘk∣Θ^ki+1/S=σ^k2i+1,μ^0,ki,μ^1,k′i,σ^0,k2i,σ^1,k2i,Θ^ki+3/S=arg⁡maxΘ⁡lΘk∣Θ^ki+2/S=σ^k2i+1,μ^0,ki+1,μ^1,k′i,σ^0,k2i,σ^1,k2i,Θ^ki+4/S=arg⁡maxΘ⁡lΘk∣Θ^ki+3/S=σ^k2i+1,μ^0,ki+1,μ^1,k′i+1,σ^0,k2i,σ^1,k2i⋮Θ^ki+1=arg⁡maxΘ⁡lΘk∣Θ^ki+S-1/S=σ^k2i+1,μ^0,ki+1,μ^1,k′i+1,σ^0,k2i+1,σ^1,k2i,where S(S=5) is the dimension of vector Θ, and Θ^k(i+1) is the maximum value of Θ based on Θ^k(i+(S-1)/S) from the derivative equation ∂l(Θk∣Θ^k(i+(S-1)/S))/∂σ^1,k2(i)=0. As our focus is on the algorithm fusion process, for the theory of ECM algorithm and its convergence analysis, see Meng and Rubin [18], Van et al. [19], and Liu and Rubin [20].


## 3.2. The Collaboration between Bayesian Updating and ECM Algorithm

After we have reviewed the ECM algorithm, we begin our fusion process. Figure2 shows the procedure of the fusion process.Figure 2
The collaboration algorithm between Bayesian updating and ECM algorithm.First, we use random initial a priori distributions to start the fusion algorithm, when we collect the degradation signalS(t1) at time t1. Then, we get the posterior distribution of (θ′,β′), p(θ′,β′∣S1:k), and a residual life distribution by Bayesian updating method. Of course, the results are inaccurate.Then, we use the posterior distribution of(θ′,β′), p(θ′,β′∣S1:k), to substitute the distribution, given Θ^k(i), of hidden variables in the E-step. So, the E-step of the ECM algorithm can be rewritten as(18)lΘk∣pθ′,β′∣S1:k=Eθ′,β′∣pθ′,β′∣S1:kln⁡pS1:k,θ′,β′∣Θk.Because of the different E-step, we get a rewritten CM-step which is different from the one in the ECM algorithm. In order to get a more accurate estimated value, Θ^k, in ECM algorithm we need multiple iterations. However, we can get the optimal solution of Θk only through S(S=5) steps calculation, in our fusion algorithm. This will be proved in Theorem 4. In our fusion algorithm the rewritten CM-step can be expressed as(19)Θ^k1/5=arg⁡maxΘ⁡lΘk∣pθ′,β′∣S1:k=σk2,μ0,k,μ1,k′,σ0,k2,σ1,k2,Θ^k2/5=arg⁡maxΘ⁡lΘk∣pθ′,β′∣S1:k1/5=σ^k2,μ0,k,μ1,k′,σ0,k2,σ1,k2,Θ^k3/5=arg⁡maxΘ⁡lΘk∣pθ′,β′∣S1:k2/5=σ^k2,μ^0,k,μ1,k′,σ0,k2,σ1,k2,Θ^k4/5=arg⁡maxΘ⁡lΘk∣pθ′,β′∣S1:k3/5=σ^k2,μ^0,k,μ^1,k′,σ0,k2,σ1,k2,Θ^k=arg⁡maxΘ⁡lΘk∣pθ′,β′∣S1:k4/5=σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ1,k2.In this five-step calculation of formula (19), we can get the optimal estimations σk2=σ^k2, μ0,k=μ^0,k, μ1,k′=μ^1,k′, σ0,k2=σ^0,k2, σ1,k2=σ^1,k2, respectively. The results of the rewritten CM-step are the a priori distributions of the next Bayesian updating procedure. Each time we collect a new degradation signal S(tk) at time tk, we can recalculate the a priori distributions and the residual life distribution.


## 3.3. Estimation of Prior Information

After we have finished the fusion process, we begin to calculate the optimal variable value ofΘk by our fusion algorithm. The log-likelihood function of complete data can be expressed as(20)ln⁡pS1:k,θ′,β′∣Θk=ln⁡pS1:k∣θ′,β′,Θk+ln⁡pθ′,β′∣Θk=-k+22ln⁡2π-k2ln⁡σk2-S1-θ′-β′t122σk2t1-∑i=2kSi-β′ti-ti-122σk2ti-ti-1-12ln⁡σ0,k2-12ln⁡σ1,k2-θ′-μ0,k22σ0,k2-β′-μ1,k′22σ1,k2.According to formulas (18) and (19), we can compute the rewritten E-step as follows:(21)lΘk∣pθ′,β′∣S1:k=Eθ′,β′∣pθ′,β′∣S1:kln⁡pS1:k,θ′,β′∣Θk=-k+22ln⁡2π-k2ln⁡σk2-S12-2S1μθ′,k+μβ′,kt1+μθ′,k2+σθ′,k2+2t1ρkσθ′,kσβ′,k+μθ′,kμβ′,k+t12μβ′,k2+σβ′,k22σk2t1-∑i=2kSi2-2Siμβ′,kti-ti-1+μβ′,k2+σβ′,k2ti-ti-122σk2ti-ti-1-12ln⁡σ0,k2-12ln⁡σ1,k2-μθ′,k2+σθ′,k2-2μθ′,kμ0,k+μ0,k22σ0,k2-μβ′,k2+σβ′,k2-2μβ′,kμ1,k′+μ1,k′22σ1,k2.Then, based on the rewritten CM-step, we can get the optimal σ^k2, μ^0,k, μ^1,k′, σ^0,k2, σ^1,k2 by the following steps, respectively.Step 1.
By calculating∂l(Θk∣p(θ′,β′∣S1:k)/∂σk2=0, we can get the optimal(22)σ^k2=1kS12-2S1μθ′,k+μβ′,kt1+μθ′,k2+σθ′,k2+2t1ρkσθ′,kσβ′,k+μθ′,kμβ′,k+t12μβ′,k2+σβ′,k2t1+∑i=2kSi2-2Siμβ′,kti-ti-1+μβ′,k2+σβ′,k2ti-ti-12ti-ti-1.Step 2.
By calculating∂l(Θk∣p(θ′,β′∣S1:k)(1/5)/∂μ0,k=0, we can get the optimal μ^0,k=μθ′,k.Step 3.
By calculating∂l(Θk∣p(θ′,β′∣S1:k)(2/5)/∂μ1,k′=0, we can get the optimal μ^1,k′=μβ′,k.Step 4.
By calculating∂l(Θk∣p(θ′,β′∣S1:k)(3/5)/∂σ0,k2=0, we can get the optimal σ^0,k2=σθ′,k2.Step 5.
By calculating∂l(Θk∣p(θ′,β′∣S1:k)(4/5)/∂σ1,k2=0, we can get the optimal σ^1,k2=σβ′,k2.According to the above five steps,Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2] can be expressed as(23)σ^k2=1kS12-2S1μθ′,k+μβ′,kt1+μθ′,k2+σθ′,k2+2t1ρkσθ′,kσβ′,k+μθ′,kμβ′,k+t12μβ′,k2+σβ′,k2t1+∑i=2kSi2-2Siμβ′,kti-ti-1+μβ′,k2+σβ′,k2ti-ti-12ti-ti-1,μ^0,k=μθ′,k,σ^0,k2=σθ′,k2,μ^1,k′=μβ′,k,σ^1,k2=σβ′,k2.We get the optimal σ^k2, μ^0,k, μ^1,k′, σ^0,k2, σ^1,k2 in each step, respectively. However, it does not mean that Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2] is the only maximum point of ∂l(Θk∣p(θ′,β′∣S1:k)/∂Θk=0. So, we need a further discussion.Theorem 4.
According to the result of formula (23), the estimated Θ^k is the only maximum point of ∂l(Θk∣p(θ′,β′∣S1:k)/∂Θk=0.Proof.
First, from formula (21) we can get Θ^k which is the only solution of ∂l(Θk∣p(θ′,β′∣S1:k)/∂Θk=0. Next, we prove that Θ^k is corresponding to maximum l(Θk∣p(θ′,β′∣S1:k). We can perform the second derivative of log-likelihood function on vector Θk as (24)∂2lΘk∣pθ′,β′∣S1:k∂Θk∂ΘkT=k2σk4-τσk600000-1σ0,k20μ0,k-μθ′,kσ0,k4000-1σ1,k20μ1,k′-μβ′,kσ1,k40μ0,k-μθ′,kσ0,k4012σ0,k4-ψ1σ0,k6000μ1,k′-μβ′,kσ1,k4012σ1,k4-ψ2σ1,k6,where(25)τ=S12-2S1μθ′,k+μβ′,kt1+μθ′,k2+σθ′,k2+2t1ρkσθ′,kσβ′,k+μθ′,kμβ′,k+t12μβ′,k2+σβ′,k2t1+∑i=2kSi2-2Siμβ′,kti-ti-1+μβ′,k2+σβ′,k2ti-ti-12ti-ti-1,ψ1=μθ′,k2+σθ′,k2-2μθ′,kμ0,k+μ0,k2,ψ2=μβ′,k2+σβ′,k2-2μβ′,kμ1,k′+μ1,k′2.Then we use the estimated variable value Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2] in (23), and we can find the order principal minor of the matrix in (24) as follows:(26)Δ1Θk=Θ^k<0,Δ2Θk=Θ^k>0,Δ3Θk=Θ^k<0,Δ4Θk=Θ^k>0,Δ5Θk=Θ^k<0.Thus, we can conclude that (24) is a negative-definite matrix when Θk=Θ^k, which means that Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2] is the only maximum point.This fusion procedure between Bayesian updating and ECM algorithm can be performed each time a new degradation signal is observed. That is to say, each time the degradation signalS(tk) is observed, we can recalculate Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2] and obtain new estimates of residual life for the newly made system. What is more, the initial values of priors in the Bayesian updating for the first time are unrestricted, and once we get a new Θ^k, it would be used in the next Bayesian updating as the priors. In Section 4, we will use simulation method to further evaluate the performance of our collaboration algorithm and illustrate that the unrestricted initial values of priors have little influence on the estimated vector Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2].


## 4. Simulation and Analysis

In this section, we will adopt the simulation method to further evaluate the performance of our collaboration algorithm. First, in order to represent a degradation process, we create a set of simulation data based on the exponential degradation model in Section2. In the simulation, we assume that φ=0, θ′~N(0.02,2×10-6), β′~N(0.01,1×10-6), ε(t)~N(0,t·4×10-4). In order to observe enough degradation signals and obvious degradation trend, we let the threshold w=60 and the sampling interval (ti-ti-1)=10,i=2,3,4,…. Figure 3 shows the trajectory of the simulated degradation signal. We obtain 41 degradation samples, and the degradation reaches the standard threshold, w=60, at time t=403.5.Figure 3
Simulation of degradation path.We know thatθ′~N(μ0,σ02) and β′~N(μ1-σ2/2,σ12), and μ0 and μ1-σ2/2 dominate the degradation and RUL. We will use the estimated results of them to prove that these unrestricted initial priors have little influence on the estimated accuracy of the RUL. The different initial prior distributions are as shown in Table 1.Table 1
Different initial prior distributions.
Data
θ
′
β
′
ε
(
t
)
(1)
N
(
0.1,2
×
1
0
-
4
)
N
(
0.05,1
×
1
0
-
4
)
N
(
0
,
t
⋅
4
×
1
0
-
3
)
(2)
N
(
0.2,2
×
1
0
-
4
)
N
(
0.1,1
×
1
0
-
4
)
N
(
0
,
t
⋅
4
×
1
0
-
3
)
(3)
N
(
0.3,2
×
1
0
-
4
)
N
(
0.15,1
×
1
0
-
4
)
N
(
0
,
t
⋅
4
×
1
0
-
3
)On the other hand, in order to prove that our collaboration algorithm can get a more accurate RUL than other separate Bayesian updating methods, we will compare our fusion algorithm with the method in Gebraeel et al. [10]. by using the 2nd set of initial priors in Table 1. Figure 4 represents the estimated results of μ0. Figure 5 represents the estimated results of μ1-σ2/2.Figure 4
Estimated results ofμ0.Figure 5
Estimated results ofμ1-σ2/2.From Figures4 and 5, we know that the collaboration algorithm can estimate the mean of θ′ and β′ accurately even with inaccurate prior information. Although we have used a variety of inaccurate prior distributions, it can still be observed that this aspect has little effect on the estimation of the means of θ′ and β′, showing that the degradation measurements dominate and the collaboration algorithm is effective. We can know that the estimated means of θ′ and β′ gradually approximate to the simulated value, 0.02 and 0.01, at time about 50, in spite of different initial priors.What is more, our collaboration algorithm can get more accurate estimations compared with the separate Bayesian updating by using the 2nd set of initial priors in Table1.Given the estimated prior distributions we calculate the point estimations of RUL of the newly made component. Figure6 shows the point estimations of RUL by our algorithm and separate Bayesian updating method.Figure 6
Point estimations of RUL by collaboration algorithm and separate Bayesian Updating method.From Figure6, we can know that our collaboration algorithm can get a more accurate RUL comparing with the separate Bayesian method. And inaccurate a priori distributions have little effect on the estimation of the RUL. Furthermore, the point estimations of RUL by our collaboration algorithm can also reflect the fluctuation of degradation caused by Brownian motion error, and this can be known from Figure 3.


## 5. Conclusion

In this paper, we presented a collaboration algorithm that contains the characteristic of Bayesian updating and ECM algorithm. The difficulties of the fusion process mainly consist of two parts; the first is building the connection between ECM algorithm and Bayesian updating. It is not easy to find another substitute for theith iteration distribution in the rewritten E-step, and the posterior distribution of the Bayesian updating is an optimal one. The second is proving the optimal estimating result about the prior information. We can get a maximum estimating result, Θ^k=[σ^k2,μ^0,k,μ^1,k′,σ^0,k2,σ^1,k2], by the rewritten CM-step. However, it does not mean that the estimating result is the only maximum point of ∂l(Θk∣p(θ′,β′∣S1:k)/∂Θk=0. In this paper, we use the second derivative of log-likelihood function on vector Θk and the order principal minor of the matrix in (24) to prove it.Our fusion algorithm can predict failures of newly made systems in the absence of prior degradation knowledge. Although our fusion algorithm started with random initial a priori distributions, the simulation experiments show that the inaccurate a priori distributions have little effect on the estimation of the RUL, and our fusion algorithm can get a better prediction than the separate Bayesian method. Nevertheless, there are still some issues that needed further investigation for the estimation of the RUL on the single component without prior degradation knowledge.First, we assume that the exponential degradation signal model and the error fluctuation are adequate for the given degradation path. Actually, we can not guarantee which model is the best one by visual judgment. So we also need to focus on model selection and evaluate the goodness-of-fit of various degradation path models.Second, we assume that the stochastic parameters of the exponential degradation model are normally distributed. However, the stochastic parameters may follow a Gamma distribution or other distributions. Therefore, we also should investigate the performance of our collaboration algorithm when the underlying normal distribution assumptions are not satisfied.Third, our work started the RUL calculation in the degradation processes assuming that the exact point-in-time of the initial degradation is known. However, the functioning system would be stable within a period of time, and the point-in-time for the initial degradation is unknown and stochastic. By considering the distribution of the initial degradation point, we may be able to predict the RUL right after its installment.



---

*Source: 4375690-2017-12-03.xml*