Remaining useful life prediction of implicit linear Wiener degradation process based on multi-source information

Abstract

Accurate remaining useful life (RUL) prediction is helpful to improve the reliability and safety of complex systems. However, in practical engineering applications, it often occurs imperfect or scarce prior degradation information for the degradation system with measurement error (ME). In order to solve this problem, based on the implicit linear Wiener degradation process, a RUL prediction method which reasonably fuses failure time data or multi-source information is proposed in this paper. Firstly, based on the implicit linear Wiener degradation process, we obtain the relationship between the natures of parameters estimation and degradation data by theoretical derivation, which provides a theoretical basis regarding how to fuse multi-source information. Secondly, according to the natures of parameters estimation, we use field degradation data and historical degradation data to estimate the fixed parameters of the two prediction cases respectively, and fuse failure time data into the degradation model by the expectation maximization (EM) algorithm. Then, the Kalman filtering algorithm is used to online update the drift parameter based on field degradation data. Finally, we use some simulation experiments to further verify the natures of parameters estimation, and two practical case studies to verify the superiority of the proposed method.

Introduction

In practical engineering applications, the performance of the equipment changes cumulatively over time, which may lead to equipment shutdown and even accidents when its degradation level exceeds the failure threshold. ¹ Accurate remaining useful life (RUL) prediction helps operators make correct maintenance decisions before the equipment approaches failure, so as to effectively improve the reliability and safety and reduce the economic cost of the equipment. ² Here, RUL is defined as the length of time from the present time to the failure time of the equipment. ³ In recent years, with the rapid development of aerospace, military, and other high reliability fields, statistical data-driven RUL prediction method has become a popular research issue in the field of reliability.^4–6

Statistical data-driven RUL prediction method can make full use of the degradation data of congeneric equipment. ⁷ Its main principle is to estimate the unknown parameters representing the common characteristics of degradation processes based on historical degradation data or failure time data, then online update the random parameters representing the unit-to-unit variability among equipment based on the field degradation data detected in real time, and finally predict the RUL by the probability density function (PDF) of RUL. Among them, the failure time data reveals the reliability information on the time scale, the historical degradation data describes the degradation characteristics of the failure process, ⁸ and the field degradation data represents the individual characteristic of the assessed equipment. Although these three kinds of data contain different information, they can all be used for RUL prediction. Due to the randomness of the degradation process, stochastic processes are often used to predict RUL. ⁹ RUL prediction methods based on stochastic processes have the advantages of clear physical explanations and great mathematical properties. Gebraeel et al. ¹⁰ are the first authors using Bayesian method to reasonably fuse field degradation information and prior degradation information. This method is first used for random coefficient regression model,^10–12 and then applied to Gamma process, ¹³ two-phase Gamma process, ¹⁴ inverse Gaussian process, ¹⁵ inverse Gaussian process with random effects, ¹⁶ Wiener degradation process, ¹⁷ two-phase Wiener process, ¹⁸ Wiener process with skew-normal random effects, ¹⁹ adaptive Wiener process, ²⁰ etc. Since Wiener process model is suitable for monotonic and nonmonotonic degradation systems simultaneously, it has been widely used for many equipment in recent years.^21,22 However, in practical engineering applications, it often occurs imperfect or scarce prior degradation information, which may lead to inaccurate or impossible estimation of prior parameters,^11,23 and further lead to accuracy reduction or even failure of RUL prediction.^9,17 In the existing literature, there are mainly two ways to solve this problem.

The first way is to combine expectation maximization (EM) algorithm and Bayesian method, which was first proposed by Wang et al. ²⁴ to fit the adaptive Brownian motion model with drift parameter. Then, Si et al. ¹⁷ used Bayesian updating method and EM algorithm to update the parameters and the RUL distribution, and deduced an accurate and closed RUL distribution function based on the linear Wiener degradation process. Later, this method has been applied to implicit linear Wiener degradation process, ²⁵ nonlinear Wiener degradation process,^26,27 and implicit nonlinear Wiener degradation process.^20,28,29 The combination of EM algorithm and Bayesian method can overcome the influence of imperfect or scarce prior information on RUL prediction. However, the reason why this method can obtain better prediction results than the traditional Bayesian method is still unclear. By assuming that the drift parameter is fixed, Tang et al. ³⁰ deduced the analytical expression of parameters based on the implicit linear Wiener degradation process, and proved that the combined algorithm has the same estimation results as the traditional maximum likelihood estimation (MLE) method. It indicates that EM algorithm completely gets rid of the influence of prior information and only depends on field degradation data. Based on this conclusion, Tang et al. ³⁰ proposed a heuristic algorithm to reasonably fuse prior information and field degradation data for linear degradation model. Then, Changhao et al. ²³ generalized it to the nonlinear Wiener degradation process, and Han et al. ³¹ further generalized it to the implicit nonlinear Wiener degradation process.

The second way is to fuse multi-source information. Generally, there are three kinds of data used for parameters estimation, that is, historical degradation data of congeneric equipment (including accelerated degradation data ³² ), failure time data of congeneric equipment, and field degradation data of the assessed equipment. How to use the degradation information contained in these three kinds of data to predict the RUL is called the RUL prediction method with fusing multi-source information. Tang et al. ³³ obtained the relationship between the parameters estimation and the characteristics of degradation data based on the basic linear Wiener degradation process, which provides a theoretical basis regarding how to fuse multi-source information. However, in practical engineering applications, there may be scarce failure time data or historical degradation data especially for new equipment. Therefore, it is necessary to reasonably fuse these two kinds of data to estimate the unknown prior parameters. At present, some scholars have used the method with fusing multi-source information to predict RUL.^34–37 For example, due to the different stress between accelerated degradation test data and field degradation data, Wang et al. ³⁴ proposed a Bayesian evaluation method based on joint likelihood function to fuse accelerated degradation test data and failure time data. Pang et al. ³⁵ further studied the method with fusing accelerated degradation test data and failure time data based on nonlinear Wiener degradation process. Liu et al. ³⁶ fused a small amount of failure time data and degradation data to predict the RUL of equipment under the Bayesian framework. Zhang et al. ³⁷ proposed a new Bayesian framework to fuse binary degradation data and failure time data for reliability analysis. However, these literatures did not consider the influence of random effects. Tang et al. ³⁸ proved that the RUL prediction with considering random effects could obtain better results in theory. However, Tang et al. ³⁸ only used the historical degradation data of congeneric equipment when estimating unknown parameters, without fusing the failure time data. Recently, Tang et al. ³³ proposed a RUL prediction method based on the basic linear Wiener degradation process and the nonlinear Wiener degradation process respectively, which fuses failure time data and considering random effects, and verified the superiority of this method by comparing it with the traditional MLE method. Wang et al. ³⁹ further proposed a RUL prediction method based on the basic linear Wiener degradation process and the nonlinear Wiener degradation process respectively, which fuses multi-source information and considering random effects, and verified the superiority of the method by comparing it with the RUL prediction method using only degradation data or failure time data. However, due to the measuring instrument, the technical level of the operators and environmental factors, the measured degradation data may have measurement error (ME), which could not directly reflect the actual degradation state. ⁴⁰ Considering that the RUL prediction results were influenced by three sources of variability, that is, temporal variability, unit-to-unit variability, and measurement variability, Si et al. ⁴¹ proposed an online updating method of random parameters based on Kalman filter algorithm, and obtained the analytical expression of PDF of RUL. For the degradation equipment with ME, Tang et al. ⁴² proved that if the influence of ME is not considered, which may lead to premature maintenance, thus increase the maintenance cost. Therefore, the influence of ME needs to be considered when predicting RUL. However, in the existing literature, the method with fusing multi-source information is only used in the basic Wiener degradation process, and has not been generalized to the implicit Wiener degradation process.

From the above review over related works, the RUL prediction method under imperfect or scarce prior information has not been studied thoroughly. Therefore, this paper mainly studies the RUL prediction method based on the implicit linear Wiener degradation process. Firstly, we obtain some natures regarding parameters estimation of the implicit linear Wiener degradation process model, which has not been reported before. Secondly, according to these natures of parameters estimation, we propose a RUL prediction method with fusing multi-source information based on the implicit linear Wiener degradation process model. When the prior degradation information is imperfect, this method uses field degradation data to estimate the fixed parameters, and then uses failure time data to estimate the prior drift parameter. When the prior degradation information is scarce, this method estimates the prior parameters by fusing failure time data and historical degradation data. In the above two cases, the drift parameter and the actual degradation state was online updated based on field degradation data. Finally, we use some simulation experiments and two practical case studies to verify the effectiveness of the proposed method.

The remainders of this paper are presented as follow. Section “Natures of parameters estimation for the implicit linear Wiener degradation process” analyzes the natures of parameters estimation based on the implicit linear Wiener degradation process. Section “RUL prediction of the implicit linear Wiener degradation process” introduces the RUL prediction method of the implicit linear Wiener degradation process. In order to solve the problem of imperfect or scarce prior degradation information, a RUL prediction method with fusing multi-source information based on the implicit linear Wiener degradation is proposed in Section “Prior parameters estimation of the implicit linear Wiener degradation process based on multi-source information”. Section “Experimental studies” verifies the effectiveness of the proposed RUL prediction method by some simulation experiments and two practical case studies.

Natures of parameters estimation for the implicit linear Wiener degradation process

The prior parameters describe the common degradation characteristics of equipment. Accurate estimation of the prior parameters is helpful to improve the accuracy of RUL prediction. In practical engineering applications, the observed degradation data often have ME. If the influence of ME is not considered, the prediction accuracy of RUL could be reduced, and it could even obtain wrong estimation results. ⁴² Therefore, this section analyzes the natures of parameters estimation based on the implicit linear Wiener degradation process.

Let $Y(t)$ denote the measured degradation data of the unit at time $t$ .Then, it can be expressed as follows

$$Y(t)=x_0+\lambda t+{\sigma }_{B}B(t)+\epsilon$$

(1)

where $x_0$ is the initial degradation state, $\lambda$ is the drift parameter denoting the degradation rate, ${\sigma }_B$ is the diffusion parameter, $B(t)$ is the standard Brownian motion, and $\epsilon$ is the ME. Without loss of generality, we set $x_0=0$ . In order to characterize individual differences among equipment, $\lambda$ is assumed to be random and follow normal distribution, that is, $\lambda ~N({\mu }_{\lambda },{\sigma }_{\lambda }^2)$ . Moreover, $\lambda$ is assumed to be s-independent with $B(t)$ and $\epsilon$ . Based on equation (1), the parameters to be estimated are $\Theta =\{{\mu }_{\lambda },{\sigma }_{\lambda }^2,{\sigma }_B^2,{\sigma }_{\epsilon }^2\}$ .

Suppose that there are $n$ units for testing and each unit are measured at the same times $t_1,t_2,···,t_m$ , then the measured degradation state of the ith unit at time $t_j$ can be expressed as

$$y_{i,j}=x_0+{\lambda }_i{t}_j+{\sigma }_{B}B(t_j)+{\epsilon }_{i,j}$$

(2)

Let $\mathrm{\Delta }{y}_{i,j}=y_i(t_j)-y_i(t_{j-1})$ , where $1\le i\le n$ , $1\le j\le m$ , and ${\lambda }_i$ is the degradation rate of the ith unit, then $\mathrm{\Delta }{y}_{i,j}$ can be expressed as

$$\mathrm{\Delta }{y}_{i,j}=y_i(t_j)-y_i(t_{j-1})={\lambda }_i\mathrm{\Delta }{t}_j+{\sigma }_{B}B(\mathrm{\Delta }{t}_j)+{\epsilon }_{i,j}-{\epsilon }_{i,j-1}$$

(3)

