# A Nonlinear Prognostic Model Based on the Wiener Process with Three Sources of Uncertainty

Authors: Huifang Niu; Jianchao Zeng; Hui Shi; Bin Wang; Tianye Liu
Journal: Shock and Vibration
 (2021)

---

## Abstract

Estimation of the remaining useful life (RUL) is an important component of prognostics and health management (PHM). The accuracy of the RUL estimation for complex systems is mainly affected by three sources of uncertainty, i.e., the temporal uncertainty, the product-to-product uncertainty, and measurement errors. To improve PHM and account for the effects of the three sources of uncertainty, a nonlinear prognostic model with three sources of uncertainty is presented here. An approximated analytical expression for the probability density function (PDF) of the RUL is obtained based on the concept of first hitting time (FHT). Model parameters are then obtained by the expectation maximization (EM) algorithm, and the drift parameter is estimated adaptively using a Bayesian procedure. Finally, in order to illustrate the practical applications of the presented approach, a comparative study of real data on fatigue crack propagation is presented. Results demonstrate that our method improves model fit and increases the accuracy of the lifetime estimation.


---

## Body

## 1. Introduction

The development of condition monitoring technology makes it easier to obtain the degradation of a system [1]. In the past decade, research about estimation of the RUL based on condition monitor data had made great progress. Accurate RUL estimation allows decision makers to make maintenance decisions as planned, improving overall system reliability and reducing the occurrence of failures and maintenance costs [2–5]. However, for most complex systems or components, such as gyros [6], bearing [7–12], LED lamps [13–16], and battery systems [2, 17–19], the degradation process is stochastic, and thus, the accurate RUL estimation becomes very challenging.In practical applications, it is indispensable to establish a degradation model to determine the RUL due to stochastic degradation of the system. To achieve this goal, we notice that the RUL of a system is mainly affected by three sources of uncertainty [3, 20, 21], i.e., the temporal uncertainty, the product-to-product uncertainty, and measurement errors. Therefore, the estimated accuracy of the RUL may be improved by fully accounting for the three sources of uncertainty in the degrading system.Most of the existing approaches to Wiener process prediction take into account the influence of one or two sources of uncertainty on the RUL of a degrading system. For example, in [22], a Wiener process with adaptive drift for the RUL estimation was presented, and the impact of the temporal uncertainty on the PDF of the RUL was considered. The approach from [22] was improved [23] by introducing a Wiener process model based on recursive filtering algorithm, while considering the influence of the temporal uncertainty and the product-to-product uncertainty on the PDF of the RUL. In [20], the diffusion process corresponding to a nonlinear drift coefficient with constant threshold was transformed by a time-space transformation [24] into a standard Brownian motion with a variable threshold. An analytic approximation was also obtained for the RUL distribution, and the impact of the temporal uncertainty and the product-to-product uncertainty on the PDF of the RUL was considered. The uncertainty of measurement error on the RUL was ignored, and the parameters were estimated by historical degradation data without real-time online update. Similar study can also be found in other papers [25–35]. However, there are several models based on a Wiener process that take into account the influence of three sources of uncertainty on the RUL. In [36], a Wiener process model with three sources of uncertainty has been presented, and the RUL distribution has been derived. However, the model is limited to linear degradation process, while in practical engineering applications, the degradation of a system is nonlinear, in particular, at the late stage, when degradation is accelerating. In other words, a linear Wiener process is not suitable to describe the actual degradation of a system. In [21, 37, 38], models were presented which include three sources of uncertainty of degradation, which were used to derive the RUL distribution. However, the unknown parameters are obtained by the MLE method with unobserved hidden variables, which leads to inaccuracies in the estimation results. Hao et al. [39] presented a nonlinear stepped stress degradation model with three sources of uncertainty, which has the advantage of requiring small sample size and short test time. However, the method is unable to reproduce the actual failure process of the system, and it is thus obviously not suitable for complex systems with high reliability.The conclusion to be drawn from the above overview is that current Wiener process models for prediction only partially includes the effects of the temporal uncertainty, the product-to-product uncertainty, and measurement errors, and thus, the true dynamics of system degradation cannot be completely described. In addition, in the existing literature about Wiener process prognostic models with three sources of uncertainty, the unknown parameters of the model have been obtained by the MLE method, which is not suitable when hidden variables are contained in the likelihood function. Therefore, a more accurate RUL prognostic model is needed that takes into account the impact of three sources of uncertainty on the PDF of the RUL.Inspired by the above studies, in this paper, we present a nonlinear prognostic model with three sources of uncertainty and derive an approximate analytic expression for the PDF of the RUL. By combining a Bayesian approach with the EM algorithm, unknown parameters of the model are estimated. In particular, the drift parameter is adaptively estimated by a Bayesian method, and other parameters are obtained by the EM algorithm.The main innovations of this article are the following. (1) A nonlinear prognostic model which takes into account simultaneously the temporal uncertainty, the product-to-product uncertainty, and the measurement errors is proposed. (2)An approximate analytic expression for the RUL distribution, which incorporates the effects of the three sources of uncertainty is obtained. (3)The unknown parameters of the model are adaptively estimated using a suitable combination of Bayesian methods and the EM algorithm.The rest of this paper is organized as follows. In Section2, we present the nonlinear prognostic model with three sources of uncertainty and derive the approximate analytic form of the PDF of the RUL. In Section 3, the parameter estimation approach of the prognostic model is presented. In Section 4, we study a practical case to demonstrate the effectiveness of the presented method. Section 5 concludes with a summary.


