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

#### 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 Section 2, 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

It is further assumed that ; then, we havewhere is the drift coefficient, which is used to represent the product-to-product uncertainty, is a function of time and the parameter vector , is the diffusion coefficient, is the standard Brownian motion (BM) process, used to describe the temporal uncertainty of the degradation process, and is the initial degeneration state of the system. In the following, we assume .

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., , and and 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 data cannot reflect the actual degradation state of the system accurately. The relationship between them may be expressed aswhere describes the random measurement errors and . , , 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 may be expressed as

Similarly, the RUL of a system at time is given by

At first, for the degradation data without measurement errors, we have the following lemma [20].

Lemma 1. *For the degradation process in (1), if the function is continuous in , under mild assumptions [20], then the PDF of the lifetime of crossing the constant boundary is given bywhere **Based on Lemma 1, the PDF of the RUL of a system at time and degradation state iswhere and .**If one considers measurement errors, the FHT of crossing threshold may be defined aswhere and .**Similarly, the RUL of a system with measurement errors at time is expressed as follows:where follows a normal distribution [13, 14], i.e., .**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 is obtained using the law of total probability, i.e.,**To calculate the integral in equation (10) exactly, the following lemma holds [20].*

Lemma 2. *If and , , then we have**Now, based on Lemma 1 and 2, we can obtain the PDF of the FHT.*

Theorem 1. *If and , then the PDF of the FHT for a degradation model considering measurement errors may be expressed aswhere .*

*Proof. *By replacing in *E*, (6) with , and considering that , we haveNow, according to Lemma 2, the result is proved by setting , , , and . Furthermore, when the drift parameter is taken into account, we haveOn the basis of Lemma 2, Theorem 1 is proved by setting , , , and .

In the same way, we obtain the following theorem.

Theorem 2. *If and , given the current measured data , then the PDF of the RUL for the established model with measurement errors at time can be expressed aswhere .**The proof of Theorem 2 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

Before determining parameters of the model, we need to know the function . Following the approach [37], we know that the form of 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 . Indeed, the function 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 may be obtained easily as .

The specific steps to determine the function are the following. At first, the degradation data for all units are averaged at each sampling time. Afterwards, the sample average is fitted with a large class of possible functional forms, including the exponential function and power laws. The fitting function represents , where and is the measurement data of unit at the th time, and .

##### 3.2. The Parameter Estimation

Let us now address the estimation of the parameters , and in equations (2) and (3). We denote by the vector of all unknown parameters. We assume that there are independent units for each testing with denoting the measurement data for the th unit at time , . Here, we assume that the degradation data of each unit is measured in the same time sequence . According to equation (2) and (3), the degradation measurement of the unit can be expressed as , where is the drift coefficient of the th unit. The variables are i.i.d of and are i.i.d realizations of . The increment of the degradation measurement for th unit is expressed as and the increments of the degradation measurements for the test units are denoted as . According to the nature of the Wiener processes, when is fixed, , and the joint PDF of may be expressed aswhere the elements of the covariance matrix are given by

In order to simplify the inference procedure and the estimation of parameters, we use the re-parameterization and , and denotes the vector of drift coefficients. Therefore, given the measurements and , the log-likelihood function of iswhere represents a generic constant. Then, the MLE of is represented by , which is given by .

It is impossible to maximize (18) directly since the drift parameters 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 the estimate at the th step based on .

As mentioned above, . According to a Bayesian approach, the posterior distribution of can be written as follows:

Using the above derivation, it can be seen that the drift parameter of the th 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 expectation of can be written as

Setting , , and , we have

It is easily found that depends on and . Next, by substituting , , and in equation (20), we can obtain the profile likelihood function:

The genetic algorithm (GA) is used to maximize (23) and obtain the MLE of and . A more detailed overview of GA may be found in [43]. Next, the estimate of in the th step may be obtained by substituting and into equation (22). Furthermore, can be estimated as . 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 aswhere and are the estimated RUL and the actual RUL at time , respectively, and is the conditional PDF of the RUL estimated by (15).

In a situation where there are measurements, the TMSE of the lifetime can be written as

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 Case 2, 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 [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 Figure 1, 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 Section 3.1, we determine the function 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 . The calculated values of the correlation coefficients are shown in Table 1, from which an exponential is found. Assuming that the degradation measurement at is zero, we have . Then, from equations (2) and (3), we have that the estimated mean degradation path is . 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 Figure 3, 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 Figure 4.

From Figure 4, 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 Figure 5.

In Figure 5, 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.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

#### 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].

#### Data Availability

The data used to support the findings of this study have been deposited in the (Remaining Useful Life Estimation Based on a Nonlinear Diffusion Degradation Process) repository (DOI:10.1109/TR.2011.2182221).

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This research was supported by the National Natural Foundation of China (Grant no. 72071183), the Youth Program of National Natural Science Foundation of China (Grant no. 61703297), the Natural Science Foundation of Shanxi Province (Grant nos. 201801D121188 and 20210302123206) and Shanxi Scholarship Council of China (2021-135). The authors would like to express their gratitude to EditSprings (https://www.editsprings.com/) for providing expert linguistic services provided.