where, $\mathrm{\Delta }{t}_j=t_j-t_{j-1}$ , ${\epsilon }_{i,j}$ is the ME of the ith unit at time $t_j$ . Let $\mathrm{\Delta }{\mathit{y}}_i=(\mathrm{\Delta }{y}_{i,1},\mathrm{\Delta }{y}_{i,2},···,\mathrm{\Delta }{y}_{i,m}{)}^{\prime }$ , $\mathrm{\Delta }\mathit{t}=(\mathrm{\Delta }{t}_1,\mathrm{\Delta }{t}_2,···,\mathrm{\Delta }{t}_m{)}^{\prime }$ , and $\mathit{Y}=(\mathrm{\Delta }{\mathit{y}}_1^{\prime },\mathrm{\Delta }{\mathit{y}}_2^{\prime },···,\mathrm{\Delta }{\mathit{y}}_n^{\prime }{)}^{\prime }$ . Then, $\mathrm{\Delta }{\mathit{y}}_i$ follows multivariate normal distribution with mean ${\mu }_{\lambda }\mathrm{\Delta }\mathit{t}$ , and variance $\mathbf{\Sigma }\mathbf{=}{\mathbf{\sigma }}_{\mathbf{\lambda }}^2\mathbf{\Delta }\mathit{t}\mathbf{\Delta }{\mathit{t}}^{\mathbf{’}}\mathbf{+}{\mathbf{\sigma }}_{\mathbf{B}}^2\mathbf{\Omega }\mathbf{+}{\mathbf{\sigma }}_{\mathbf{\epsilon }}^2\mathit{F}$ , where $\mathbf{\Omega }\mathbf{=}\mathrm{diag}\mathbf{(}\mathbf{\Delta }{\mathbf{t}}_{\mathbf{j}}\mathbf{)}$ and the elements of $\mathit{F}$ are shown in equation (4). Generally, since the initial state of the unit can be considered as a fixed value and there is no ME, it is set that $f_{1,1}=1$ . ⁴³

$$f_{i,j}=f_{j,i}=\{\begin{array}{cc}0\hfill & i>j+1\hfill \\ \multicolumn{1}{c}{-1}& i=j+1\hfill \\ \multicolumn{1}{c}{1}& i=j=1\hfill \\ \multicolumn{1}{c}{2}& i=j>1\hfill \end{array}$$

(4)

Let $\mathit{D}=\mathbf{\Omega }\mathbf{+}\mathbf{\varphi }\mathit{F}$ , where $\varphi ={\sigma }_{\epsilon }^2/{\sigma }_B^2$ , then $\mathrm{\Delta }{\mathit{y}}_i~N({\mu }_{\lambda }\mathrm{\Delta }\mathit{t},{\sigma }_{\lambda }^2\mathrm{\Delta }\mathit{t}\mathrm{\Delta }{\mathit{t}}^{\prime }+{\sigma }_B^2\mathit{D})$ . And the expression of unbiased parameters estimation regarding $\mathbf{\Theta }\mathbf{=}\mathbf{\{}{\mathbf{\mu }}_{\mathbf{\lambda }}\mathbf{,}{\mathbf{\sigma }}_{\mathbf{\lambda }}^2\mathbf{,}{\mathbf{\sigma }}_{\mathbf{B}}^2\mathbf{,}{\mathbf{\sigma }}_{\mathbf{\epsilon }}^2\mathbf{\}}$ for the implicit linear Wiener degradation process can be written as ³⁸

$${\hat{\mu }}_{\lambda }(\hat{\varphi })=\frac{1}{n\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}}\sum _{i=1}^n\mathrm{\Delta }{\mathit{y}}_i^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}$$

(5)

$${\widehat{\sigma }}_{\lambda }^2(\widehat{\varphi })=\frac{1}{n-1}\sum _{i=1}^n{\left({\mu }_{\lambda }-{\widehat{\lambda }}_i\right)}^2-\frac{1}{n(m-1)\Delta t^{\prime }{D}^{-1}\Delta t}\sum _{i=1}^n{\left(\Delta y_i-{\widehat{\lambda }}_i\Delta t\right)}^{\prime }{D}^{-1}\left(\Delta y_i-{\widehat{\lambda }}_i\Delta t\right)$$

(6)

$${\hat{\sigma }}_B^2(\hat{\varphi })=\frac{1}{n(m-1)}\sum _{i=1}^n{(\mathrm{\Delta }{\mathit{y}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})}^{\prime }{\mathit{D}}^{-1}(\mathrm{\Delta }{\mathit{y}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})$$

(7)

where ${\hat{\lambda }}_i=\mathrm{\Delta }{\mathit{y}}_i^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}/(\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})$ and the specific derivation process can be found in Tang et al. ³⁸

The estimation of $\varphi$ can be obtained by maximizing equation (8) by the MATLAB function “FMINSEARCH.” Then, ${\mu }_{\lambda }$ , ${\sigma }_{\lambda }^2$ , and ${\sigma }_B^2$ can be obtained by substituting $\varphi$ into equations (5) to (7), and further obtains ${\sigma }_{\epsilon }^2$ .

$$\begin{array}{cc}\multicolumn{1}{c}{\ln L(\varphi |\mathbf{Y})}& =-\frac{\mathrm{mn}}{2}\ln 2\pi -\frac{n(m-1)}{2}\hfill \\ \multicolumn{1}{c}{}& \ln \left(\frac{1}{n(m-1)}\sum _{i=1}^n{(\mathrm{\Delta }{\mathit{y}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})}^{\prime }{\mathit{D}}^{-1}(\mathrm{\Delta }{\mathit{y}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})\right)\hfill \\ \multicolumn{1}{c}{}& -\frac{n}{2}\ln \left[\frac{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}}{n}\sum _{i=1}^n{({\hat{\mu }}_{\lambda }-\frac{\mathrm{\Delta }{\mathit{y}}_i^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}}{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}})}^2\right]\hfill \\ \multicolumn{1}{c}{}& -\frac{n}{2}\ln |\mathit{D}|-\frac{\mathrm{mn}}{2}\hfill \end{array}$$

(8)

At present, some natures of the unbiased parameters estimation have been proved by scholars, as shown in Lemma 1. ³⁸ Since it is very complicated to obtain the expression of all parameters at the same time, here it is assumed that ${\hat{\sigma }}_{\epsilon }^2$ is the actual ${\sigma }_{\epsilon }^2$ .

• Lemma 1 ³⁸ : the expectation of unbiased estimation parameters ${\hat{\mu }}_{\lambda }$ , ${\hat{\sigma }}_{\lambda }^2$ , and ${\hat{\sigma }}_B^2$ of the implicit linear Wiener degradation process are

$$\begin{array}{cc}\multicolumn{1}{c}{}& E({\hat{\mu }}_{\lambda })={\mu }_{\lambda },\hfill \\ \multicolumn{1}{c}{}& E({\hat{\sigma }}_B^2)={\sigma }_B^2,\hfill \\ \multicolumn{1}{c}{}& E({\hat{\sigma }}_{\lambda }^2)={\sigma }_{\lambda }^2.\hfill \end{array}$$

(9)

More details about the proof process can be found in Tang et al. ³⁸ Lemma 1 gives the expectation of unbiased parameters estimation. Based on Lemma 1, we can know that the results of unbiased estimation parameters can converge to the actual value, which reflects the effectiveness of unbiased parameter estimation method. However, which factors are related to the accuracy of parameters estimation? It cannot be obtained from Lemma 1. Therefore, we further deduce the variance of unbiased parameters estimation, and the specific results are shown in Theorem 1.

• Theorem 1: The variance of unbiased parameters estimation ${\hat{\mu }}_{\lambda }$ , ${\hat{\sigma }}_{\lambda }^2$ , and ${\hat{\sigma }}_B^2$ of the implicit linear Wiener degradation process are

$$D({\hat{\mu }}_{\lambda })=\frac{1}{n}({\sigma }_{\lambda }^2+\frac{{\sigma }_B^2}{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}})$$

(10)

$$D\left({\hat{\sigma }}_B^2\right)=\frac{2{({\sigma }_B^2)}^2}{n(m-1)}$$

(11)

$$D\left({\hat{\sigma }}_{\lambda 1}^2\right)=\frac{2}{n-1}{({\sigma }_{\lambda }^2+\frac{{\sigma }_B^2}{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}})}^2$$

(12)

$$D\left({\hat{\sigma }}_{\lambda 2}^2\right)=\frac{2{({\sigma }_B^2)}^2}{{(\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})}^{2}n(m-1)}$$

(13)

where ${\hat{\sigma }}_{\lambda 1}^2=\frac{1}{n-1}\sum _{i=1}^n{({\hat{\mu }}_{\lambda }-{\hat{\lambda }}_i)}^2$ and ${\hat{\sigma }}_{\lambda 2}^2=\frac{1}{n(m-1)\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}}\sum _{i=1}^n{(\mathrm{\Delta }{\mathit{z}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})}^{\prime }{\mathit{D}}^{-1}(\mathrm{\Delta }{\mathit{z}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})$ .

Since it is very complicated to directly calculate the variance of ${\hat{\sigma }}_{\lambda }^2$ , we split the expression of ${\hat{\sigma }}_{\lambda }^2$ into ${\hat{\sigma }}_{\lambda 1}^2$ and ${\hat{\sigma }}_{\lambda 2}^2$ to calculate the variance separately. It can be seen that ${\hat{\sigma }}_{\lambda 1}^2$ and ${\hat{\sigma }}_{\lambda 2}^2$ have the same characteristics, which are inversely proportional to the number of units and the length of detection time of a single unit.

Proof: According to equation (1) and the natures of Wiener degradation process, we can obtain $\mathrm{\Delta }{\mathit{y}}_i~\mathrm{MVN}({\mu }_{\lambda }\mathrm{\Delta }\mathit{t},{\sigma }_{\lambda }^2\mathrm{\Delta }\mathit{t}\mathrm{\Delta }{\mathit{t}}^{\prime }+{\sigma }_B^2\mathit{D})$ and $\mathrm{\Delta }{\mathit{y}}_i~N({\lambda }_i\mathrm{\Delta }\mathit{t},{\sigma }_B^2\mathit{D})$ . On account of $\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }{\mathit{y}}_i=\mathrm{\Delta }{\mathit{y}}_i^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}$ , we can obtain $\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }{\mathit{y}}_i~N({\mu }_{\lambda }\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t},{\sigma }_{\lambda }^2{(\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})}^2+{\sigma }_B^2\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})$ , $\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }{\mathit{y}}_i/(\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})~N({\mu }_{\lambda },{\sigma }_{\lambda }^2+{\sigma }_B^2/(\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}))$ , and $\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }{\mathit{y}}_i/(\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})~N({\lambda }_i,{\sigma }_B^2/(\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}))$ . Then ${\hat{\mu }}_{\lambda }$ follows normal distribution, which can be written as follows.

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\hat{\mu }}_{\lambda }=\frac{1}{n\mathrm{\Delta }\mathit{t}{\mathit{D}}^{-1}\mathrm{\Delta }{\mathit{t}}^{\prime }}\sum _{i=1}^n\mathrm{\Delta }\mathit{t}{\mathit{D}}^{-1}\mathrm{\Delta }{\mathit{y}}_i^{\prime }~\hfill \\ \multicolumn{1}{c}{}& N({\mu }_{\lambda },\frac{1}{n}({\sigma }_{\lambda }^2+\frac{{\sigma }_B^2}{\mathrm{\Delta }\mathit{t}{\mathit{D}}^{-1}\mathrm{\Delta }{\mathit{t}}^{\prime }}))\hfill \end{array}$$

(14)

According to equation (14), we have

$$D({\hat{\mu }}_{\lambda })=\frac{1}{n}({\sigma }_{\lambda }^2+\frac{{\sigma }_B^2}{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}})$$

(15)

In order to calculate the expectation and the variance of ${\hat{\sigma }}_B^2$ , $\mathrm{\Delta }{y}_i=(\mathrm{\Delta }{y}_{i,1},\mathrm{\Delta }{y}_{i,2},···,\mathrm{\Delta }{y}_{i,m}{)}^{\prime }$ , and $\mathrm{\Delta }\mathit{t}=(\mathrm{\Delta }{t}_1,\mathrm{\Delta }{t}_2,···,\mathrm{\Delta }{t}_m{)}^{\prime }$ are defined. According to the natures of Wiener process, we can obtain $\mathrm{\Delta }{y}_i~\mathrm{MVN}({\lambda }_i\mathrm{\Delta }\mathit{t},{\sigma }_B^2\mathit{D})$ , where ${\sigma }_B^2\mathit{D}={\sigma }_B^2\mathbf{\Omega }\mathbf{+}{\mathbf{\sigma }}_{\mathbf{\epsilon }}^2\mathit{F}$ . Then, we have