## 2. Nonlinear Wiener Process Model and RUL Estimation

### 2.1. Nonlinear Wiener Process Model

Generally speaking, a nonlinear Wiener process can be expressed as(1)Xt=X0+λ∫0tμτ;θdτ+σBt.It is further assumed that∫0tμτ;θdτ=Λt;θ; then, we have(2)Xt=X0+λΛt;θ+σBt,where λ is the drift coefficient, which is used to represent the product-to-product uncertainty, Λt;θ is a function of time t and the parameter vector θ, σ is the diffusion coefficient, Bt is the standard Brownian motion (BM) process, used to describe the temporal uncertainty of the degradation process, and X0=x0 is the initial degeneration state of the system. In the following, we assume X0=x0=0.Due to the differences in the material and production processes of the system, the rate of degradation in each unit degradation process is different. Usually, a random parameterλ is called rate of degradation of different units, whereas θ and σ are fixed parameters, representing the common degradation features of all systems in the population. We assume that λ follows a normal distribution, i.e., λ∼Nμλ,σλ2, and λ and Bt are mutually independent. Similar assumptions were made in [20] and [25, 26, 40].In addition, since measurements are affected by fluctuations due to the external environmental noise and to nonideal behavior of detectors, the measurement dataYt cannot reflect the actual degradation state Xt of the system accurately. The relationship between them may be expressed as(3)Yt=Xt+ε,where ε describes the random measurement errors and ε∼N0,σε2. λ, Bt, and ε are assumed mutually independent. Equation (3) has been used in literature, see e.g., [21, 36, 39]. Equations (2) and (3) are used to describe a nonlinear prognostic model with three sources of uncertainty.


### 2.2. RUL Estimation

As in other prognostic models [3, 15, 36] and according to the definition of the FHT, if one does not consider measurement errors, the system life T may be expressed as(4)T=inft:Xt≥ω|X0<ω.Similarly, the RUL of a system at timetk is given by(5)Lk=inflk:Xtk+lk≥ω|Xt<ω,∀t≤tk.At first, for the degradation data without measurement errors, we have the following lemma [20].Lemma 1.
For the degradation processXt in (1), if the function μt;θ is continuous in t∈0,+∞, under mild assumptions [20], then the PDF of the lifetime of Xt crossing the constant boundary ω is given by(6)ft|λ≅12πtMtt+λσμt;θexp−M2t2t,where Mt=1/σω−λΛt;θ
Based on Lemma1, the PDF of the RUL of a system at time tk and degradation state xk is(7)flk|λ,xk=ωk−Λ∗lk;θλ+lkμlk+tk;θλσ2πlk3exp−ωk−Λ∗lk;θλ22lkσ2,where Λ∗lk;θ=Λlk+tk;θ−Λtk;θ and ωk=ω−xk.
If one considers measurement errors, the FHTTe of Yt crossing threshold ω may be defined as(8)Tε=inft:Yt≥ω|Y0<ω=inft:Xt≥ωε|X0<ωε,where ωε=ω−ε and ωε∼Nω,σε2.
Similarly, the RULLk of a system with measurement errors at time tk is expressed as follows:(9)Lεk=inflk:Ytk+lk≥ω|Yt<ω,∀t≤tk=inflk:Xtk+lk≥ωεXt<ωε∀t≤tk,where ωε,k=ωε−xk follows a normal distribution [13, 14], i.e., ωε,k∼Nω−yk,2σε2.
In the following, we are going to deduce the unconditional PDF of the FHT. As mentioned above,λ and ωε are random variables that are independent of each other. By considering the impact of λ and ωε on the PDF of the lifetime, the PDF of Tε is obtained using the law of total probability, i.e.,(10)ft=∫−∞+∞∫−∞+∞fTε|λ,ωεtpωε|λpλdωεdλ=EλEωε|λfTε|λ,ωεt.
To calculate the integral in equation (10) exactly, the following lemma holds [20].Lemma 2.
Ifα∼Nμα,σα2 and a,b,c,d∈R, r∈R+, then we have(11)Eαa−bα⋅exp−c−dα22r=rd2σα2+r×a−bdσα2c+μαrd2σα2+rexp−c−dμα22d2σα2+r.
Now, based on Lemma1 and 2, we can obtain the PDF of the FHT.Theorem 1.
Ifλ∼Nμλ,σλ2 and ωε∼Nω,σε2, then the PDF of the FHT for a degradation model considering measurement errors may be expressed as(12)ft=EλEωε|λfTε|λ,ωεt=1ηt2πt2Λ2t;θσλ2+ηt⋅exp−ω−Λt;θμλ22Λ2t;θσλ2+ηtωσ2t+ηttμt;θ−Λt;θ+σε2Λt;θ⋅Λt;θσλ2ω+μληtΛ2t;θσλ2+ηt,where ηt=σε2+σ2t.Proof.
By replacingω in E, (6) with ωε, and considering that ωe∼Nω,σε2, we have(13)Eωε|λfTε|λ,ωεt=Eωε|λ1σ2πt3ωε−λΛt;θ+μt;θλt⋅exp−ωε−Λt;θλ22σ2t=1σ2πt3Eωε|λωε−λΛt;θ+μt;θλt⋅exp−ωε−Λt;θλ22σ2t=12πt2σε2+σ2t⋅μt;θt−Λt;θλ+σε2Λt;θλ+ωσ2tσε2+σ2t⋅exp−ωε−Λt;θλ22σε2+σ2t.
Now, according to Lemma2, the result is proved by setting a=Λt;θ−μt;θtλ, b=d=1, c=Λt;θλ, and r=σ2t. Furthermore, when the drift parameter λ∼Nμλ,σλ2 is taken into account, we have(14)EλEωε|λfTε|λ,ωεt=Eλ12πt2σε2+σ2tμt;θt−Λt;θλ+σε2Λt;θλ+ωσ2tσε2+σ2t⋅exp−Λt;θλ−ω22σε2+σ2t=12πt2σε2+σ2tEλ|Y1:kμt;θt−Λt;θλ+σε2Λt;θλ+ωσ2tσε2+σ2t⋅exp−Λt;θλ−ω22σε2+σ2t=1σε2+σ2t2πt2Λ2t;θσλ2+σε2+σ2tωσ2t+σε2+σ2ttμt;θ−Λt;θ+σε2Λt;θ⋅Λt;θσλ2ω+μλσε2+σ2tΛ2t;θσλ2+σε2+σ2texp−ω−Λt;θμλ22Λ2t;θσλ2+σε2+σ2t.
On the basis of Lemma2, Theorem 1 is proved by setting a=ωσ2t, b=σε2+σ2tΛt;θ−tμt;θ−σε2Λt;θ, d=Λt;θ, and r=Λ2t;θσλ2+σε2+σ2t.
In the same way, we obtain the following theorem.Theorem 2.
Ifλ∼Nμi,σi2 and ωε,k∼Nω−yk,2σε2, given the current measured data yk, then the PDF of the RUL for the established model with measurement errors at time tk can be expressed as(15)fLk|Y1:klk|Y1:k=1δlk2πlk2Λ∗lk;θ2σi2+δlklkσ2ω−yk−δlkΛ∗lk;θ−lkμlk+t;kθ−2σε2Λ∗lk;θ⋅Λ∗lk;θσi2ω−yk+μiδlkΛ∗lk;θ2σi2+δlkexp−ω−yk−Λ∗lk;θμi22Λ∗lk;θ2σi2+δlk,where δlk=2σε2+lkσ2.
The proof of Theorem2 is similar to Theorem 1 and is omitted here.