$$\frac{1}{{\sigma }_B^2}{(\mathrm{\Delta }{\mathit{z}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})}^{\prime }{\mathit{D}}^{-1}(\mathrm{\Delta }{\mathit{z}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})~{\chi }^2(m-1)$$

(16)

According to equation (16) and the natures of chi square distribution, we can obtain the expectation and variance of ${\hat{\sigma }}_B^2$ as follows.

$$\begin{array}{cc}\multicolumn{1}{c}{D\left({\hat{\sigma }}_B^2\right)}& =D\left(\frac{1}{n(m-1)}\sum _{i=1}^n{(\mathrm{\Delta }{\mathit{y}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})}^{\prime }{\mathit{D}}^{-1}(\mathrm{\Delta }{\mathit{y}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})\right)\hfill \\ \multicolumn{1}{c}{}& =\frac{1}{n^2{(m-1)}^2}\sum _{i=1}^{n}D\left[{\sigma }_B^2\frac{1}{{\sigma }_B^2}{(\mathrm{\Delta }{\mathit{y}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})}^{\prime }{\mathit{D}}^{-1}(\mathrm{\Delta }{\mathit{y}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})\right]\hfill \\ \multicolumn{1}{c}{}& =\frac{1}{n^2{(m-1)}^2}\sum _{i=1}^n{({\sigma }_B^2)}^{2}2(m-1)\hfill \\ \multicolumn{1}{c}{}& =\frac{2{({\sigma }_B^2)}^2}{n(m-1)}\hfill \end{array}$$

(17)

Additionally, since ${\hat{\lambda }}_i=\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }{\mathit{y}}_i/(\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})~N({\mu }_{\lambda },{\sigma }_{\lambda }^2+{\sigma }_B^2/\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})$ , $\sum _{i=1}^n{({\hat{\mu }}_{\lambda }-{\hat{\lambda }}_i)}^2/(({\sigma }_{\lambda }^2+{\sigma }_B^2/\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}))~{\chi }^2(n-1)$ can be obtained. Then, we can obtain the expectation and variance of ${\hat{\sigma }}_{\lambda 1}^2$ as follows.

$$\begin{array}{cc}\multicolumn{1}{c}{E({\hat{\sigma }}_{\lambda 1}^2)}& =E\left(\frac{1}{n-1}\sum _{i=1}^n{({\hat{\mu }}_{\lambda }-{\hat{\lambda }}_i)}^2\right)\hfill \\ \multicolumn{1}{c}{}& =\frac{1}{n-1}({\sigma }_{\lambda }^2+\frac{{\sigma }_B^2}{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}})\hfill \\ \multicolumn{1}{c}{}& E\left(\frac{1}{({\sigma }_{\lambda }^2+{\sigma }_B^2/\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})}\sum _{i=1}^n{({\hat{\mu }}_{\lambda }-{\hat{\lambda }}_i)}^2\right)\hfill \\ \multicolumn{1}{c}{}& =\frac{1}{n-1}({\sigma }_{\lambda }^2+\frac{{\sigma }_B^2}{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}})(n-1)={\sigma }_{\lambda }^2+\frac{{\sigma }_B^2}{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}}\hfill \end{array}$$

(18)

$$\begin{array}{cc}\multicolumn{1}{c}{D({\hat{\sigma }}_{\lambda 1}^2)}& =D\left(\frac{1}{n-1}\sum _{i=1}^n{({\hat{\mu }}_{\lambda }-{\hat{\lambda }}_i)}^2\right)\hfill \\ \multicolumn{1}{c}{}& ={\left(\frac{1}{n-1}\right)}^2({\sigma }_{\lambda }^2+\frac{{\sigma }_B^2}{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}})\hfill \\ \multicolumn{1}{c}{}& \left(\frac{1}{({\sigma }_{\lambda }^2+{\sigma }_B^2/\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})}\sum _{i=1}^n{({\hat{\mu }}_{\lambda }-{\hat{\lambda }}_i)}^2\right)\hfill \\ \multicolumn{1}{c}{}& ={\left(\frac{1}{n-1}\right)}^2({\sigma }_{\lambda }^2+\frac{{\sigma }_B^2}{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}})2(n-1)\hfill \\ \multicolumn{1}{c}{}& =\left(\frac{2}{n-1}\right)({\sigma }_{\lambda }^2+\frac{{\sigma }_B^2}{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}})\hfill \end{array}$$

(19)

Based on equation (16), we can obtain the expectation and variance of ${\hat{\sigma }}_{\lambda 2}^2$ as follows.

$$\begin{array}{l}E({\widehat{\sigma }}_{\lambda 2}^2)=E\left(\frac{1}{n(m-1)\Delta t^{\prime }{D}^{-1}\Delta t}\sum _{i=1}^n{\left(\Delta z_i-{\widehat{\lambda }}_i\Delta t\right)}^{\prime }{D}^{-1}\left(\Delta z_i-{\widehat{\lambda }}_i\Delta t\right)\right)\\ =\frac{{\sigma }_B^2}{n(m-1)\Delta t^{\prime }{D}^{-1}\Delta t}\sum _{i=1}^{n}E\left(\frac{1}{{\sigma }_B^2}{\left(\Delta z_i-{\widehat{\lambda }}_i\Delta t\right)}^{\prime }{D}^{-1}\left(\Delta z_i-{\widehat{\lambda }}_i\Delta t\right)\right)\\ =\frac{{\sigma }_B^2}{n(m-1)\Delta t^{\prime }{D}^{-1}\Delta t}n(m-1)=\frac{{\sigma }_B^2}{\Delta t^{\prime }{D}^{-1}\Delta t}\end{array}$$

(20)

$$\begin{array}{cc}\multicolumn{1}{c}{D({\hat{\sigma }}_{\lambda 2}^2)}& =D(\frac{1}{n(m-1)\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}}\hfill \\ \multicolumn{1}{c}{}& \sum _{i=1}^n{(\mathrm{\Delta }{\mathit{z}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})}^{\prime }{\mathit{D}}^{-1}(\mathrm{\Delta }{\mathit{z}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t}))\hfill \\ \multicolumn{1}{c}{}& =\frac{{({\sigma }_B^2)}^2}{{(n(m-1)\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})}^2}\hfill \\ \multicolumn{1}{c}{}& \sum _{i=1}^{n}D\left(\frac{1}{{\sigma }_B^2}{(\mathrm{\Delta }{\mathit{z}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})}^{\prime }{\mathit{D}}^{-1}(\mathrm{\Delta }{\mathit{z}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})\right)\hfill \\ \multicolumn{1}{c}{}& =\frac{{({\sigma }_B^2)}^2}{{(n(m-1)\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})}^2}2n(m-1)\hfill \\ \multicolumn{1}{c}{}& =\frac{2{({\sigma }_B^2)}^2}{n(m-1){(\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})}^2}\hfill \end{array}$$

(21)

This completes the proof.

According to Theorem 1, we have obtained the variance of unbiased parameters estimation ${\hat{\mu }}_{\lambda }$ , ${\hat{\sigma }}_{\lambda }^2$ , and ${\hat{\sigma }}_B^2$ of the implicit linear Wiener degradation process. In addition, ${\hat{\sigma }}_{\epsilon }^2$ also needs to be analyzed. Similar to the above derivation process, let $\mathit{P}=\phi \mathbf{\Omega }\mathbf{+}\mathit{F}$ , where $\phi ={\sigma }_B^2/{\sigma }_{\epsilon }^2$ , then $\mathrm{\Delta }{\mathit{y}}_i~N({\mu }_{\lambda }\mathrm{\Delta }\mathit{t},{\sigma }_{\lambda }^2\mathrm{\Delta }\mathit{t}\mathrm{\Delta }{\mathit{t}}^{\prime }+{\sigma }_{\epsilon }^2\mathit{P})$ , and then the expression of ${\hat{\sigma }}_{\epsilon }^2$ can be expressed as

$${\hat{\sigma }}_{\epsilon }^2(\hat{\phi })=\frac{1}{n(m-1)}\sum _{i=1}^n{(\mathrm{\Delta }{\mathit{y}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})}^{\prime }{\mathit{P}}^{-1}(\mathrm{\Delta }{\mathit{y}}_i-{\hat{\lambda }}_i\mathrm{\Delta }\mathit{t})$$

(22)

According to equation (22), we can observe that the expressions of ${\hat{\sigma }}_B^2$ and ${\hat{\sigma }}_{\epsilon }^2$ have the same structure and are proportional to the number of units and the detection times of a single unit. Therefore, the expectation and variance of ${\hat{\sigma }}_B^2$ and ${\hat{\sigma }}_{\epsilon }^2$ have the same derivation process, and the results are given in Theorem 2.

• Theorem 2: The expectation and variance of unbiased parameter estimation ${\hat{\sigma }}_{\epsilon }^2$ of the implicit linear Wiener degradation process can be written as follows.

$$E({\hat{\sigma }}_{\epsilon }^2)={\sigma }_{\epsilon }^2$$

(23)

$$D({\hat{\sigma }}_{\epsilon }^2)=\frac{2{({\sigma }_{\epsilon }^2)}^2}{n(m-1)}$$

(24)

Similarly, ${\hat{\sigma }}_B^2$ is assumed as the actual value. The specific derivation process is omitted here.

• Remark 1: Based on above two theorems, the following conclusions can be obtained.

• (1) From equations (10), (12) and (13), we can observe that the estimation accuracy of ${\mu }_{\lambda }$ and ${\sigma }_{\lambda }^2$ is mainly affected by the number of units and the length of detection time, that is, the more the number of units $n$ or the longer the length of time $t_m$ , the closer the parameters estimation results are to the actual value. Therefore, when the historical degradation data is scarce or even absent, fusing failure time data into the degradation model can improve the estimation accuracy of ${\mu }_{\lambda }$ and ${\sigma }_{\lambda }^2$ from a certain extent.

• (2) According to equations (11) and (24), the estimation accuracy of ${\sigma }_B^2$ and ${\sigma }_{\epsilon }^2$ denoting the degradation uncertainty is directly proportional to the number of units $n$ and the detection times $m-1$ . In other words, when estimating the unknown parameters, the more degradation data are used, the closer the estimated ${\sigma }_B^2$ and ${\sigma }_{\epsilon }^2$ are to the actual value.

RUL prediction of the implicit linear Wiener degradation process

Online updating of random parameters based on Kalman filter algorithm

As mentioned above, the randomness of the drift parameter describes the unit-to-unit variability. In order to adapt to the individual degradation characteristics of the assessed unit, it is necessary to online update the drift parameter using the field degradation data. Since the detected degradation data have ME, it cannot directly reflect the actual degradation state. Therefore, the Kalman filter algorithm can be used to online update the drift parameter and estimate the current degradation state in real time. ⁴¹

For the implicit linear Wiener degradation process as shown in equation (1), when the degradation state ${\mathit{Y}}_{1:k}$ of the assessed unit at time $t_k$ is detected, the state space equation at time $t_k$ can be expressed as

$$\{\begin{array}{c}{x}_k=x_{k-1}+{\lambda }_{k-1}(t_k-t_{k-1})+v_k\hfill \\ \multicolumn{1}{c}{{\lambda }_k={\lambda }_{k-1}}\\ \multicolumn{1}{c}{y_k=x_k+{\epsilon }_k}\end{array}$$

(25)

where $v_k={\sigma }_{B}B(t_k)-{\sigma }_{B}B(t_{k-1})$ , ${\lambda }_k$ is drift parameter, $x_k$ denotes the actual degradation state at time $t_k$ , and ${\epsilon }_k$ denotes the ME at time $t_k$ . Here, ${\left\{v_k\right\}}_{k\ge 1}$ and ${\left\{{\epsilon }_k\right\}}_{k\ge 1}$ are independent and identically distributed.

Due to the influence of ME, the actual degradation state and drift parameter are hidden state. Hence, the state space equation (25) can be further expressed as ⁴¹

$$\{\begin{array}{c}{\mathit{z}}_k={\mathit{A}}_k{\mathit{z}}_{k-1}+{\mathbf{\eta }}_{\mathbf{k}}\hfill \\ \multicolumn{1}{c}{y_k=\mathit{C}{\mathit{z}}_k+{\epsilon }_k}\end{array}$$

(26)

According to Si et al., ⁴¹ the state equation can be updated by Kalman filter algorithm. Let ${\hat{x}}_{k|k}=E\left(x_k|{\mathit{Y}}_{1:k}\right)$ , ${\hat{\lambda }}_{k|k}=E\left({\lambda }_k|{\mathit{Y}}_{1:k}\right)$ , ${\kappa }_{x,k}^2=\mathrm{var}\left(x_k|{\mathit{Y}}_{1:k}\right)$ , ${\kappa }_{\lambda ,k}^2=\mathrm{var}\left({\lambda }_k|{\mathit{Y}}_{1:k}\right)$ , ${\kappa }_{x\lambda ,k}^2=\mathrm{cov}\left(x_k{\lambda }_k|{\mathit{Y}}_{1:k}\right)$ . Then the online update results of the drift parameter and the estimated current degradation state can be expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\lambda }_k|{\mathit{Y}}_{1:k}~N({\hat{\lambda }}_{k|k},{\kappa }_{\lambda ,k}^2),\hfill \\ \multicolumn{1}{c}{}& x_k|{\mathit{Y}}_{1:k}~N({\hat{x}}_{k|k},{\kappa }_{x,k}^2),\hfill \\ \multicolumn{1}{c}{}& x_k|{\lambda }_k,{\mathit{Y}}_{1:k}~N({\mu }_{x_k|\lambda ,k},{\sigma }_{x_k|\lambda ,k}^2).\hfill \end{array}$$

(27)

where

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\mu }_{x_k|\lambda ,k}={\hat{x}}_{k|k}+{\rho }_k\frac{{\kappa }_{x,k}}{{\kappa }_{\lambda ,k}}({\lambda }_k-{\hat{\lambda }}_{k|k}),\hfill \\ \multicolumn{1}{c}{}& {\sigma }_{x_k|\lambda ,k}^2={\kappa }_{x,k}^2(1-{\rho }_k^2),\hfill \\ \multicolumn{1}{c}{}& {\rho }_k={\kappa }_{x\lambda ,k}^2/({\kappa }_{\lambda ,k}{\kappa }_{x,k}).\hfill \end{array}$$

(28)

RUL prediction

RUL refers to the length of time from the present time to the failure time. Define the failure threshold as $w$ , the observed data of the assessed unit at time $t_k$ as ${\mathit{Y}}_{1:k}=\{y_1,y_2,···,y_k\}$ , and the actual degradation state as ${\mathit{X}}_{1:k}=\{x_1,x_2,···,x_k\}$ , then, according to the concept of first hitting time (FHT),^33,44 once the actual degradation state ${\mathit{X}}_{1:k}$ at time $t_k$ is obtained, RUL can be expressed as

$$L_k=inf\{l_k:X(l_k+t_k)\ge w|X_{1:k}\}$$

(29)

Based on equation (1), when the actual degradation state $x_k$ is obtained, the PDF of the RUL can be expressed as ⁴¹

$$f_{L_k}(l_k)=\frac{w-x_k}{\sqrt{2\pi {l_k}^3({\sigma }_{\lambda ,k}^2{l}_k+{\sigma }_B^2)}}\exp (-\frac{{(w-x_k-{\mu }_{\lambda ,k}{l}_k)}^2}{2l_k({\sigma }_{\lambda ,k}^2{l}_k+{\sigma }_B^2)})$$

(30)

However, due to the influence of ME, $x_k$ is an implicit state and there is uncertainty in the estimated $x_k$ , which should be considered in RUL prediction. To solve this problem, Si et al. ⁴¹ deduced the PDF of RUL with considering the state estimation as follows.

$$\begin{array}{cc}\multicolumn{1}{c}{f_{L_k|{\mathit{Y}}_{1:k}}(l_k|{\mathit{Y}}_{1:k})}& =\frac{w_{1,k}(B_k^2{\kappa }_{\lambda ,k}^2+C_k)-A_k{B}_k{w}_{2,k}{\kappa }_{\lambda ,k}^2-A_k{C}_k{\hat{\lambda }}_{k|k}}{C_k\sqrt{2\pi {(B_k^2{\kappa }_{\lambda ,k}^2+C_k)}^3}}\hfill \\ \multicolumn{1}{c}{}& \exp (-\frac{{(w-{\hat{x}}_{k|k}-{\hat{\lambda }}_{k|k}{l}_k)}^2}{2(B_k^2{\kappa }_{\lambda ,k}^2+C_k)})\hfill \end{array}$$

(31)

where

$$\begin{array}{cc}\multicolumn{1}{c}{w_{1,k}}& =(w-{\hat{x}}_{k|k}+{\rho }_k\frac{{\kappa }_{x,k}}{{\kappa }_{\lambda ,k}}{\hat{\lambda }}_{k|k}){\sigma }_B^2,\hfill \\ \multicolumn{1}{c}{w_{2,k}}& =\frac{w_{1,k}}{{\sigma }_B^2},A_k={\rho }_k\frac{{\kappa }_{x,k}}{{\kappa }_{\lambda ,k}}{\sigma }_B^2-{\sigma }_{x_k|\lambda ,k}^2,\hfill \\ \multicolumn{1}{c}{B_k}& ={\rho }_k\frac{{\kappa }_{x,k}}{{\kappa }_{\lambda ,k}}+l_k,C_k={\sigma }_{x_k|\lambda ,k}^2+{\sigma }_B^2{l}_k\hfill \end{array}$$

Prior parameters estimation of the implicit linear Wiener degradation process based on multi-source information

Prior parameters estimation based on failure time data

In practical engineering applications, we often encounter imperfect prior information. The traditional MLE method may not obtain the correct prior parameters estimation results, or even estimate the prior parameters. In order to solve this problem, Tang et al. ³³ proposed a RUL prediction method that reasonably fuses failure time data and field degradation information. However, Tang et al. ³³ did not consider the influence of ME, which may increase the uncertainty of RUL prediction. Therefore, based on Tang et al., ³³ this subsection further proposes a method with fusing failure time data for the implicit linear Wiener degradation process. The flow chart is shown in Figure 1.

Figure 1.

The flow chart of RUL prediction based on failure time data.

The flow chart can be explained as follows. Firstly, the fixed parameters ${\sigma }_B^2$ and ${\sigma }_{\epsilon }^2$ are estimated based on field degradation data. Then, the prior distribution of ${\mu }_{\lambda }$ and ${\sigma }_{\lambda }^2$ are estimated by the EM algorithm based on failure time data. Next, the drift parameters ${\mu }_{\lambda }$ and ${\sigma }_{\lambda }^2$ are online updated by Kalman filter algorithm based on field degradation data, as shown in equation (28). Finally, the RUL is predicted, as shown in equation (31). According to the natures of parameters estimation obtained in Section “Natures of parameters estimation for the implicit linear Wiener degradation process,” the estimation accuracy of unknown parameters can be improved by increasing the number of field degradation data. During the implementation of this prediction method, once a new field degradation data is detected, it needs to carry out the estimation process of ${\sigma }_B^2$ and ${\sigma }_{\epsilon }^2$ in Figure 1 again.

Define ${\mathit{y}}_{0:m}=\{y_0,y_1,···,y_m\}$ as field degradation data, ${\mathit{T}}_{1:n^{\prime }}=\{T_1,T_2,···,T_{n^{\prime }}\}$ as the failure time data of $n^{\prime }$ units, that is, the time when the degradation value first hits the failure threshold $w$ , then, the prior parameters estimation for RUL prediction method with fusing failure time data mainly includes the following two steps.

• Step 1: Estimating the fixed parameters ${\sigma }_B^2$ and ${\sigma }_{\epsilon }^2$ by using the traditional MLE method based on field degradation data.

For the degradation process modeling of the single assessed unit, it is not necessary to consider the influence of random effects when establishing the degradation model of the implicit linear Wiener degradation process. Since there is no individual difference, $\lambda$ can be regarded as a fixed parameter. Based on equation (1), let $\mathrm{\Delta }{y}_j=y_j-y_{j-1}$ , $\mathrm{\Delta }{t}_j=t_j-t_{j-1}$ , $1\le j\le m$ , $\mathrm{\Delta }\mathit{y}=(\mathrm{\Delta }{y}_1,\mathrm{\Delta }{y}_2,···,\mathrm{\Delta }{y}_m{)}^{\prime }$ , and $\mathrm{\Delta }\mathit{t}=(\mathrm{\Delta }{t}_1,\mathrm{\Delta }{t}_2,···,\mathrm{\Delta }{t}_m{)}^{\prime }$ . Then $\mathrm{\Delta }\mathit{y}~\mathrm{MVN}(\lambda \mathrm{\Delta }\mathit{t},{\sigma }_B^2\mathbf{\Omega }\mathbf{+}{\mathbf{\sigma }}_{\mathbf{\epsilon }}^2\mathit{F})$ , where $\mathbf{\Omega }\mathbf{=}\mathrm{diag}(\mathbf{\Delta }{\mathbf{t}}_1,\mathbf{\Delta }{\mathbf{t}}_2,\mathbf{·}\mathbf{·}\mathbf{·},\mathbf{\Delta }{\mathbf{t}}_{\mathbf{m}})$ and the elements of $\mathit{F}$ are shown in equation (4). Let $D=\mathbf{\Omega }\mathbf{+}\mathbf{\varphi }\mathit{F}$ , where $\varphi ={\sigma }_{\epsilon }^2/{\sigma }_B^2$ . Then the log likelihood function of a single unit can be expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{}& \ln L(\lambda ,{\sigma }_B^2,\varphi |\mathit{Y})=-\frac{m}{2}\ln 2\pi -\frac{1}{2}m\ln ({\sigma }_B^2)\hfill \\ \multicolumn{1}{c}{}& -\frac{1}{2}\ln \left|\mathit{D}\right|-\frac{(\mathrm{\Delta }\mathit{y}-\lambda \mathrm{\Delta }\mathit{t}){\mathit{D}}^{-1}(\mathrm{\Delta }\mathit{y}-\lambda \mathrm{\Delta }\mathit{t}{)}^{\prime }}{2{\sigma }_B^2}\hfill \end{array}$$

(32)

Taking the partial derivatives of $\lambda$ and ${\sigma }_B^2$ for equation (32) gives

$$\frac{\partial \ln L(\lambda ,{\sigma }_B^2,\varphi |\mathbf{Y})}{\partial \lambda }=-\frac{-2\mathrm{\Delta }\mathit{y}{\mathit{D}}^{-1}\mathrm{\Delta }{\mathit{t}}^{\prime }+2\lambda \mathrm{\Delta }\mathit{t}{\mathit{D}}^{-1}\mathrm{\Delta }{\mathit{t}}^{\prime }}{2{\sigma }_B^2}$$

(33)

$$\frac{\partial \ln L(\lambda ,{\sigma }_B^2,\varphi |\mathbf{Y})}{\partial {\sigma }_B^2}=-\frac{1}{2{\sigma }_B^2}m+\frac{(\mathrm{\Delta }\mathit{y}-\lambda \mathrm{\Delta }\mathit{t}){\mathit{D}}^{-1}(\mathrm{\Delta }\mathit{y}-\lambda \mathrm{\Delta }\mathit{t}{)}^{\prime }}{2{({\sigma }_B^2)}^2}$$

(34)