## 2.1. Nonlinear Wiener Process Model

Generally speaking, a nonlinear Wiener process can be expressed as(1)Xt=X0+λ∫0tμτ;θdτ+σBt.It is further assumed that∫0tμτ;θdτ=Λt;θ; then, we have(2)Xt=X0+λΛt;θ+σBt,where λ is the drift coefficient, which is used to represent the product-to-product uncertainty, Λt;θ is a function of time t and the parameter vector θ, σ is the diffusion coefficient, Bt is the standard Brownian motion (BM) process, used to describe the temporal uncertainty of the degradation process, and X0=x0 is the initial degeneration state of the system. In the following, we assume X0=x0=0.Due to the differences in the material and production processes of the system, the rate of degradation in each unit degradation process is different. Usually, a random parameterλ is called rate of degradation of different units, whereas θ and σ are fixed parameters, representing the common degradation features of all systems in the population. We assume that λ follows a normal distribution, i.e., λ∼Nμλ,σλ2, and λ and Bt are mutually independent. Similar assumptions were made in [20] and [25, 26, 40].In addition, since measurements are affected by fluctuations due to the external environmental noise and to nonideal behavior of detectors, the measurement dataYt cannot reflect the actual degradation state Xt of the system accurately. The relationship between them may be expressed as(3)Yt=Xt+ε,where ε describes the random measurement errors and ε∼N0,σε2. λ, Bt, and ε are assumed mutually independent. Equation (3) has been used in literature, see e.g., [21, 36, 39]. Equations (2) and (3) are used to describe a nonlinear prognostic model with three sources of uncertainty.


## 2.2. RUL Estimation

As in other prognostic models [3, 15, 36] and according to the definition of the FHT, if one does not consider measurement errors, the system life T may be expressed as(4)T=inft:Xt≥ω|X0<ω.Similarly, the RUL of a system at timetk is given by(5)Lk=inflk:Xtk+lk≥ω|Xt<ω,∀t≤tk.At first, for the degradation data without measurement errors, we have the following lemma [20].Lemma 1.
For the degradation processXt in (1), if the function μt;θ is continuous in t∈0,+∞, under mild assumptions [20], then the PDF of the lifetime of Xt crossing the constant boundary ω is given by(6)ft|λ≅12πtMtt+λσμt;θexp−M2t2t,where Mt=1/σω−λΛt;θ
Based on Lemma1, the PDF of the RUL of a system at time tk and degradation state xk is(7)flk|λ,xk=ωk−Λ∗lk;θλ+lkμlk+tk;θλσ2πlk3exp−ωk−Λ∗lk;θλ22lkσ2,where Λ∗lk;θ=Λlk+tk;θ−Λtk;θ and ωk=ω−xk.
If one considers measurement errors, the FHTTe of Yt crossing threshold ω may be defined as(8)Tε=inft:Yt≥ω|Y0<ω=inft:Xt≥ωε|X0<ωε,where ωε=ω−ε and ωε∼Nω,σε2.
Similarly, the RULLk of a system with measurement errors at time tk is expressed as follows:(9)Lεk=inflk:Ytk+lk≥ω|Yt<ω,∀t≤tk=inflk:Xtk+lk≥ωεXt<ωε∀t≤tk,where ωε,k=ωε−xk follows a normal distribution [13, 14], i.e., ωε,k∼Nω−yk,2σε2.
In the following, we are going to deduce the unconditional PDF of the FHT. As mentioned above,λ and ωε are random variables that are independent of each other. By considering the impact of λ and ωε on the PDF of the lifetime, the PDF of Tε is obtained using the law of total probability, i.e.,(10)ft=∫−∞+∞∫−∞+∞fTε|λ,ωεtpωε|λpλdωεdλ=EλEωε|λfTε|λ,ωεt.
To calculate the integral in equation (10) exactly, the following lemma holds [20].Lemma 2.
Ifα∼Nμα,σα2 and a,b,c,d∈R, r∈R+, then we have(11)Eαa−bα⋅exp−c−dα22r=rd2σα2+r×a−bdσα2c+μαrd2σα2+rexp−c−dμα22d2σα2+r.
Now, based on Lemma1 and 2, we can obtain the PDF of the FHT.Theorem 1.
Ifλ∼Nμλ,σλ2 and ωε∼Nω,σε2, then the PDF of the FHT for a degradation model considering measurement errors may be expressed as(12)ft=EλEωε|λfTε|λ,ωεt=1ηt2πt2Λ2t;θσλ2+ηt⋅exp−ω−Λt;θμλ22Λ2t;θσλ2+ηtωσ2t+ηttμt;θ−Λt;θ+σε2Λt;θ⋅Λt;θσλ2ω+μληtΛ2t;θσλ2+ηt,where ηt=σε2+σ2t.Proof.
By replacingω in E, (6) with ωε, and considering that ωe∼Nω,σε2, we have(13)Eωε|λfTε|λ,ωεt=Eωε|λ1σ2πt3ωε−λΛt;θ+μt;θλt⋅exp−ωε−Λt;θλ22σ2t=1σ2πt3Eωε|λωε−λΛt;θ+μt;θλt⋅exp−ωε−Λt;θλ22σ2t=12πt2σε2+σ2t⋅μt;θt−Λt;θλ+σε2Λt;θλ+ωσ2tσε2+σ2t⋅exp−ωε−Λt;θλ22σε2+σ2t.
Now, according to Lemma2, the result is proved by setting a=Λt;θ−μt;θtλ, b=d=1, c=Λt;θλ, and r=σ2t. Furthermore, when the drift parameter λ∼Nμλ,σλ2 is taken into account, we have(14)EλEωε|λfTε|λ,ωεt=Eλ12πt2σε2+σ2tμt;θt−Λt;θλ+σε2Λt;θλ+ωσ2tσε2+σ2t⋅exp−Λt;θλ−ω22σε2+σ2t=12πt2σε2+σ2tEλ|Y1:kμt;θt−Λt;θλ+σε2Λt;θλ+ωσ2tσε2+σ2t⋅exp−Λt;θλ−ω22σε2+σ2t=1σε2+σ2t2πt2Λ2t;θσλ2+σε2+σ2tωσ2t+σε2+σ2ttμt;θ−Λt;θ+σε2Λt;θ⋅Λt;θσλ2ω+μλσε2+σ2tΛ2t;θσλ2+σε2+σ2texp−ω−Λt;θμλ22Λ2t;θσλ2+σε2+σ2t.
On the basis of Lemma2, Theorem 1 is proved by setting a=ωσ2t, b=σε2+σ2tΛt;θ−tμt;θ−σε2Λt;θ, d=Λt;θ, and r=Λ2t;θσλ2+σε2+σ2t.
In the same way, we obtain the following theorem.Theorem 2.
Ifλ∼Nμi,σi2 and ωε,k∼Nω−yk,2σε2, given the current measured data yk, then the PDF of the RUL for the established model with measurement errors at time tk can be expressed as(15)fLk|Y1:klk|Y1:k=1δlk2πlk2Λ∗lk;θ2σi2+δlklkσ2ω−yk−δlkΛ∗lk;θ−lkμlk+t;kθ−2σε2Λ∗lk;θ⋅Λ∗lk;θσi2ω−yk+μiδlkΛ∗lk;θ2σi2+δlkexp−ω−yk−Λ∗lk;θμi22Λ∗lk;θ2σi2+δlk,where δlk=2σε2+lkσ2.
The proof of Theorem2 is similar to Theorem 1 and is omitted here.


## 3. Statistical Inference

In this section, we describe the methods we used to determine the function and to estimate the unknown parameters and in equations (2) and (3) based on the available degradation data.


### 3.1. The Determination of the FunctionΛt;θ