By setting equations (33) and (34) to zero, the restricted estimation for $\lambda$ and ${\sigma }_B^2$ can be obtained as follows.

$$\hat{\lambda }(\hat{\varphi })=\frac{\mathrm{\Delta }\mathit{y}{\mathit{D}}^{-1}\mathrm{\Delta }{\mathit{t}}^{\prime }}{\mathrm{\Delta }\mathit{t}{\mathit{D}}^{-1}\mathrm{\Delta }{\mathit{t}}^{\prime }}$$

(35)

$${\hat{\sigma }}_B^2(\hat{\varphi })=\frac{(\mathrm{\Delta }\mathit{y}-\lambda \mathrm{\Delta }\mathit{t}){\mathit{D}}^{-1}(\mathrm{\Delta }\mathit{y}-\lambda \mathrm{\Delta }\mathit{t}{)}^{\prime }}{m}$$

(36)

By substituting equation (36) into equation (32), the expression of the profile log-likelihood function regarding $\varphi$ can be obtained as follows

$$\begin{array}{cc}\multicolumn{1}{c}{\ln L(\varphi |\mathbf{Y})}& =-\frac{m}{2}\ln 2\pi -\frac{1}{2}\ln \left|\mathit{D}\right|-\frac{1}{2}m\hfill \\ \multicolumn{1}{c}{}& \ln \left(\frac{(\mathrm{\Delta }\mathit{y}-\lambda \mathrm{\Delta }\mathit{t}){\mathit{D}}^{-1}(\mathrm{\Delta }\mathit{y}-\lambda \mathrm{\Delta }\mathit{t}{)}^{\prime }}{m}\right)-\frac{m}{2}\hfill \end{array}$$

(37)

The estimation of $\varphi$ can be obtained by maximizing equation (37) by the MATLAB function “FMINSEARCH.” Then, $\lambda$ and ${\sigma }_B^2$ can be obtained by substituting $\varphi$ into equations (35) and (36), and ${\sigma }_{\epsilon }^2$ can be further obtained.

• Step 2: Transforming the failure time data into prior distribution of ${\mu }_{\lambda }$ and ${\sigma }_{\lambda }^2$ by using EM algorithm.

According to conclusions in Chien-Yu and Sheng-Tsaing, ⁴⁴ the lifetime distribution function of the implicit linear Wiener degradation process is not affected by ME. Then, the prior distribution of drift parameter $\lambda$ under the estimated ${\sigma }_B^2$ and ${\mathit{T}}_{1:n^{\prime }}=\{T_1,T_2,···,T_{n^{\prime }}\}$ can be obtained, as shown in the Lemma 3.

• Lemma 3 ³³ : Define $T_v$ as the FHT of vth unit. According to Bayesian theory, given the prior distribution of the $\lambda$ , the posterior drift parameter $\lambda$ under $T_v$ follows normal distribution, that is, $\mathbf{\lambda }\mathbf{|}{\mathbf{T}}_{\mathbf{v}}\mathbf{~}\mathbf{N}\mathbf{(}{\mathbf{\mu }}_{\mathbf{\lambda }\mathbf{,}\mathbf{v}}\mathbf{,}{\mathbf{\sigma }}_{\mathbf{\lambda }\mathbf{,}\mathbf{v}}^2\mathbf{)}$ . Then, we have

$${\mu }_{\lambda ,v}=\frac{w{\sigma }_{\lambda }^2+{\mu }_{\lambda }{\sigma }_B^2}{T_v{\sigma }_{\lambda }^2+{\sigma }_B^2},{\sigma }_{\lambda ,v}^2=\frac{{\sigma }_{\lambda }^2{\sigma }_B^2}{T_v{\sigma }_{\lambda }^2+{\sigma }_B^2}$$

(38)

Based on Lemma 3, given ${\mathit{T}}_{1:n^{\prime }}=\{T_1,T_2,···,T_{n^{\prime }}\}$ and $w$ , where $1\le v\le n^{\prime }$ . Then, the profile log-likelihood function of failure time ${\mathit{T}}_{1:n^{\prime }}$ can be expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{\ln L({\mu }_{\lambda },{\sigma }_{\lambda }^2|{\mathit{T}}_{1:n^{\prime }},\mathbf{\lambda }\mathbf{)}}& \mathbf{=}\mathbf{-}\frac{{\mathbf{n}}^{\mathbf{’}}\mathbf{ln}2\mathbf{\pi }}{2}\mathbf{+}{\mathbf{n}}^{\mathbf{’}}\mathbf{ln}\mathbf{w}\mathbf{-}\frac{{\mathbf{n}}^{\mathbf{’}}}{2}\mathbf{ln}{\mathbf{\sigma }}_{\mathbf{B}}^2\hfill \\ \multicolumn{1}{c}{}& \mathbf{-}\frac{3}{2}\sum _{\mathbf{v}\mathbf{=}1}^{{\mathbf{n}}^{\mathbf{’}}}\mathbf{ln}{\mathbf{T}}_{\mathbf{v}}\mathbf{-}\sum _{\mathbf{v}\mathbf{=}1}^{{\mathbf{n}}^{\mathbf{’}}}\frac{{\mathbf{(}\mathbf{w}\mathbf{-}{\mathbf{\lambda }}_{\mathbf{v}}{\mathbf{T}}_{\mathbf{v}}\mathbf{)}}^2}{2{\mathbf{\sigma }}_{\mathbf{B}}^2{\mathbf{T}}_{\mathbf{v}}}\hfill \\ \multicolumn{1}{c}{}& \mathbf{-}\frac{{\mathbf{n}}^{\mathbf{’}}}{2}\mathbf{ln}{\mathbf{\sigma }}_{\mathbf{\lambda }}^2\mathbf{-}\frac{1}{2{\mathbf{\sigma }}_{\mathbf{\lambda }}^2}\sum _{\mathbf{v}\mathbf{=}1}^{{\mathbf{n}}^{\mathbf{’}}}{\mathbf{(}{\mathbf{\lambda }}_{\mathbf{v}}\mathbf{-}{\mathbf{\mu }}_{\mathbf{\lambda }}\mathbf{)}}^2\hfill \end{array}$$

(39)

Suppose that ${\hat{\Theta }}^{(i)}=\{{\hat{\mu }}_{\lambda }^{(i)},{\hat{\sigma }}_{\lambda }^{2(i)}\}$ are the parameter estimation results based on EM algorithm at the ith step, then the (i + 1)th step of EM algorithm can be divided into E-step and M-step as follows.

• E-step: Calculating the expectation of equation (39) as follows.

$$\begin{array}{cc}\multicolumn{1}{c}{\ln L(\mathbf{\Theta }\mathbf{|}{\mathit{T}}_{1\mathbf{:}{\mathbf{n}}^{\mathbf{’}}}\mathbf{,}{\hat{}\mathbf{\Theta }}^{\mathbf{(}\mathbf{i}\mathbf{)}}\mathbf{)}}& \mathbf{=}{\mathbf{E}}_{\mathbf{\lambda }\mathbf{|}{\mathbf{T}}_{1\mathbf{:}{\mathbf{n}}^{\mathbf{’}}}\mathbf{,}{\hat{}\mathbf{\Theta }}^{\mathbf{(}\mathbf{i}\mathbf{)}}}\mathbf{[}\mathbf{ln}\mathbf{L}\mathbf{(}\mathbf{\Theta }\mathbf{|}{\mathit{T}}_{1\mathbf{:}{\mathbf{n}}^{\mathbf{’}}}\mathbf{)}\mathbf{]}\hfill \\ \multicolumn{1}{c}{}& \mathbf{=}\mathbf{-}\frac{{\mathbf{n}}^{\mathbf{’}}\mathbf{ln}2\mathbf{\pi }}{2}\mathbf{+}{\mathbf{n}}^{\mathbf{’}}\mathbf{ln}\mathbf{w}\mathbf{-}\frac{{\mathbf{n}}^{\mathbf{’}}}{2}\mathbf{ln}{\mathbf{\sigma }}_{\mathbf{B}}^2\hfill \\ \multicolumn{1}{c}{}& \mathbf{-}\frac{3}{2}\sum _{\mathbf{v}\mathbf{=}1}^{{\mathbf{n}}^{\mathbf{’}}}\mathbf{ln}{\mathbf{T}}_{\mathbf{v}}\mathbf{-}\sum _{\mathbf{v}\mathbf{=}1}^{{\mathbf{n}}^{\mathbf{’}}}\frac{{\mathbf{(}\mathbf{w}\mathbf{-}{\mathbf{\mu }}_{\mathbf{\lambda }\mathbf{,}\mathbf{v}}{\mathbf{T}}_{\mathbf{v}}\mathbf{)}}^2\mathbf{+}{\mathbf{\sigma }}_{\mathbf{\lambda }\mathbf{,}\mathbf{v}}^2{\mathbf{T}}_{\mathbf{v}}^2}{2{\mathbf{\sigma }}_{\mathbf{B}}^2{\mathbf{T}}_{\mathbf{v}}}\hfill \\ \multicolumn{1}{c}{}& \mathbf{-}\frac{{\mathbf{n}}^{\mathbf{’}}}{2}\mathbf{ln}{\mathbf{\sigma }}_{\mathbf{\lambda }}^2\mathbf{-}\frac{1}{2{\mathbf{\sigma }}_{\mathbf{\lambda }}^2}\sum _{\mathbf{v}\mathbf{=}1}^{{\mathbf{n}}^{\mathbf{’}}}\mathbf{(}{\mathbf{(}{\mathbf{\mu }}_{\mathbf{\lambda }\mathbf{,}\mathbf{v}}\mathbf{-}{\mathbf{\mu }}_{\mathbf{\lambda }}\mathbf{)}}^2\mathbf{+}{\mathbf{\sigma }}_{\mathbf{\lambda }\mathbf{,}\mathbf{v}}^2\mathbf{)}\hfill \end{array}$$

(40)

where ${\mu }_{\lambda ,v}$ and ${\sigma }_{\lambda ,v}^2$ are the posterior distribution of $\lambda$ under ${\mathit{T}}_{1:n^{\prime }}$ and $\hat{}{\Theta }^{(i)}$ .

• M-step: Maximizing equation (40).

$${\hat{\Theta }}^{(i+1)}=\arg \underset{\Theta }{max}L(\Theta |{\hat{\Theta }}^{(i)})$$

Taking the partial derivatives of ${\mu }_{\lambda }$ and ${\sigma }_{\lambda }^2$ for equation (40) and setting them to zero, we can obtain the expression of (i + 1)th iteration as follows.

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\hat{\mu }}_{\lambda }=\frac{1}{n^{\prime }}\sum _{v=1}^{n^{\prime }}{\mu }_{\lambda ,v},\hfill \\ \multicolumn{1}{c}{}& {\hat{\sigma }}_{\lambda }^2=\frac{1}{n^{\prime }-1}\sum _{v=1}^{n^{\prime }}[{({\mu }_{\lambda ,v}-{\hat{\mu }}_{\lambda })}^2+{\sigma }_{\lambda ,v}^2].\hfill \end{array}$$

(41)

Then, the above E-step and M-step are iterated until $||{\mathbf{\Theta }}^{(i+1)}-{\mathbf{\Theta }}^{(i)}||$ is sufficiently small.

Prior parameters estimation based on multi source information

When prior degradation information is scarce and there are both prior historical degradation data and failure time data, according to the natures of parameters estimation obtained in Section “Natures of parameters estimation for the implicit linear Wiener degradation process,” we can fuse as much degradation data as possible to further improve the accuracy of unknown parameters estimation. Recently, to solve this problem, Wang et al. ³⁹ proposed a RUL prediction method with fusing multi-source information and considering random effects. However, Wang et al. ³⁹ did not consider the influence of ME. In this subsection, we further extend this method to the implicit linear Wiener degradation process, and the specific flow chart is shown in Figure 2.

Figure 2.

The flow chart of RUL prediction based on multi source information.