Before determining parameters of the model, we need to know the functionΛt;θ. Following the approach [37], we know that the form of Λt;θ can be obtained based on the prior knowledge or the degradation physics. Even when we cannot understand the degradation physics about the system, degradation data usually contain enough information to determine the form of the function Λt;θ. Indeed, the function Λt;θ may be obtained by analyzing the degenerate measurement data. For example, in this paper, exploiting the relationship between the cycle number and the fatigue crack size, the function Λt;θ may be obtained easily as EYt=λΛt;θ.The specific steps to determine the functionΛt;θ are the following. At first, the degradation data for all units are averaged at each sampling time. Afterwards, the sample average y¯1,y¯2,…,y¯m is fitted with a large class of possible functional forms, including the exponential function and power laws. The fitting function represents Λt;θ, where y¯j=1/n∑i=1nYitj and Yitj is the measurement data of unit i at the jth time, i=1,2,…,n and j=1,2,…,mi.


### 3.2. The Parameter Estimation

Let us now address the estimation of the parametersμλ,σλ2,σ2,σε2, and θ in equations (2) and (3). We denote by Θ=μ0,σ02,σ2,σε2,θ the vector of all unknown parameters. We assume that there are n independent units for each testing with Yitj denoting the measurement data for the ith unit at time tj, i=1,2,…,nandj=1,2,…mi. Here, we assume that the degradation data of each unit is measured in the same time sequence t<1t2<…<tm. According to equation (2) and (3), the degradation measurement of the unit i can be expressed as Yitj=λiΛtj;θ+σBtj+εij , where λi is the drift coefficient of the ith unit. The variables λi are i.i.d of λ and εij are i.i.d realizations of ε. The increment of the degradation measurement for ith unit is expressed as ΔYi=yi1,yi2−yi1,…,yij−yi,j−1,…,i=1,2,…,n,j=1,2,…,mT and the increments of the degradation measurements for the n test units are denoted as ΔY=ΔY1,ΔY2,…,ΔYn. According to the nature of the Wiener processes, when λi is fixed, ΔYi∼NλiΛ,Σ, and the joint PDF of ΔYi may be expressed as(16)pΔYi|λi,Θ=2π−m2Σ−12exp−12ΔYi−λiΛTΣ−1ΔYi−λiΛ,where the elements of the covariance matrix Σ are given by(17)Σij=σ2Δti+σε2,i=j=1,σ2Δti+2σε2,i=j>1,−σε2,i−j=1,0,otherwise.In order to simplify the inference procedure and the estimation of parameters, we use the re-parameterizationσ˜ε2=σε2/σ2 and Σ˜σ˜ε2,θ=Σ/σ2, and Ω=λ1,λ2,…,λm denotes the vector of drift coefficients. Therefore, given the measurements Y and Ω, the log-likelihood function of Θ is(18)ℓΘ|ΔY,Ω=∑i=1nlnfΔYi|λi,Θ+lnpλi|Θ=C−n2lnΣ˜σ˜ε2,θ−nm2lnσ2−n2lnσλ2−12∑i=1nΔYi−λiΛTΣ˜−1σ˜ε2,θΔYi−λiΛσ2+λi−μλ2σλ2,where C represents a generic constant. Then, the MLE of Θ is represented by Θ^, which is given by Θ^=argmaxΘℓΘ.It is impossible to maximize (18) directly since the drift parameters λi are an unobserved hidden variable. However, the EM algorithm provides a way to cope with hidden variables [14, 25, 41]. The EM algorithm is a (converging) iteration process, which incorporates the update of Θ at each step. Let us denote by Θk=μλk,σλk,σ2k,σε2k,θk the estimate at the kth step based on ΔY.As mentioned above,λi∼Nμλ,σλ2. According to a Bayesian approach, the posterior distribution of λi can be written as follows:(19)pλi|ΔYi,Θk=pΔYi|λi,Θkpλi|ΘkpΔYi|Θk∝pΔYi|λi,Θkpλi|Θk∝exp−12σ2kΔYi−λiΛTΣ˜−1σ˜ε2k,θkΔYi−λiΛexp−λi−μλk22σλ2k∝exp−12ΛTΣ˜−1σ˜ε2k,θkΛσ2k+1σλ2kλi2−2ΔYiTΣ˜−1σ˜ε2k,θkΛσ2k+μλkσλ2kλi∝exp−λi−σλ2kΔYiTΣ˜−1σ˜ε2k,θkΛ+μλkσ2k/σλ2kΛTΣ˜−1σ˜ε2k,θkΛ+σ2k22σ2kσλ2k/σλ2kΛTΣ˜−1σ˜ε2k,θkΛ+σ2k∼Nμik,σi2k,μik=σλ2kΔYiTΣ˜−1σ˜ε2k,θkΛ+μλkσ2kσλ2kΛTΣ˜−1σ˜ε2k,θkΛ+σ2k,σi2k=σ2kσλ2kσλ2kΛTΣ˜−1σε2k,θkΛ+σ2k,uik=Eλi2|ΔYi=σλ2k+μλk2.Using the above derivation, it can be seen that the drift parameterμi,σi2 of the ith unit can be adjusted adaptively to make its uncertainty smaller. Next, we obtain the MLE of Θ using the EM algorithm iteratively. The detailed description of the EM algorithm may be found [42].The expectationQΘ|ΔY,Θk of ℓΘ|ΔY,Ω can be written as(20)QΘ|ΔY,Θk=EΩ|ΘkℓΘ|ΔY,Ω=C−12∑i=1nuik−2μikμλ+μλ2σλ2+lnσλ2−12∑i=1nΔYiTΣ˜−1σ˜ε2,θΔYi−2μikΔYiTΣ˜−1σ˜ε2,θΛ+uikΛTΣ˜−1σ˜ε2,θΛσ2+lnΣ˜σ˜ε2,θ+mlnσ2.Setting∂QΘ|ΔY,Θk/∂μλ=0, ∂QΘ|ΔY,Θk/∂σλ2=0, and ∂QΘ|ΔY,Θk/∂σ2=0, we have(21)μλk+1=1n∑i=1nμik,σλ2k+1=1n∑i=1nuik−μλ2k+1,(22)σ2k+1=1nm∑i=1nΔYiTϕ−1σ˜ε2,θΔYi−2μikΔYiTϕ−1σ˜ε2,θΛ+uikΛTϕ−1σ˜ε2,θΛ.It is easily found thatσ2k+1 depends on σ˜ε2 and θ. Next, by substituting μλk+1, σλ2k+1, and σ2k+1 in equation (20), we can obtain the profile likelihood function:(23)ℓσ˜ε2,θ|Θk=C−12∑i=1nlnΣ˜σ˜ε2,θ+mln1nm∑i=1nΔYiTΣ˜−1σ˜ε2,θΔYi−2μikΔYiTΣ˜−1σ˜ε2,θΛ+uikΛTΣ˜−1σ˜ε2,θΛ.The genetic algorithm (GA) is used to maximize (23) and obtain the MLE of σ˜ε2k+1 and θk+1. A more detailed overview of GA may be found in [43]. Next, the estimate of σ2 in the k+1th step may be obtained by substituting σ˜ε2k+1 and θk+1 into equation (22). Furthermore, σε2 can be estimated as σ˜ε2=σε2/σ2. The specific content of iteration termination was included in [42].


## 3.1. The Determination of the FunctionΛt;θ

Before determining parameters of the model, we need to know the functionΛt;θ. Following the approach [37], we know that the form of Λt;θ can be obtained based on the prior knowledge or the degradation physics. Even when we cannot understand the degradation physics about the system, degradation data usually contain enough information to determine the form of the function Λt;θ. Indeed, the function Λt;θ may be obtained by analyzing the degenerate measurement data. For example, in this paper, exploiting the relationship between the cycle number and the fatigue crack size, the function Λt;θ may be obtained easily as EYt=λΛt;θ.The specific steps to determine the functionΛt;θ are the following. At first, the degradation data for all units are averaged at each sampling time. Afterwards, the sample average y¯1,y¯2,…,y¯m is fitted with a large class of possible functional forms, including the exponential function and power laws. The fitting function represents Λt;θ, where y¯j=1/n∑i=1nYitj and Yitj is the measurement data of unit i at the jth time, i=1,2,…,n and j=1,2,…,mi.


## 3.2. The Parameter Estimation

Let us now address the estimation of the parametersμλ,σλ2,σ2,σε2, and θ in equations (2) and (3). We denote by Θ=μ0,σ02,σ2,σε2,θ the vector of all unknown parameters. We assume that there are n independent units for each testing with Yitj denoting the measurement data for the ith unit at time tj, i=1,2,…,nandj=1,2,…mi. Here, we assume that the degradation data of each unit is measured in the same time sequence t<1t2<…<tm. According to equation (2) and (3), the degradation measurement of the unit i can be expressed as Yitj=λiΛtj;θ+σBtj+εij , where λi is the drift coefficient of the ith unit. The variables λi are i.i.d of λ and εij are i.i.d realizations of ε. The increment of the degradation measurement for ith unit is expressed as ΔYi=yi1,yi2−yi1,…,yij−yi,j−1,…,i=1,2,…,n,j=1,2,…,mT and the increments of the degradation measurements for the n test units are denoted as ΔY=ΔY1,ΔY2,…,ΔYn. According to the nature of the Wiener processes, when λi is fixed, ΔYi∼NλiΛ,Σ, and the joint PDF of ΔYi may be expressed as(16)pΔYi|λi,Θ=2π−m2Σ−12exp−12ΔYi−λiΛTΣ−1ΔYi−λiΛ,where the elements of the covariance matrix Σ are given by(17)Σij=σ2Δti+σε2,i=j=1,σ2Δti+2σε2,i=j>1,−σε2,i−j=1,0,otherwise.In order to simplify the inference procedure and the estimation of parameters, we use the re-parameterizationσ˜ε2=σε2/σ2 and Σ˜σ˜ε2,θ=Σ/σ2, and Ω=λ1,λ2,…,λm denotes the vector of drift coefficients. Therefore, given the measurements Y and Ω, the log-likelihood function of Θ is(18)ℓΘ|ΔY,Ω=∑i=1nlnfΔYi|λi,Θ+lnpλi|Θ=C−n2lnΣ˜σ˜ε2,θ−nm2lnσ2−n2lnσλ2−12∑i=1nΔYi−λiΛTΣ˜−1σ˜ε2,θΔYi−λiΛσ2+λi−μλ2σλ2,where C represents a generic constant. Then, the MLE of Θ is represented by Θ^, which is given by Θ^=argmaxΘℓΘ.It is impossible to maximize (18) directly since the drift parameters λi are an unobserved hidden variable. However, the EM algorithm provides a way to cope with hidden variables [14, 25, 41]. The EM algorithm is a (converging) iteration process, which incorporates the update of Θ at each step. Let us denote by Θk=μλk,σλk,σ2k,σε2k,θk the estimate at the kth step based on ΔY.As mentioned above,λi∼Nμλ,σλ2. According to a Bayesian approach, the posterior distribution of λi can be written as follows:(19)pλi|ΔYi,Θk=pΔYi|λi,Θkpλi|ΘkpΔYi|Θk∝pΔYi|λi,Θkpλi|Θk∝exp−12σ2kΔYi−λiΛTΣ˜−1σ˜ε2k,θkΔYi−λiΛexp−λi−μλk22σλ2k∝exp−12ΛTΣ˜−1σ˜ε2k,θkΛσ2k+1σλ2kλi2−2ΔYiTΣ˜−1σ˜ε2k,θkΛσ2k+μλkσλ2kλi∝exp−λi−σλ2kΔYiTΣ˜−1σ˜ε2k,θkΛ+μλkσ2k/σλ2kΛTΣ˜−1σ˜ε2k,θkΛ+σ2k22σ2kσλ2k/σλ2kΛTΣ˜−1σ˜ε2k,θkΛ+σ2k∼Nμik,σi2k,μik=σλ2kΔYiTΣ˜−1σ˜ε2k,θkΛ+μλkσ2kσλ2kΛTΣ˜−1σ˜ε2k,θkΛ+σ2k,σi2k=σ2kσλ2kσλ2kΛTΣ˜−1σε2k,θkΛ+σ2k,uik=Eλi2|ΔYi=σλ2k+μλk2.Using the above derivation, it can be seen that the drift parameterμi,σi2 of the ith unit can be adjusted adaptively to make its uncertainty smaller. Next, we obtain the MLE of Θ using the EM algorithm iteratively. The detailed description of the EM algorithm may be found [42].The expectationQΘ|ΔY,Θk of ℓΘ|ΔY,Ω can be written as(20)QΘ|ΔY,Θk=EΩ|ΘkℓΘ|ΔY,Ω=C−12∑i=1nuik−2μikμλ+μλ2σλ2+lnσλ2−12∑i=1nΔYiTΣ˜−1σ˜ε2,θΔYi−2μikΔYiTΣ˜−1σ˜ε2,θΛ+uikΛTΣ˜−1σ˜ε2,θΛσ2+lnΣ˜σ˜ε2,θ+mlnσ2.Setting∂QΘ|ΔY,Θk/∂μλ=0, ∂QΘ|ΔY,Θk/∂σλ2=0, and ∂QΘ|ΔY,Θk/∂σ2=0, we have(21)μλk+1=1n∑i=1nμik,σλ2k+1=1n∑i=1nuik−μλ2k+1,(22)σ2k+1=1nm∑i=1nΔYiTϕ−1σ˜ε2,θΔYi−2μikΔYiTϕ−1σ˜ε2,θΛ+uikΛTϕ−1σ˜ε2,θΛ.It is easily found thatσ2k+1 depends on σ˜ε2 and θ. Next, by substituting μλk+1, σλ2k+1, and σ2k+1 in equation (20), we can obtain the profile likelihood function:(23)ℓσ˜ε2,θ|Θk=C−12∑i=1nlnΣ˜σ˜ε2,θ+mln1nm∑i=1nΔYiTΣ˜−1σ˜ε2,θΔYi−2μikΔYiTΣ˜−1σ˜ε2,θΛ+uikΛTΣ˜−1σ˜ε2,θΛ.The genetic algorithm (GA) is used to maximize (23) and obtain the MLE of σ˜ε2k+1 and θk+1. A more detailed overview of GA may be found in [43]. Next, the estimate of σ2 in the k+1th step may be obtained by substituting σ˜ε2k+1 and θk+1 into equation (22). Furthermore, σε2 can be estimated as σ˜ε2=σε2/σ2. The specific content of iteration termination was included in [42].