The flow chart can be explained as follows. Firstly, the fixed parameters ${\sigma }_B^2$ and ${\sigma }_{\epsilon }^2$ are estimated based on historical degradation data. Secondly, the prior distribution of ${\mu }_{\lambda }$ and ${\sigma }_{\lambda }^2$ are estimated by EM algorithm based on historical degradation data and failure time data. Then, the drift parameter ${\mu }_{\lambda }$ and ${\sigma }_{\lambda }^2$ are online updated by Kalman filter algorithm based on field degradation data, as shown in equation (28). Finally, the RUL is predicted, according to equation (31).

Suppose that there are $n$ units for testing and each unit are measured at the same times $t_1,t_2,···,t_m$ . Let $y_{i,j}$ denote the degradation state measured by the ith unit at time $t_j$ , $\mathrm{\Delta }{y}_{i,j}=y_i(t_j)-y_i(t_{j-1})$ , $\mathrm{\Delta }{t}_j=t_j-t_{j-1}$ , $\mathrm{\Delta }{\mathit{y}}_i=(\mathrm{\Delta }{y}_{i,1},\mathrm{\Delta }{y}_{i,2},···,\mathrm{\Delta }{y}_{i,m}{)}^{\prime }$ , $\mathrm{\Delta }\mathit{t}=(\mathrm{\Delta }{t}_1,\mathrm{\Delta }{t}_2,···,\mathrm{\Delta }{t}_m{)}^{\prime }$ , and $\mathit{Y}=(\mathrm{\Delta }{\mathit{y}}_1^{\prime },\mathrm{\Delta }{\mathit{y}}_2^{\prime },···,\mathrm{\Delta }{\mathit{y}}_n^{\prime }{)}^{\prime }$ , where $1\le i\le n$ , $1\le j\le m$ , and ${\lambda }_i$ denotes the degradation rate of the ith unit. Then, it can be obtained that $\mathrm{\Delta }{\mathit{y}}_i$ follows multivariate normal with mean ${\mu }_{\lambda }\mathrm{\Delta }\mathit{t}$ and variance ${\sigma }_{\lambda }^2\mathrm{\Delta }\mathit{t}\mathrm{\Delta }{\mathit{t}}^{\prime }+{\sigma }_B^2\mathbf{\Omega }\mathbf{+}{\mathbf{\sigma }}_{\mathbf{\epsilon }}^2\mathit{F}$ . In addition, it is assumed that there are $n^{\prime }$ groups of failure time data ${\mathit{T}}_{1:n^{\prime }}=\{T_1,T_2,···,T_{n^{\prime }}\}$ and the failure time data of the vth unit is $T_v$ , where, $1\le v\le n^{\prime }$ . Let $\mathit{D}=\mathbf{\Omega }\mathbf{+}\mathbf{\varphi }\mathit{F}$ , where $\varphi ={\sigma }_{\epsilon }^2/{\sigma }_B^2$ . Inspired by Wang et al., ³⁹ the log likelihood function of $\mathbf{\Theta }\mathbf{=}\mathbf{(}{\mathbf{\mu }}_{\mathbf{\lambda }}\mathbf{,}{\mathbf{\sigma }}_{\mathbf{\lambda }}^2\mathbf{)}$ can be obtained as follows.

$$\begin{array}{cc}\multicolumn{1}{c}{\ln L(\mathbf{\Theta }\mathbf{|}\mathbf{Y}\mathbf{,}{\mathit{T}}_{1\mathbf{:}{\mathbf{n}}^{\mathbf{’}}}\mathbf{,}\mathbf{\lambda }\mathbf{)}}& \mathbf{=}\mathbf{-}\frac{{\mathbf{n}}^{\mathbf{’}}\mathbf{ln}2\mathbf{\pi }}{2}\mathbf{+}{\mathbf{n}}^{\mathbf{’}}\mathbf{ln}\mathbf{w}\mathbf{-}\frac{{\mathbf{n}}^{\mathbf{’}}}{2}\mathbf{ln}{\mathbf{\sigma }}_{\mathbf{B}}^2\hfill \\ \multicolumn{1}{c}{}& \mathbf{-}\frac{3}{2}\sum _{\mathbf{v}\mathbf{=}1}^{{\mathbf{n}}^{\mathbf{’}}}\mathbf{ln}{\mathbf{T}}_{\mathbf{v}}\mathbf{-}\sum _{\mathbf{v}\mathbf{=}1}^{{\mathbf{n}}^{\mathbf{’}}}\frac{{\mathbf{(}\mathbf{w}\mathbf{-}{\mathbf{\lambda }}_{\mathbf{v}}{\mathbf{T}}_{\mathbf{v}}\mathbf{)}}^2}{2{\mathbf{\sigma }}_{\mathbf{B}}^2{\mathbf{T}}_{\mathbf{v}}}\hfill \\ \multicolumn{1}{c}{}& \mathbf{-}\frac{{\mathbf{n}}^{\mathbf{’}}\mathbf{ln}2\mathbf{\pi }}{2}\mathbf{-}\frac{{\mathbf{n}}^{\mathbf{’}}}{2}\mathbf{ln}{\mathbf{\sigma }}_{\mathbf{\lambda }}^2\mathbf{-}\frac{1}{2{\mathbf{\sigma }}_{\mathbf{\lambda }}^2}\sum _{\mathbf{v}\mathbf{=}1}^{{\mathbf{n}}^{\mathbf{’}}}{\mathbf{(}{\mathbf{\lambda }}_{\mathbf{v}}\mathbf{-}{\mathbf{\mu }}_{\mathbf{\lambda }}\mathbf{)}}^2\hfill \\ \multicolumn{1}{c}{}& \mathbf{-}\frac{\mathbf{mn}\mathbf{ln}2\mathbf{\pi }}{2}\mathbf{-}\frac{\mathbf{mn}\mathbf{ln}{\mathbf{\sigma }}_{\mathbf{B}}^2}{2}\mathbf{-}\frac{\mathbf{n}}{2}\mathbf{ln}\mathbf{\left|}\mathit{D}\mathbf{\right|}\hfill \\ \multicolumn{1}{c}{}& \mathbf{-}\frac{1}{2{\mathbf{\sigma }}_{\mathbf{B}}^2}\sum _{\mathbf{i}\mathbf{=}1}^{\mathbf{n}}{\mathbf{(}\mathbf{\Delta }{\mathit{y}}_{\mathbf{i}}\mathbf{-}{\mathbf{\lambda }}_{\mathbf{i}}\mathbf{\Delta }\mathit{t}\mathbf{)}}^{\mathbf{’}}{\mathit{D}}^{\mathbf{-}1}\mathbf{(}\mathbf{\Delta }{\mathit{y}}_{\mathbf{i}}\mathbf{-}{\mathbf{\lambda }}_{\mathbf{i}}\mathbf{\Delta }\mathit{t}\mathbf{)}\hfill \\ \multicolumn{1}{c}{}& \mathbf{-}\frac{\mathbf{n}}{2}\mathbf{ln}2\mathbf{\pi }\mathbf{-}\frac{\mathbf{n}}{2}\mathbf{ln}{\mathbf{\sigma }}_{\mathbf{\lambda }}^2\mathbf{-}\frac{1}{2{\mathbf{\sigma }}_{\mathbf{\lambda }}^2}\sum _{\mathbf{i}\mathbf{=}1}^{\mathbf{n}}{\mathbf{(}{\mathbf{\lambda }}_{\mathbf{i}}\mathbf{-}{\mathbf{\mu }}_{\mathbf{\lambda }}\mathbf{)}}^2\hfill \end{array}$$

(42)

Suppose that ${\hat{\Theta }}^{(k)}=\{{\hat{\mu }}_{\lambda ,i}^{(k)},{\hat{\sigma }}_{\lambda ,i}^{2(k)},{\hat{\mu }}_{\lambda ,v}^{(k)},{\hat{\sigma }}_{\lambda ,v}^{2(k)}\}$ are the parameters estimation results based on EM algorithm at the kth step, then the (k + 1)th step of EM algorithm can be divided into E-step and M-step as follows.

• E-step: Calculating the expectation of equation (42)

$$\begin{array}{cc}\multicolumn{1}{c}{E(\ln L(\mathbf{\Theta }\mathbf{|}\mathbf{Y}\mathbf{,}{\mathit{T}}_{1\mathbf{:}{\mathbf{n}}^{\mathbf{’}}}\mathbf{,}\mathbf{\lambda }\mathbf{)})}& =-\frac{n^{\prime }\ln 2\pi }{2}+n^{\prime }\ln w-\frac{n^{\prime }}{2}\ln {\sigma }_B^2-\frac{3}{2}\sum _{v=1}^{n^{\prime }}\ln T_v-\sum _{v=1}^{n^{\prime }}\frac{{(w-{\mu }_{\lambda ,v}{T}_v)}^2+{\sigma }_{\lambda ,v}^2{T}_v^2}{2{\sigma }_B^2{T}_v}\hfill \\ \multicolumn{1}{c}{}& -\frac{n^{\prime }\ln 2\pi }{2}-\frac{n^{\prime }}{2}\ln {\sigma }_{\lambda }^2-\frac{1}{2{\sigma }_{\lambda }^2}\sum _{v=1}^{n^{\prime }}({({\mu }_{\lambda ,v}-{\mu }_{\lambda })}^2+{\sigma }_{\lambda ,v}^2)-\frac{\mathrm{mn}\ln 2\pi }{2}-\frac{\mathrm{mn}}{2}\ln {\sigma }_B^2\hfill \\ \multicolumn{1}{c}{}& -\frac{n}{2}\ln \left|\mathit{D}\right|-\frac{1}{2{\sigma }_B^2}\sum _{i=1}^n[{(\mathrm{\Delta }{\mathit{y}}_i-{\mu }_{\lambda ,i}\mathrm{\Delta }\mathit{t})}^{\prime }{\mathit{D}}^{-1}(\mathrm{\Delta }{\mathit{y}}_i-{\mu }_{\lambda ,i}\mathrm{\Delta }\mathit{t})+{\sigma }_{\lambda ,i}^2\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}]\hfill \\ \multicolumn{1}{c}{}& -\frac{1}{2{\sigma }_{\lambda }^2}\sum _{i=1}^n[{({\mu }_{\lambda ,i}-{\mu }_{\lambda })}^2+{\sigma }_{\lambda ,i}^2]-\frac{n}{2}\ln 2\pi -\frac{n}{2}\ln {\sigma }_{\lambda }^2\hfill \end{array}$$

(43)

where

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\mu }_{\lambda ,i}=E({\lambda }_i|\mathrm{\Delta }{\mathit{y}}_i,{\hat{}\Theta }^{(k)})=\frac{\mathrm{\Delta }{{\mathit{y}}_i}^{\prime }{\mathit{D}}^{(k)-1}\mathrm{\Delta }\mathit{t}·\sigma {{}_{\lambda }^2}^{(k)}+{{\mu }_{\lambda }}^{(k)}{\sigma }_B^{2(k)}}{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{(k)-1}\mathrm{\Delta }\mathit{t}·\sigma {{}_{\lambda }^2}^{(k)}+{\sigma }_B^{2(k)}},\hfill \\ \multicolumn{1}{c}{}& {\sigma }_{\lambda ,i}^2=\mathrm{var}({\lambda }_i|\mathrm{\Delta }{\mathit{y}}_i,{\hat{}\Theta }^{(k)})=\frac{\sigma {{}_{\lambda }^2}^{(k)}{\sigma }_B^{2(k)}}{\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{(k)-1}\mathrm{\Delta }\mathit{t}·\sigma {{}_{\lambda }^2}^{(k)}+{\sigma }_B^{2(k)}}.\hfill \end{array}$$

(44)

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\mu }_{\lambda ,v}=E({\lambda }_v|T_v,{\hat{}\Theta }^{(k)})=\frac{w\sigma {{}_{\lambda }^2}^{(k)}+{{\mu }_{\lambda }}^{(k)}\sigma {{}_B^2}^{(k)}}{T_v\sigma {{}_{\lambda }^2}^{(k)}+\sigma {{}_B^2}^{(k)}},\hfill \\ \multicolumn{1}{c}{}& {\sigma }_{\lambda ,v}^2=\mathrm{var}({\lambda }_v|T_v,{\hat{}\Theta }^{(k)})=\frac{\sigma {{}_B^2}^{(k)}\sigma {{}_{\lambda }^2}^{(k)}}{T_v\sigma {{}_{\lambda }^2}^{(k)}+\sigma {{}_B^2}^{(k)}}.\hfill \end{array}$$