## 4. A Case Study

To assess the usefulness of our approach in practical applications, we now present a case study based on the fatigue crack data given in [20]. The results coming from our model are also compared to those obtained using some other existing models. We use the overall mean squared errors (TMSE) [20] of the fitting model to compare the fitting of our method to that of alternative modeling approaches.The mean squared error (MSE) which is usually used to directly evaluate the fit to the data is defined as(24)MSEk=∫0∞lk−l˜k2fLk|Y0:klk|Y0:kdlk,where lk and l˜k are the estimated RUL and the actual RUL at time tk, respectively, and fLk|Y0:klk|Y0:k is the conditional PDF of the RUL estimated by (15).In a situation where there arem measurements, the TMSE of the lifetime can be written as(25)TMSE=∑k=1mMSEk.According to the above criteria, we have that the smaller is the TMSE, the better is the accuracy of the model.To illustrate the effectiveness of our model, the following two competitors are considered.Case 1.
We refer to Si’s method, see [20], as Case 1. It considers both the temporal uncertainty and the product-to-product uncertainty. A detailed description can be found in [20].Case 2.
We refer to Tang’s method as Case2, see [38], for details. It considers the temporal uncertainty, the product-to-product uncertainty, and the measurement errors simultaneously. However, Tang’s method estimates the parameters of the model by MLE.
For simplicity, our method is denoted as Case 3.
The model fitting of these three cases is compared using the case of fatigue crack data for the 2017-T4 aluminum alloy [20]. The obtained data represent the crack length propagation of four 2017-T4 aluminum alloy specimens at the stress level of 200 MPa, simulating the actual operating conditions (see Figure 1). The fatigue crack length of all specimen are measured from 1.5 million to 2.4 million cycles at steps of 0.1 million cycles. The failed threshold is ω=5.6[20]. Generally speaking, the time at which the threshold is crossed cannot be measured exactly. Therefore, the FHT is approximated using the data from the first time at which the measured value exceeds the threshold. In this experiment, the average failure time was about 2.4 million cycles. In the case of small samples, the estimated results of our method are satisfactory in terms of the TMSE.
From Figure1, we can see that the degradation paths of the fatigue crack for four 2017-T4 aluminum alloy may be rather different, that is, there is a significant individual variability. Therefore, the stochastic process model of equation (2), driven by BM with nonlinear drift is appropriate to characterize the degradation process of fatigue crack propagation.
In the following, according to the method in Section3.1, we determine the function Λt;θ from the degradation measurements. At first, the mean of the degradation measurements for fatigue crack propagation data at each measuring point is calculated, as depicted in Figure 2. Then, we fit the sample mean by using the possible functions of Λt;θ. The calculated values of the correlation coefficients are shown in Table 1, from which an exponential Λt;θ is found. Assuming that the degradation measurement at t=0 is zero, we have Λt;θ=expθt−1. Then, from equations (2) and (3), we have that the estimated mean degradation path is μ^λΛt;θ. Intuitively, if the established model is correct, the estimated mean degradation path should be almost the same as the sample mean curve. From Figure 2, we can see that these two curves indeed nearly coincide, indicating the fitting of the established model.
From Figure3, it is apparent that the estimates of all the parameters in the present model quickly converge to a constant value, which illustrates that the EM algorithm proposed in Section 3.2 is effective. Additionally, it indicates that the parameter estimation method adopted in this paper is feasible. As expected, we can see from Table 2 that the estimated value of θ is greater than 1, which verifies the nonlinearity of the degradation path. In addition, results from Table 2 show that the fitting of case 3 is better than Cases 1-2 in terms of TMSE.
We illustrate the PDF of the lifetimes for Cases 1–3 with the fatigue crack data in Figure4.
From Figure4, it is clear that the PDF of the lifetimes of the fatigue crack data is completely different between Case 1 and Case 2. This is mostly due to the inability of Case 1 to fully take into account the uncertainty (measurement errors) in the degradation data. Similarly, a few differences of the PDF of the lifetimes are present between Case 2 and Case 3. The main reason here is the unobservability of λ in the likelihood function. When there are hidden variables in the likelihood function that cannot be observed, the EM algorithm used in case 3 provides better results.
To further demonstrate the practicability of the proposed method, the degradation data of the 3rd unit is used to obtain the PDFs of the RULs for case 1–case 3 at run time. Accordingly, for the convenience of comparison, the actual RULs and the PDFs of the RULs are also shown in Figure5.
In Figure5, the actual RULs are shown with solid line and asterisk marks. We see that the PDFs of the RULs in Case 3 are more peaked around the actual RULs than Case 1 and Case 2, and this illustrates that RUL estimation is much less uncertain than in the other two cases. In addition, we can see that the estimated results produced by Case 3 significantly outperform Case 1, which denotes the necessity of considering measurement error of degradation data. Finally, the estimated results in Case 3 are more accurate than Case 2. The main reason is that the parameters of the model obtained by the EM algorithm in Case 3 are closer to the actual values.
Based on the above results, we conclude that the model and the method of Case 3, which includes three source uncertainties, are more accurate than those of other cases.Figure 1
Degradation measurements of fatigue crack propagation data.Figure 2
Comparison of the sample mean curve and the estimated mean degradation path.Table 1
Correlation coefficients of different fitting functions.
Λt;θ typeCorrelation coefficientsExponential function0.9924Power law function0.9867Linear function0.8063Figure 3
Illustration of the iterative convergence using the EM parameter estimation of case 3.
(a)(b)(c)(d)(e)Table 2
Comparison of three cases with the fatigue crack propagation data.
μλσλ2σ2σε2θTMSECase 11.4955e-41.2468e-96.0652e-5—4.44020.0215Case 27.566e-46.8498e-228.8904e-40.19153.76500.0078Case 34.304e-41.8912e-112.4405e-30.19393.99250.0043Figure 4
Comparison of the PDF of the lifetime for Cases 1–3 of the fatigue crack propagation data.Figure 5
Comparison of the PDFs of the RULs for Cases 1–3 of the fatigue crack propagation data.


## 5. Conclusion

In this paper, we have proposed a nonlinear prognostic model with three sources of uncertainty. An analytical approximate expression for the PDF of the lifetime for the degradation model has been derived to provide decision makers with a reliable tool for making planned maintenance decisions. The unknown parameters of the model are effectively estimated using the EM algorithm. Comparing our results with those obtained using existing approaches, we can see that our model with three sources of uncertainty improves the accuracy of the PDF of lifetime and that the MLE obtained by the EM algorithm increases the accuracy of the parameter estimation. A practical case study for fatigue crack propagation data is analyzed and discussed, showing that our model provides more accurate estimates and has great potential advantages in preventive maintenance scheduling [44].



---

*Source: 6475056-2021-12-29.xml*