(45)

• M-step: Maximizing equation (43)

Taking the partial derivatives of ${\mu }_{\lambda }$ and ${\sigma }_{\lambda }^2$ in equation (43) gives

$$\begin{array}{cc}\multicolumn{1}{c}{\frac{\partial E(\ln L(\mathbf{\Theta }\mathbf{|}\mathbf{Y}\mathbf{,}{\mathit{T}}_{1\mathbf{:}{\mathbf{n}}^{\mathbf{’}}}\mathbf{,}\mathbf{\lambda }\mathbf{)})}{\partial {\mu }_{\lambda }}}& =-\sum _{v=1}^{n^{\prime }}{\mu }_{\lambda ,v}+n^{\prime }\hfill \\ \multicolumn{1}{c}{}& -\sum _{i=1}^n{\mu }_{\lambda ,i}+n\hfill \end{array}$$

(46)

$$\begin{array}{cc}\multicolumn{1}{c}{}& \frac{\partial E(\ln L(\mathbf{\Theta }\mathbf{|}\mathbf{Y}\mathbf{,}{\mathit{T}}_{1\mathbf{:}{\mathbf{n}}^{\mathbf{’}}}\mathbf{,}\mathbf{\lambda }\mathbf{)})}{\partial {\sigma }_{\lambda }^2}=-\frac{n^{\prime }}{2{\sigma }_{\lambda }^2}-\frac{n}{2{\sigma }_{\lambda }^2}\hfill \\ \multicolumn{1}{c}{}& +\frac{\sum _{v=1}^{n^{\prime }}({({\mu }_{\lambda ,v}-{\mu }_{\lambda })}^2+{\sigma }_{\lambda ,v}^2)+\sum _{i=1}^n[{({\mu }_{\lambda ,i}-{\mu }_{\lambda })}^2+{\sigma }_{\lambda ,i}^2]}{2{\left({\sigma }_{\lambda }^2\right)}^2}\hfill \end{array}$$

(47)

Then, by setting equations (46) and (47) to zero, we can obtain the expression of the (k + 1)th iteration as follows.

$$\begin{array}{cc}\multicolumn{1}{c}{}& {\hat{\mu }}_{\lambda }=\frac{\sum _{v=1}^{n^{\prime }}{\mu }_{\lambda ,v}+\sum _{i=1}^n{\mu }_{\lambda ,i}}{n^{\prime }+n},\hfill \\ \multicolumn{1}{c}{}& {\hat{\sigma }}_{\lambda }^2=\frac{\sum _{v=1}^{n^{\prime }}[{({\mu }_{\lambda ,v}-{\mu }_{\lambda })}^2+{\sigma }_{\lambda ,v}^2]+\sum _{i=1}^n[{({\mu }_{\lambda ,i}-{\mu }_{\lambda })}^2+{\sigma }_{\lambda ,i}^2]}{n^{\prime }+n}.\hfill \end{array}$$

(48)

Similarly, the above E-step and M-step are iterated until $||{\mathbf{\Theta }}^{(k+1)}-{\mathbf{\Theta }}^{(k)}||$ is sufficiently small.

• Remark 2: The fusion of multi-source information refers to the fusion of historical degradation data information, failure time data information, and field degradation data information of the assessed equipment. However, when the prior degradation information is imperfect, the estimation results of the parameters may be far from the actual value, thus lead to the low accuracy of RUL prediction, if the fixed parameters is estimated by the traditional MLE method. For this special case, when using the method with fusing multi-source information to predict RUL, the historical degradation data can be converted into failure time data, and then the RUL can be predicted by fusing the failure time data information and the field degradation data information of the assessed unit.

Experimental studies

In this section, the effectiveness and superiority of the proposed method are verified by some simulation experiments and two practical case studies. Firstly, the Monte Carlo algorithm is used to generate the simulation data of the implicit linear Wiener degradation process to verify the natures of parameters estimation obtained in Section “Natures of parameters estimation for the implicit linear Wiener degradation process.” Secondly, the lithium battery degradation data is used to verify the superiority of the method with fusing failure time data, and the laser degradation data is used to verify the superiority of the method with fusing multi-source information.

Simulation experiments

Tang et al. ³³ analyzed the natures of parameters estimation based on the basic linear Wiener degradation process model, and concluded that the accuracy of parameters estimation is proportional to the number of units and the detection times of degradation data. According to the conclusions in Section “Natures of parameters estimation for the implicit linear Wiener degradation process,” when the degradation data have ME, the accuracy of parameters estimation of the implicit Wiener degradation process model also has similar natures. In this subsection, we further use simulation experiments to verify these natures.

The parameters are assumed as ${\mu }_{\lambda }=1$ , ${\sigma }_{\lambda }^2=0.04$ , ${\sigma }_B^2=0.25$ , and ${\sigma }_{\epsilon }^2=0.49$ . Then, let the number of units $n=80$ , detection times $m=40$ , detection interval $\mathrm{\Delta }t=1$ , where each unit is detected at the same time. The simulated degradation paths of some units are shown in Figure 3. In order to verify the influence of the number of units on the accuracy of parameters estimation, we use 5, 20, and 80 groups of the degradation data to estimate parameters. In addition, in order to verify the influence of the length of detection time and detection times on the estimation accuracy, once a new degradation data is generated, the parameters are estimated again. Significantly, since the simulation data is generated randomly, if the parameters are estimated only once with the generated data, the results may not fully reflect the relationship between the degradation data and the accuracy of parameters estimation. Therefore, we use the same method to generate degradation data by 20 times, and use the same parameters estimation method every time. The changing curves of the estimated parameters are shown in Figures 4 to 7.

Figure 3.

Some degradation paths.

Figure 4.

The estimation of ${\mu }_{\lambda }$ with the change of the number $n$ of units: (a) n = 5, (b) n = 20, and (c) n = 80.

Figure 5.

The estimation of ${\sigma }_{\lambda }^2$ with the change of the number $n$ of units: (a) n = 5, (b) n = 20, and (c) n = 80.

Figure 6.

The estimation of ${\sigma }_B^2$ with the change of the number $n$ of units: (a) n = 5, (b) n = 20, and (c) n = 80.

Figure 7.

The estimation of ${\sigma }_{\epsilon }^2$ with the change of the number $n$ of units: (a) n = 5, (b) n = 20, and (c) n = 80.

It can be seen from Figures 4 and 5 that when the number $n$ of units is the same, the longer the length of detection time $t_m$ of degradation data, the closer the parameters estimation results are to the actual value. This is because when there is no ME, $\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}$ in equations (10), (12) and (13) is $t_m$ , and when there is no Brownian motion, $\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}$ in equations (10), (12) and (13) is $(t_1:t_m)(t_1:t_m{)}^{\prime }$ , that is, $\sum _1^m{t}_i^2$ . For the implicit linear Wiener degradation process, the value of $\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}$ is between them. Therefore, the longer the length of detection time, the bigger the value of $\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t}$ , and then the estimation results of ${\mu }_{\lambda }$ and ${\sigma }_{\lambda }^2$ can be more accurate. In addition, it can also be seen that when the detection time exceeds 10, the parameters estimation results do not change significantly with the increase of the length of time. Significantly, the accuracy of parameters estimation can be further improved by increasing the number of units $n$ . That is because ${\sigma }_B^2/(\mathrm{\Delta }{\mathit{t}}^{\prime }{\mathit{D}}^{-1}\mathrm{\Delta }\mathit{t})$ in equations (10), (12) and (13) are very small at this time, and the estimation accuracy of ${\mu }_{\lambda }$ and ${\sigma }_{\lambda }^2$ is mainly effected by the number of units $n$ . Therefore, when the prior degradation information is imperfect, using failure time data to estimate the prior drift parameter can improve its estimation accuracy. It is because the failure time data reflects the number of units $n$ . When there are both failure time data and historical degradation data, the estimation accuracy of prior drift parameter can further improve by reasonably fusing these two kinds degradation data.

Additionally, from Figures 6 and 7, it can be observed that when the number of units $n$ is the same, the changing curves of the estimated parameters gradually converge to the actual value with the increase of detection times. And when the detection times is the same, the more units are used, the more accurate the results of parameters estimation. The reason is that the values of equations (11) and (24) are inversely proportional to $n(m-1)$ . In other words, the estimation accuracy of ${\sigma }_B^2$ and ${\sigma }_{\epsilon }^2$ is proportional to the times of detection $m$ and the number of units $n$ . Therefore, when predicting the RUL in the case of imperfect prior degradation information, if the number of detection times of the field degradation data is large, such as the lithium battery degradation data described in Subsection “Predicting the RUL based on failure time data,” the fixed parameters ${\sigma }_B^2$ and ${\sigma }_{\epsilon }^2$ can be estimated without using the historical degradation data.

To sum up, simulation experiments verify the natures of parameters estimation obtained in Section “Natures of parameters estimation for the implicit linear Wiener degradation process,” and provide a theoretical basis for RUL prediction based on the method with fusing multi-source information and considering ME.

Predicting the RUL based on failure time data

In this subsection, we apply the degradation data published by NASA to verify the effectiveness of the RUL prediction method proposed in this paper when the prior degradation information is imperfect. The degradation data includes the capacity degradation data of four lithium batteries of the same type, which are repeatedly cycled under three operating modes of charging, discharging, and impedance measurement at room temperature. The capacity of lithium battery may regenerate after a long rest, which is also called relaxation effect. ⁴⁵ In practical operation, since the rest time of lithium battery is uncertain, RUL may not be able to be predicted accurately. Therefore, it is necessary to eliminate the relaxation effect in the original data. The degradation path after eliminating the relaxation effect is shown in Figure 8. For more details about how to eliminate the relaxation effect, see references Tang et al. ³⁰ and Jin et al. ⁴⁵ and references therein. Generally, when the battery capacity degrades to $70\mathrm{\%}~80\mathrm{\%}$ of the rated capacity, the battery is considered as failure. ³⁰ In this experiment, we assume that the failure threshold of lithium battery is 70% (1.4 Ahr) of the rated capacity, and select B0018 battery as the assessed unit, and other batteries are used for estimating the prior parameters.

Figure 8.

The degradation path of lithium battery capacity.

From Figure 8, we can observe that the capacity degradation curves of B0005 and B0007 have nonlinear characteristics, while the capacity degradation curves of B0006 and B0018 have linear characteristics, if the parameters are estimated directly using historical degradation data, the estimation results will be far from the actual value, thus the accuracy of RUL prediction will be reduced, ³³ which can be regarded as the case that the prior degradation information is imperfect. In fact, these four groups of degradation data not only contain degradation data information, but also contain failure time information, and their failure time information can be used to predict the RUL. ³³ For simplicity, the proposed RUL prediction method with fusing failure time data is referred to as ${\mathrm{M}}_0$ , the RUL prediction method based on the basic linear Wiener degradation process with fusing failure time data proposed in Tang et al. ³³ is referred to as ${\mathrm{M}}_1$ , the common RUL prediction method that only uses the historical degradation data to estimate the prior parameters of implicit linear Wiener degradation process model ⁴¹ is referred to as ${\mathrm{M}}_2$ . Based on the data with eliminating the relaxation effect, we can calculate that the failure times of B0005, B0006, and B00007 are 85.55, 69.5, and 110.3 cycles respectively, which are used as failure time data to calculate the prior drift parameter.

Figure 9 shows the RUL distributions of the assessed battery B00018 at some times, and Figure 10 shows the RUL distribution in the 44th cycle. It can be observed that the RUL distributions of ${\mathrm{M}}_0$ , ${\mathrm{M}}_1$ , and ${\mathrm{M}}_2$ can all cover the actual RUL, and the distributions of ${\mathrm{M}}_0$ is more concentrated on the actual RUL than ${\mathrm{M}}_1$ and ${\mathrm{M}}_2$ , which indicates that ${\mathrm{M}}_0$ has higher prediction accuracy. And it indicates that considering ME can improve the accuracy of RUL prediction by comparing the distributions of ${\mathrm{M}}_0$ and ${\mathrm{M}}_1$ . In addition, we can observe that although ${\mathrm{M}}_2$ consider ME, its RUL prediction accuracy is still poor by comparing the distributions of ${\mathrm{M}}_1$ and ${\mathrm{M}}_2$ , it indicates that the imperfect prior degradation information may lead to greater error in RUL prediction results.

Figure 9.

The estimated RULs at some cycle times.

Figure 10.

The estimated RULs at 44 cycle.

Then, we analyze the RUL prediction accuracy of the three prediction methods by comparing the results of parameters estimation. The parameters estimated in real time are shown in Figure 11. From Figure 11, it can be observed that ${\sigma }_B^2$ and ${\sigma }_{\epsilon }^2$ estimated by ${\mathrm{M}}_2$ are always bigger than ${\mathrm{M}}_0$ . The reason is that the capacity degradation curves of B0005 and B0007 have nonlinear degradation characteristics for a period of time, and the linear degradation model cannot effectively track the nonlinear degradation part, thus the nonlinear part is submerged into the parameters ${\sigma }_B^2$ and ${\sigma }_{\epsilon }^2$ estimated by ${\mathrm{M}}_2$ . In addition, it can be observed that ${\sigma }_B^2$ estimated by ${\mathrm{M}}_1$ are always bigger than ${\mathrm{M}}_0$ . The reason is that ${\mathrm{M}}_1$ does not consider the influence of ME, the degradation and measurement uncertainty of B0018 is all reflected by ${\sigma }_B^2$ . On the contrary, the uncertainty of B0018 is shared by ${\sigma }_B^2$ and ${\sigma }_{\epsilon }^2$ for ${\mathrm{M}}_0$ . Generally, the smaller the ${\sigma }_B^2$ , the sharper the RUL distribution. And then the uncertainty of RUL prediction results is smaller, which further indicates the superiority of the ${\mathrm{M}}_0$ method.

Figure 11.

The posterior parameters of ${\mathrm{M}}_0$ , ${\mathrm{M}}_1$ , and ${\mathrm{M}}_2$ at some cycle times: (a) ${\mu }_{\lambda }$ , (b) ${\sigma }_{\lambda }^2$ , (c) ${\sigma }_B^2$ , and (d) ${\sigma }_{\epsilon }^2$ .

In order to further distinguish the prediction accuracy of three prediction methods, the mean square errors (MSEs) of each model in each cycle are calculated. The calculation expression of MSEs at the kth time cycle can be written as follows.

$$\mathrm{MS}{E}_k={\int }_0^{\infty }{(l_k+t_k-T)}^2{f}_{L_k|{\mathbf{Y}}_{0:k}}(l_k|{\mathbf{Y}}_{1:k})d{l}_k$$

(49)

Figure 12 shows the MSEs of RULs of three prediction methods. It can be observed that the MSEs of ${\mathrm{M}}_0$ is always smaller than ${\mathrm{M}}_1$ and ${\mathrm{M}}_2$ , which indicates that the prediction accuracy of ${\mathrm{M}}_0$ is higher than ${\mathrm{M}}_1$ and ${\mathrm{M}}_2$ . Significantly, when the time of cycles of the assessed unit is less than 45, the MSEs of ${\mathrm{M}}_1$ are always smaller than ${\mathrm{M}}_2$ , while when the time of cycles of the assessed unit is more than 45, the MSEs of ${\mathrm{M}}_1$ are bigger than ${\mathrm{M}}_2$ , it indicates that it is necessary to consider the influence of ME. In general, when the prior degradation information is imperfect, the method with fusing failure time data can improve the accuracy of parameters estimation and RUL prediction.

Figure 12.

The MSEs of RULs of ${\mathrm{M}}_0$ , ${\mathrm{M}}_1$ , and ${\mathrm{M}}_2$ .

Predicting the RUL based on multi source information

In this subsection, we apply the laser degradation data provided by Meeker et al. ⁴⁶ to verify the effectiveness of the RUL prediction method with fusing multi-source information and considering the influence of ME. The degradation data includes 15 units, and the degradation state of each unit is detected at the same time. In addition, the detection interval is 250 h and each unit includes 16 measurement points as shown in Figure 13. In order to verify the superiority of the method with fusing multi-source information, it is assumed that among the 15 groups of degradation data, 1 group of degradation data is the field degradation data of the assessed unit, 11 groups of degradation data are historical degradation data, and 3 groups of degradation data can only obtain its failure time. ³⁹ For simplicity, the RUL prediction method with fusing multi-source information and considering the influence of ME proposed in this paper is referred to as ${\mathrm{M}}_3$ .

Figure 13.

Laser data degradation path.

From Figure 13, it can be observed that the degradation data has good linear degradation characteristics. In this subsection, we assume that the red solid line in Figure 13 is the assessed unit, the blue dotted line is the historical degradation data, which is used to estimate the unknown parameters of ${\mathrm{M}}_2$ , and the failure threshold is set as $w=10$ . ⁴⁴ Then the failure time of the three groups of black dotted lines in Figure 13 are 3611.7, 3707.4, and 3273.8 h respectively, which are used as failure time data for ${\mathrm{M}}_0$ and ${\mathrm{M}}_3$ .

Additionally, Figure 14 shows the RUL distributions of the assessed unit at some points, and Figure 15 shows the RUL distributions of the assessed unit at 3750 h. It can be seen that although the RUL distributions of ${\mathrm{M}}_0$ is higher than ${\mathrm{M}}_2$ and ${\mathrm{M}}_3$ , the RUL distributions of ${\mathrm{M}}_0$ is far from the actual RULs, and the RULs is even failed to be predicted at some points. However, the RUL distributions of ${\mathrm{M}}_2$ and ${\mathrm{M}}_3$ can all cover the actual RULs, and the RUL distributions of ${\mathrm{M}}_3$ is higher than ${\mathrm{M}}_2$ . It indicates that the prediction accuracy of ${\mathrm{M}}_3$ is better than ${\mathrm{M}}_2$ .

Figure 14.

The estimated RULs at some points.

Figure 15.

The estimated RULs at 3750 h.

Figure 16 shows the posterior parameters of three prediction methods. It can be observed that ${\sigma }_B^2$ estimated by ${\mathrm{M}}_0$ is close to zero before 3000 h and far from the estimated values of ${\mathrm{M}}_2$ and ${\mathrm{M}}_3$ . This is because the uncertainty of the field degradation data is small and the detection times $m$ in equation (11) is small, which may cause ${\sigma }_B^2$ to be incorrectly estimated, thus affect the online updating of the drift parameter of ${\mathrm{M}}_0$ under the Bayesian framework, and further leading to failure prediction of RUL. It indicates that for the units with less field degradation data, the method with only fusing failure time data may not be able to accurately predict the RUL. Figure 17 displays the MSEs of RULs of three prediction methods at different points. It can be observed that the MSEs calculated by ${\mathrm{M}}_3$ is all less than ${\mathrm{M}}_2$ . Although the MSEs calculated by ${\mathrm{M}}_0$ is less than ${\mathrm{M}}_2$ and ${\mathrm{M}}_3$ , it fails to predict RUL at many points. Therefore, the RUL prediction results of ${\mathrm{M}}_0$ is unstable and the reliability is poor. Overall ${\mathrm{M}}_3$ has the best prediction effect. In summary, the accuracy of parameters estimation can be improved by reasonably fusing as much prior degradation information as possible, and then the accuracy of RUL prediction can be improved.

Figure 16.

The posterior parameters of ${\mu }_{\lambda }$ , ${\sigma }_{\lambda }^2$ , and ${\sigma }_B^2$ at some points: (a) ${\mu }_{\lambda }$ , (b) ${\sigma }_{\lambda }^2$ , and (c) ${\sigma }_B^2$ .

Figure 17.

The MSEs of RULs at some points.

In order to further verify the RUL prediction method proposed in this paper, which fuses multi-source information and considering ME, can improve the accuracy of RUL prediction for degradation systems with ME, the RUL prediction method with fusing multi-source information based on the basic linear Wiener degradation process proposed in Wang et al. ³⁹ is referred to as ${\mathrm{M}}_4$ . Figure 18 shows the RUL distributions of the assessed unit at some points, and Figure 19 shows the RUL distributions of the assessed unit at 3500 h. From Figure 18, we can observe that due to the ME of laser degradation data is very small, the RUL distributions of ${\mathrm{M}}_3$ and ${\mathrm{M}}_4$ can all cover the actual RULs, and the RUL distributions of ${\mathrm{M}}_3$ and ${\mathrm{M}}_4$ are very close. In addition, from Figure 19, we can see that the RUL distributions of ${\mathrm{M}}_3$ is slightly higher than ${\mathrm{M}}_4$ , which shows that the RUL distributions of ${\mathrm{M}}_3$ is more concentrated and the uncertainty of prediction results is smaller.

Figure 18.

The estimated RULs at some points.

Figure 19.

The estimated RULs at 3500 h.

In order to significantly compare the effect of considering ME, we appropriately add a ME following normal distribution ${\epsilon }^{\prime }~N(0,0.2)$ to the original laser degradation data, and regard the original degradation data as the actual value. Its degradation path is displayed in Figure 20.

Figure 20.

The degradation path of simulated laser data.

Similarly, the failure time of black dotted lines in Figure 20 are 3450.29, 3571.71, and 3298.04 h, respectively. Figure 21 shows the RUL distributions of ${\mathrm{M}}_3$ and ${\mathrm{M}}_4$ at some points. From Figure 21, it can be observed that the RUL distributions of ${\mathrm{M}}_3$ is narrower than ${\mathrm{M}}_4$ , and the RUL distributions of ${\mathrm{M}}_3$ is more concentrated on actual RULs. This result shows that when the degradation data have ME, ${\mathrm{M}}_3$ can get higher prediction accuracy than ${\mathrm{M}}_4$ . Then, we calculate the MSEs at different points, as shown in Figure 22. The MSEs of ${\mathrm{M}}_3$ are always smaller than ${\mathrm{M}}_4$ , which further reflects the superiority of the proposed method. In summary, it is necessary to consider the influence of ME for degradation systems with ME.

Figure 21.

The estimated RULs at some points.

Figure 22.

The MSEs of RULs at some points.

Conclusions

Accurate parameters estimation is the premise to ensure the accuracy of RUL prediction. In this paper, we first analyze the natures of parameters estimation of the implicit linear Wiener degradation process. Then, based on these natures, a RUL prediction method with fusing multi-source information is proposed for the implicit linear Wiener degradation process. Finally, we use some simulation experiments and two practical cases studies to verify the natures of parameters estimation and the effectiveness of the proposed RUL prediction method. The main contributions of this paper can be summarized as follows:

• (1) Based on the unbiased parameter estimation results of the implicit linear Wiener degradation process, the variances of the estimated parameters are deduced. And the relationship between the accuracy of parameters estimation and the number of units, the length of detection time and the number of detection times are analyzed based on the implicit linear Wiener degradation process. It provides a theoretical basis regarding how to reasonably fuse multi-source information for RUL prediction.

• (2) For stochastic degradation systems with ME, a RUL prediction method with fusing multi-source information is proposed to deal with the problem that the prior degradation information is imperfect or scarce.

At present, there are many methods that can be used to solve the problem of imperfect and scarce prior degradation information. For example, the RUL prediction method based on Gaussian regression process ⁴⁷ can use a small amount of historical degradation data to predict the RUL, and the RUL prediction method based on particle filter algorithm ⁴⁸ can use the field degradation data of the assessed equipment to predict the RUL, and they can all get good prediction results. In order to further improve accuracy of RUL prediction, the comparison of these methods is worth further study in the future.