Remaining Useful Life Prediction of Lithium-Ion Batteries Based on the Wiener Process with Measurement Error

Abstract

Remaining useful life (RUL) prediction is central to the prognostics and health management (PHM) of lithium-ion batteries. This paper proposes a novel RUL prediction method for lithium-ion batteries based on the Wiener process with measurement error (WPME). First, we use the truncated normal distribution (TND) based modeling approach for the estimated degradation state and obtain an exact and closed-form RUL distribution by simultaneously considering the measurement uncertainty and the distribution of the estimated drift parameter. Then, the traditional maximum likelihood estimation (MLE) method for population based parameters estimation is remedied to improve the estimation efficiency. Additionally, we analyze the relationship between the classic MLE method and the combination of the Bayesian updating algorithm and the expectation maximization algorithm for the real time RUL prediction. Interestingly, it is found that the result of the combination algorithm is equal to the classic MLE method. Inspired by this observation, a heuristic algorithm for the real time parameters updating is presented. Finally, numerical examples and a case study of lithium-ion batteries are provided to substantiate the superiority of the proposed RUL prediction method.

1. Introduction

Lithium-ion batteries have been widely used in many fields, e.g., consumer electronics, electric vehicles, marine systems, aircrafts, satellites, etc., due to their high power density, low weight, long lifetime, low self-discharge rate, no memory effect and other advantages [1,2]. The demand for lithium-ion batteries demonstrates the necessity to evaluate their reliability. Failure of lithium-ion batteries could lead to performance degradation, operational impairment, and even catastrophic failure [3–5]. For example, in 2006, the National Aeronautics and Space Administration's Mars Global Surveyor stopped working due to the failure of batteries [6]. In 2013, all Boeing 787 Dreamliners were indefinitely grounded due to battery failures that occurred on two planes [7]. Therefore, monitoring the degradation process, evaluating the state of health and predicting the remaining useful life (RUL) have become increasingly important for lithium-ion batteries.

Prognostics and health management (PHM) has emerged as one of the key enablers to improve system safety, increase system operations reliability and mission availability, decrease unnecessary maintenance actions, and reduce system life-cycle costs [8,9]. As a very important step of PHM, the RUL prediction based on the condition monitoring (CM) information plays an important role in maintenance strategy selection, inspection optimization, and spare parts provision [9,10]. The RUL of lithium-ion batteries is defined as the length of time from present time to the end of useful life. Since the failure data is scarce in reality, the degradation information is often chosen to describe the health status through a degradation model [11]. For the lithium-ion batteries, the capacity induced by the charge-discharge operational cycle is a suitable feature to characterize the long-term degradation process [12,13]. Then, the issue of estimating the battery's RUL could be transformed to predict the time when its capacity crosses a predefined failure threshold.

There are two main approaches for prognostics in PHM, i.e., physics-of-failure (PoF) and data-driven. As the PoF based prognostic methods depend on the knowledge of a battery's life cycle loading condition, material properties, failure mechanisms, etc., it is difficult to apply for complicated systems with unclear physical failure mechanism. However, this is most frequently encountered in practice. In contrast, data-driven techniques extract features from performance data such as current, voltage, capacity and impedance, and thus they are less complex than the PoF-based approaches. The current research about the RUL prediction of lithium-ion batteries focuses mainly on data-driven approaches. Saha et al. [14] presented an empirical model to describe battery behavior and predicted the RUL by the particle filter (PF) algorithm. Dalal et al. [15] provided details on how the PF algorithm could be used for prognostics. To improve the accuracy of the PF algorithm, Rao-Blackwellized PF and unscented PF algorithm are used [16,17]. To fully utilize the degradation data of congeneric batteries, He et al. [18] applied the Dempster-Shafer theory to evaluate the initial model parameters. Other new prognostics-related methods and models, e.g., autoregressive model, Verhulst model, fusion prognostic algorithm, adaptive bathtub-shaped function, relevance vectors, etc., can be found in [19–25].

In the above works, the probability distribution function (PDF) of the RUL is approximated by the Monte Carlo method, which is time-consuming. To solve this problem, the Wiener process has been reported to predict the RUL of batteries [12,13]. The Wiener process can provide a good description of system's dynamic characteristic due to its non-monotonic property, infinite divisibility property and physical interpretations. As such, it has been widely used to model degradation process, such as bridge beams [26], milling machines [27], light-emitting diodes (LEDs) [28,29], micro electro mechanical systems (MEMS) [30], laser generators [31], continuous stirred tank reactors [32], and gyros [33–35]. The RUL prediction based on the Wiener process, which is a type of statistics-based data driven methods, has gained much attention in recent years [36]. In this paper, we address the issue of applying the Wiener process to predict the RUL of batteries, with an emphasis on the effect of measurement error (ME) in the measured data [12,13].

In the literature, the current research about the Wiener process with ME (WPME) focuses on the following three aspects. The first aspect is studying the influence of the ME to the RUL prediction and finding the PDF of the RUL or the lifetime. Si et al. [37] presented the analytical form of the RUL distribution for the linear WPME by considering the uncertain of the ME. Feng et al. [13] presented an analytical form of the RUL distribution for the nonlinear case. Wang et al. [38] and Tang et al. [39] proposed the RUL distributions for the WPME by incorporating the uncertainties of the ME and the degradation trend simultaneously. However, when the ME is involved, the degradation state must be estimated from the measurements with ME. Suppose that the failure threshold is w and the actual CM data at the current time is x_k which is a random value, it should be satisfied that w > x_k if the item is not failed at the current time. Thus, w − x_k has a nonzero probability to be negative. Particularly, when the degradation state approaches the failure threshold, this nonzero probability will increase and cannot be ignored. Otherwise, if we use the traditional normal distribution to model w − x_k, an underestimated probability density function of the estimated RUL will be generated and thus lead to unsuccessful prognostics. Therefore, the task of this paper is to ensure the condition that w > x_k is satisfied by introducing a truncated normal distribution (TND) based modeling approach.

The second aspect is the offline parameters estimation. The maximum likelihood estimation (MLE) method is most commonly used to estimate the fixed model parameters of the WPME. Whitmore [40] did a pioneering work for the modeling and parameters estimation of the WPME for a single item. Peng and Tseng [31] incorporated the random effects into the modeling of the WPME and presented a MLE method for a type of items. Then, this method has been applied to the nonlinear Wiener process [33] and the WPME when using first differences of the observations [12,39,41,42]. However, the drift parameter is assumed to be random for a specific item, which is not consistent with the modeling assumptions, as detailed in Section 3.

The last aspect is updating the random parameters of the WPME. The updating algorithm of the WPME is developed from the updating algorithm for the basic linear Wiener process. In this paper, we review the updating algorithm for the basic linear Wiener process and WPME together. A classic work about the updating of the random parameters is proposed by Gebraeel et al. [43], whose model established a linkage between the past and current degradation data of the congeneric items by a Bayesian mechanism. Following Gebraeel et al. [43], some related issues and many variants and applications have been studied and reported [39,44–46]. The work about this Bayesian updating mechanism for the WPME is also been studied in [41]. Another updating algorithm is combining the Bayesian updating and expectation maximization (EM) algorithm, where the Bayesian updating includes the basic Bayesian algorithm, Kalman filtering, strong tracking filter, extended Kalman filter, etc., [13,32,34,35,47]. In this combination algorithm, the Bayesian updating is used to update the random parameters of the model and the EM algorithm is used to estimate the fixed parameters. However, the relationship between the Bayesian method and the combination algorithm is still unclear.

To address the above issues, we firstly use the TND to model w − x_k and in this case an exact and analytical RUL distribution is obtained. This leads to the first contribution of the paper. Secondly, we present a two-step MLE method to estimate the parameters for satisfying the assumption that the random drift parameter is fixed for a specific item. This is another contribution of the paper since this two-step method makes the estimated variance of the drift parameter positive. Then, we find an interesting result about the relationship between the Bayesian method and the combination algorithm, and based on this result we present a heuristic algorithm for the parameters updating. This is the third contribution of this paper, which is not fully explored before. Finally, some numerical examples and a case study regarding the lithium-ion batteries are presented to verify the results derived in this paper and to illustrate the application and superiority of the proposed method.

The remainder of this paper is organized as follows: Section 2 develops the ME model and derives the RUL distribution; in Section 3, we present a two-step MLE method to estimate the fixed parameters; in Section 4, we discuss the parameters updating for the WPME and propose a heuristic algorithm; numerical examples and a case study are provided in Section 5; and Section 6 draws the main conclusions.

2. Degradation Modeling and RUL Prediction

2.1. Degradation Model

The Wiener process with a linear drift is typically used for modeling the degradation process. Let X (t) denote the degradation value at time t; the degradation process can be represented as follows:

$$X(t)=X(0)+\lambda t+{\sigma }_{B}B(t)$$

(1)

where λ is the drift parameter; σ_B is the diffusion parameter; and B(t) is the standard Brownian motion representing the stochastic dynamics of the degradation process. Without loss of generality, X(0) is assumed to be zero. As the mean of the degradation is governed by the drift parameter λ, λ is assumed as a random parameter to represent the heterogeneity among different items, and σ_B is assumed to be fixed to describe common character shared by all items. Moreover, λ is assumed to be s-independent with B(t).

Generally speaking, the degradation data obtained from routine CM are inevitably contaminated by the uncertainty during the measurement process [48,49]. If the measurement process of X(t) contains ME, only the variable with ME can be observed, which is represented as follows [40]:

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

(2)

where ε denotes the ME and is normally distributed with zero mean and standard deviation σ_ε. Additionally, ε is assumed to be s-independent with λ and B(t).

Generally, the lifetime of the system is defined by the concept of the first hitting time (FHT) of the degradation process as:

$$T=inf\{t:X(t)\ge w|X(0)\le w\}$$

(3)

where w denotes the failure threshold.

In the sense of the FHT, it is well-known that the Wiener process crossing a constant threshold w obeys an inverse Gaussian distribution [50]. Accordingly, the PDF, mean and variance of the lifetime can be obtained as follows:

$$f_T(t)=\frac{w}{\sqrt{2\pi t^3{\sigma }_B^2}}exp\left(-\frac{{(w-\lambda t)}^2}{2{\sigma }_B^{2}t}\right)$$

(4)

$$E(T)=\frac{w}{\lambda },var(T)=\frac{w{\sigma }_B^2}{{\lambda }^3}$$

(5)

2.2. Modelling the ME

For the degradation process with ME, the lifetime T_e can be defined as:

$$T_{\text{e}}=inf\left\{t:X(t)\ge w|X(0)\le w\right\}=inf\left\{t:Y(t)\ge w+\mathit{\epsilon }|X(0)\le w\right\}$$

(6)

Therefore, the lifetime T_e can be calculated by the time of {Y(t), t ≥ 0} hitting the threshold w_ε = w + ε. Given the probability distribution of the ME, such as normal distribution, the PDF of the RUL can be derived by the law of total probability [31,40]. However, the restriction that w − x_k > 0 is not satisfied in References [31,40]. As mentioned above, we use the TND to solve this problem. If Z ∼ N(μ, σ²) and Z is truncated by Z > 0, then this type of TND can be written as Z ∼ TN(μ, σ²). Suppose that Z ∼ TN(μ, σ²), the PDF and mean of Z can be written as [51]:

$$f(Z)=\frac{1}{\sqrt{2{\pi \sigma }^2}\Phi (\mu /\sigma )}exp\left(-\frac{{(z-\mu )}^2}{2{\sigma }^2}\right)$$

(7)

$$E(Z)=\frac{\sigma }{\sqrt{2\pi }\Phi (\mu /\sigma )}exp\left(-{\mu }^2/(2{\sigma }^2)\right)+\mu$$

(8)

where Φ() is the cumulative distribution function (CDF) of the standard normal distribution.

If μ ≫ σ, then we have Φ(μ/σ) ≈ 1 and exp(−μ²/σ²) ≈ 0. Thus, the TND approximately turns into the traditional normal distribution. The key issue of the RUL prediction with ME is to derive the PDF of the RUL, which is addressed in the following subsection.

2.3. Real Time RUL Prediction with ME

In this subsection, the RUL prediction method at a particular point of time t_k is proposed. Define Y_0:k = {y₀,y₁,y₂,…,y_k} as the observed degradation and X_0:k = {x₀,x₁,x₂,…,x_k} as the actual degradation at CM times t₀, t₁, t₂,…,t_k, which could be irregularly spaced. Once the actual degradation process X_0:k is available at t_k, the process {X(t),t ≥ t_k} can be transformed into [34,35]:

$$\phi (l_k)=X(l_k+t_k)-X(t_k)=\lambda (l_k)+{\sigma }_{B}B(l_k)$$

(9)

where l_k = t − t_k ≥ 0, φ(0) = 0.

Therefore, the RUL at time t_k is equal to the FHT of the process {φ(l_k), l_k ≥ 0} crossing the threshold w_k = w − x_k. Accordingly, the RUL can be defined as:

$$L_k=inf\left\{l_k:X(l_k+t_k)\ge w|{\mathbf{\text{X}}}_{0:k}\right\}=inf\left\{l_k:\phi (l_k)\ge w-x_k|{\mathbf{\text{X}}}_{0:k}\right\}$$

(10)

As a result, once X_0:k is available at t_k, the PDF of the RUL with considering the uncertainty of λ in the sense of FHT can be written as [34,35]:

$$f_{L_k|X_{0:k}}(l_k|{\mathbf{\text{X}}}_{0:k})=\frac{w-x_k}{\sqrt{2\pi {l_k}^2({\sigma }_B^2{l}_k+{\sigma }_{\lambda ,k}^2{l_k}^2)}}exp\left(-\frac{{(w-x_k-{\mu }_{\lambda ,k}{l}_k)}^2}{2{\sigma }_B^2{(\sigma }_B^{2}t+{\sigma }_{\lambda ,k}^2{l_k}^2)}\right)$$

(11)

where μ_λ,k and ${\sigma }_{\lambda ,k}^2$ are the estimated mean and variance of λ conditional on X_0:k by the Bayesian method. For more details about Equation (11), see References [34,35].

However, for the degradation process with ME, X_0:k is unobserved. Only the degradation process with ME, i.e., Y_0:k, can be directly measured. Thus, the true value of x_k is a random value with uncertainty. To derive the PDF of the RUL with ME, the following lemma is given first:

Lemma 1 If Z ∼ TN(μ,σ²), and B ∈ R, C ∈ R⁺, then:

$$E_Z\left[Z\cdot exp\left(-\frac{{(Z-B)}^2}{2C}\right)\right]=\frac{1}{\Phi \left(\mu /\sigma \right)}\frac{\sqrt{C}}{C+{\sigma }^2}\left\{\sqrt{\frac{C{\sigma }^2}{2\pi }}exp\left(-\frac{B^2{\sigma }^2+{\mu }^{2}C}{2C{\sigma }^2}\right)+\frac{B{\sigma }^2+\mu C}{\sqrt{C+{\sigma }^2}}exp\left(-\frac{{(\mu -B)}^2}{2(C+{\sigma }^2)}\right)\mathrm{\Phi }\left(\frac{B{\sigma }^2+\mu C}{\sqrt{(C+{\sigma }^2)C{\sigma }^2}}\right)\right\}$$

(12)

The proof is given in the Appendix.

As the ME $\mathit{\epsilon }~N(0,{\sigma }_{\mathit{\epsilon }}^2)$, it follows that $w-x_k~N(w-y_k,{\sigma }_{\mathit{\epsilon }}^2)$. On condition that w − x_k > 0, we have $w-x_k~TN(w-y_k,{\sigma }_{\mathit{\epsilon }}^2)$. Then, based on Lemma 1 and Equation (11), the PDF and mean of the RUL based on the updated λ can be obtained by the law of total probability, which is summarized in following theorem:

Theorem 1 Once Y_0:k is available at t_k the PDF and mean of the RUL with considering the uncertainty the estimated λ and ME can be written as:

$$\begin{array}{l}{f}_{L_k|Y_{0:k}}(l_k|{\mathbf{\text{Y}}}_{0:k})=\frac{1}{\sqrt{2\pi \mathrm{\Phi }\left({\mu }_{w,k}/{\sigma }_{\mathit{\epsilon }}\right)}D}\{\sqrt{\frac{E{\sigma }_{\mathit{\epsilon }}^2}{2\pi {l_k}^2}}exp\left(-\frac{{{\mu }_{\lambda ,k}}^2{l_k}^2{\sigma }_{\mathit{\epsilon }}^2+{{\mu }_{w,k}}^{2}E}{2E{\sigma }_w^2}\right)\hfill \\ +\frac{{\mu }_{\lambda ,k}{\sigma }_{\mathit{\epsilon }}^2{l}_k+{\mu }_{w,k}E}{\sqrt{D}{l}_k}exp\left(-\frac{{({\mu }_{w,k}-{\mu }_{\lambda ,k}{l}_k)}^2}{2D}\right)\mathrm{\Phi }\left(\frac{{\mu }_{\lambda ,k}{l}_k{\sigma }_{\mathit{\epsilon }}^2+{\mu }_{w,k}E}{\sqrt{DE{\sigma }_{\mathit{\epsilon }}^2}}\right)\}\hfill \end{array}$$

(13)

$$E(L_k|{\mathbf{\text{Y}}}_{0:k})=\left[\frac{{\sigma }_{\mathit{\epsilon }}exp\left(-{{\mu }_{w,k}}^2/{\sigma }_{\mathit{\epsilon }}^2\right)}{\sqrt{2\pi }\mathrm{\Phi }\left({\mu }_{w,k}/{\sigma }_{\mathit{\epsilon }}\right)}+{\mu }_{w,k}\right]\sqrt{\frac{2}{{\sigma }_{\lambda ,k}^2}}\Lambda \left({\mu }_{\lambda ,k}/\sqrt{2{\sigma }_{\lambda ,k}^2}\right)$$

(14)

where $D={\sigma }_B^2{l}_k+{\sigma }_{\lambda ,k}^2{l_k}^2+{\sigma }_{\mathit{\epsilon }}^2$; μ_w,k = w − y_k; $E={\sigma }_B^2{l}_k+{\sigma }_{\lambda ,k}^2{l_k}^2$; and Λ() is the Dawson function, which can be expressed as follows:

$$\Lambda (x)=\{\begin{array}{ll}{e}^{-x^2}{\int }_0^x{e}^{t^2}\text{d}t\hfill & x\ge 0\hfill \\ e^{x^2}{\int }_0^x{e}^{-t^2}\text{d}t\hfill & x>0\hfill \end{array}$$

The proof is given in the Appendix.

Similarly, if the ${\sigma }_{\mathit{\epsilon }}^2$ is set to be zero, the result of Theorem 1 reduces to the result of Equation (11). This indicates that the result of Theorem 1 is a generalization of Equation (11). If μ_w,k ≫ σ_ε, we obtain that $exp\left({-{\mu }_{w,k}}^2/{\sigma }_{\mathit{\epsilon }}^2\right)\approx 0$ and Λ(z) = 1/2z based on the approximation property of the Dawson integral function. Hence, we have E(L_k|Y_0:k) ≈ _w,k/ μ_λ,k. This result is desired since the expectation of the RUL is required in some maintenance strategies [45,46]. In the following, we develop a parameters estimation algorithm to estimate and update the parameters in Equations (13) and (14).

3. Offline Parameters Estimation Method

An offline parameters estimation method is needed to estimate the fixed parameters of the WPME, which represents population-based degradation character. The results of this estimation could provide the prior information for updating the online parameters. To represent the heterogeneity among different items, the drift parameter of the Wiener process is usually assumed as a random parameter with mean μ₀ and variance ${\sigma }_0^2$ [12,31]. Strictly speaking, the drift parameter is assumed to be random for the population but fixed for a specific item. Peng and Tseng [31] presented a MLE method to estimate the unknown parameters of the WPME. This method has been further applied for the cases of nonlinear Wiener process [33], and the WPME when using the first differences of observations to develop the sample likelihood function (SLF) [12,39,41]. For more details of this MLE method, refer to Reference [31]. However, if the drift parameter λ is assumed to be random for a specific item, the estimated variance of the drift parameter may be evaluated to be negative by the MLE method presented in Reference [31]. This phenomenon is not consistent with the modeling assumptions and the actual conditions. To solve this problem, we develop a two-step MLE method. The first step is estimating the parameters $\Theta =\{{\lambda }_1{,\lambda }_2,\dots {,\lambda }_{\text{N}}{,\sigma }_{\text{B}}^2{,\sigma }_{\mathit{\epsilon }}^2\}$, where λ_n denotes the mean degradation rate for a specific item. The second is estimating the randomness of λ, i.e., $({\mu }_{\lambda },{\sigma }_{\lambda }^2)$. In the following, we present the two-step MLE method in detail.

We use the first differences of the observations to develop the SLF, similarly to References [12,41]. It is assumed that there are N tested items, and the degradation of the nth item is measured at times t_1,n, t_2,n,…, t_mn,n, where m_n denotes the available number of degradation measurements of the n-th item, and n = 1,2,…N. Let Y_n = {ΔY(t_1,n),ΔY(t_2,n),…,ΔY(t_m,n)}′ , where ΔY(t_i,n) = Y(t_i,n) − Y(t_i−1,n). In addition, define T_n = {ΔT_1,n, ΔT_2,n,,…, ΔT_mn,n}′ with ΔT_i,n = t_i,n − t_i−1,n. Then, according to the properties of the Wiener process, Y_n follows a multivariate normal distribution with mean and covariance given by:

$${\mathit{\text{Y}}}_n={\lambda }_n{\mathit{\text{T}}}_n,{\mathbf{\Sigma }}_n={\sigma }_B^2{\mathit{\text{D}}}_n+{\sigma }_{\mathit{\epsilon }}^2{\mathit{\text{P}}}_n$$

(15)

where:

$${\mathit{\text{D}}}_n=\text{diag}(\mathrm{\Delta }{T}_{1,n},\mathrm{\Delta }{T}_{2,n},\dots ,\mathrm{\Delta }{T}_{m_n{,}_n}),{\mathit{\text{P}}}_n\left[\begin{array}{ccccc}1& -1& 0& \dots & 0\\ -1& 2& -1& \dots & 0\\ 0& -1& 2& \dots & ⋮\\ ⋮& ⋮& ⋮& \ddots & -1\\ 0& 0& \dots & -1& 2\end{array}\right]$$

Note that the value in the first row and the first column of P_n, i.e., P_n(1,1), is equal to one, which is different from the result in Reference [40]. The reason is that Reference [40] ignored the observation at the first time. Recently, Ye et al. [41,42] corrected this mistake and set P_n(1,1) as 1.

To facilitate the inference, we re-parameterize the parameters by ${\tilde{\sigma }}_{\mathit{\epsilon }}^2{=\sigma }_{\mathit{\epsilon }}^2{/\sigma }_B^2$ and ${\tilde{\Sigma }}_n={\Sigma }_n/{\sigma }_B^2$. Then, the log-likelihood function of $\Theta =\{{\lambda }_1{,\lambda }_2,\dots {,\lambda }_{\text{N}}{,\sigma }_B^2,{\sigma }_{\mathit{\epsilon }}^2\}$ can be written as:

$$lnL(\Theta |\mathbf{\text{Y}})=-\frac{ln2\pi }{2}\sum _{n=1}^N{m}_n-\frac{1}{2}ln{\sigma }_B^2\sum _{n=1}^N{m}_n-\frac{1}{2}\sum _{n=1}^{N}ln\left|{\tilde{\Sigma }}_n\right|-\frac{1}{2{\sigma }_B^2}\sum _{n=1}^N{\left({\mathit{\text{Y}}}_n-{\lambda }_n{\mathit{\text{T}}}_n\right)}^{\prime }{{\tilde{\sum }}_n}^{-1}\left({\mathit{\text{Y}}}_n-{\lambda }_n{\mathit{\text{T}}}_n\right)$$

(16)

Taking the first partial derivatives of ln L(Θ|Y) with respect to $({\lambda }_1{,\lambda }_2,\dots {,\lambda }_{\text{N}}{,\sigma }_B^2)$ yields:

$$\frac{\partial }{\partial {\lambda }_n}lnL(\mathbf{\Theta }|\mathbf{\text{Y}})=-\frac{1}{{\sigma }_B^2}\left(-{{\mathit{\text{T}}}_n}^{\prime }{{\tilde{\sum }}_n}^{-1}{\mathit{\text{Y}}}_n+{\lambda }_n{{\mathit{\text{T}}}_n}^{\prime }{{\tilde{\sum }}_n}^{-1}{\mathit{\text{T}}}_n\right)$$

(17)

and:

$$\frac{\partial }{\partial {\sigma }_B^2}lnL(\mathbf{\Theta }|\mathbf{\text{Y}})=-\frac{1}{2{\sigma }_B^2}\sum _{n=1}^N{m}_n+\frac{1}{2{({\sigma }_B^2)}^2}\sum _{n=1}^N{\left({\mathit{\text{Y}}}_n-{\lambda }_n{\mathit{\text{T}}}_n\right)}^{\prime }{{\tilde{\sum }}_n}^{-1}\left({\mathit{\text{Y}}}_n-{\lambda }_n{T}_n\right)$$

(18)

Then, by setting these derivatives with respect to $({\lambda }_1{,\lambda }_2,\dots {,\lambda }_{\text{N}}{,\sigma }_B^2)$ to zeros, the results of the MLE for λ₁,λ₂,…,λ_N and ${\sigma }_B^2$ can be written as:

$${\widehat{\lambda }}_n=\frac{{{\mathit{\text{T}}}_n}^{\prime }{{\widetilde{\text{∑}}}_n}^{-1}{\mathit{\text{Y}}}_n}{{{\mathit{\text{T}}}_n}^{\prime }{{\tilde{\mathbf{\sum }}}_n}^{-1}{\mathit{\text{T}}}_n},$$

(19)

$${\widehat{\sigma }}_B^2=\frac{\text{∑}_{n=1}^N{\left({\mathit{\text{Y}}}_n-{\lambda }_n{\mathit{\text{T}}}_n\right)}^{\prime }{{\tilde{\sum }}_n}^{-1}\left({\mathit{\text{Y}}}_n-{\lambda }_n{\mathit{\text{T}}}_n\right)}{\text{∑}_{n=1}^N{m}_n}.$$

(20)

Substituting Equations (19) and (20) into Equation (16), we obtain the profile likelihood function for σ͂² in terms of the estimated $({\lambda }_1{,\lambda }_2,\dots {,\lambda }_{\text{N}}{,\sigma }_{\text{B}}^2)$ as follows:

$$lnL({\tilde{\sigma }}^2|\mathbf{Y})=-\frac{ln2\pi }{2}\sum _{n=1}^N{m}_n-\frac{1}{2}\sum _{n=1}^N{m}_{n}ln{\widehat{\sigma }}_B^2-\frac{1}{2}\sum _{n=1}^{N}ln\left|{\tilde{\Sigma }}_n\right|-\frac{1}{2}\sum _{n=1}^N{m}_n$$

(21)

The MLE for ${\tilde{\sigma }}_{\mathit{\epsilon }}^2$ can be obtained by maximizing the profile log-likelihood function in Equation (21) through a one-dimensional search. In this paper, we use the MATLAB function “FMINSEARCH” to find the estimates of σ͂². Then, the estimates of $({\lambda }_1{,\lambda }_2,\dots {,\lambda }_{\text{N}}{,\sigma }_{\text{B}}^2)$ can be obtained by substituting the estimates of ${\tilde{\sigma }}_{\mathit{\epsilon }}^2$ into Equations (19) and (20). Finally, the estimate of ${(\mu }_{\lambda }{,\sigma }_{\lambda }^2)$ can be calculated as follows:

$${\widehat{\mu }}_0=\frac{1}{N}\sum _{n=1}^N{\lambda }_n,{\widehat{\sigma }}_0^2=\frac{1}{N}\sum _{n=1}^N{({\lambda }_n-{\widehat{\mu }}_0)}^2$$

(22)

Remark 1. In the MLE method presented by Peng and Tseng [31], Y_n follows a multivariate normal distribution with the mean μ₀ΔY_n and covariance ${\sigma }_0^2{\mathit{\text{Y}}}_n{\mathit{\text{Y}}}_n^{\prime }+{\sigma }_B^2{\mathit{\text{D}}}_n+{\sigma }_{\mathit{\epsilon }}^2{\mathit{\text{P}}}_n$. It can be observed that the main diagonal of Σ_n is always positive, i.e., ${\sigma }_{\lambda }^2\mathrm{\Delta }{T_{i,n}}^2+{\sigma }_B^2\mathrm{\Delta }{T}_{i,n}+{\sigma }_{\mathit{\epsilon }}^2>0$, or ${\sigma }_{\lambda }^2\mathrm{\Delta }{T_{i,n}}^2+{\sigma }_B^2\mathrm{\Delta }{T}_{i,n}+2{\sigma }_{\mathit{\epsilon }}^2>0$. However, at the maximizing process of the lnL(Θ|Y), it cannot be ensured that ${\sigma }_{\lambda }^2>0$ on condition that ${\sigma }_{\lambda }^2\mathrm{\Delta }{T_{i,n}}^2+{\sigma }_B^2\mathrm{\Delta }{T}_{i,n}+2{\sigma }_{\mathit{\epsilon }}^2$ is the sum of ${\sigma }_{\lambda }^2\mathrm{\Delta }{T_{i,n}}^2$, ${\sigma }_B^2\mathrm{\Delta }{T}_{i,n}$ and $2{\sigma }_{\mathit{\epsilon }}^2$. If ${\sigma }_{\lambda }^2$ is not restricted that ${\sigma }_{\lambda }^2>0$, it could be negative when searching the best Θ to maximize ln L(Θ|Y) by the MLE method presented in Reference [31].

Remark 2 Whitmore [40] presented two ways to develop the SLF of the WPME for a single item, i.e., the observations at each point (referred to as SLF₁) and the first differences of the observations (referred to as SLF₂). SLF₁ has been used in References [31,33], and SLF₂ in References [12,39,41]. However, to our knowledge, the issue that which one should be selected for the MLE method has not been studied before. Here, we give a simple numerical example to illustrate this issue. For a single item, let (0, 0.9, 1.6, 4.7, 4.3, 5.6, 5.4) be the observed data at time (0, 0.8, 2, 4.2, 5, 7.5, 8.9); then the estimated log-likelihood function (log-LF), λ, ${\sigma }_{\text{B}}^2$, and ${\sigma }_{\mathit{\epsilon }}^2$ based on SLF₁ are 7.5002, 0.63424, 0.32989, and 0.16090, respectively. Interestingly, the estimated results based on SLF₂ equals to those based on SLF₁. This indicates that the SLFs developed by these two ways could derive the same results. This phenomenon also occurs for the nonlinear Wiener process and the Wiener process with random effects. Note that P_n(1,1) is equal to 1 in the numerical example, it verifies the remedy of the variance matrix by Ye et al. [41] to some extent. Additionally, since the covariance matrix in the SLF₂ is the symmetric tridiagonal matrix, it is easy to calculate its inverse [52,53]. Thus, the SLF₂, i.e., by the first differences of the observations, is suggested.

4. Online Parameter Updating

Online parameter updating is used to make the estimation adapt to the item's individual characteristic and reduce the uncertainty of the estimation. As mentioned above, there are two commonly used updating methods, i.e., the Bayesian method and the combination algorithm. From a comparative perspective, the advantage of Bayesian method is that the prior information and the in situ degradation data could be reasonably incorporated. However, it is sensitive to the prior information [34]. In contrast, the merit of the combination algorithm is its robustness over the selection of these prior parameters. Essentially, the estimation by the Bayesian method could partially reflect the population information, while the combination algorithm is more adaptable to the individual information of a specific item. In this section, we attempt to find the relationship between the Bayesian method and the combination algorithm, and then present a heuristic parameters updating algorithm. We first present the combination algorithm for the WPME.

4.1. The Combination Algorithm for Parameters Updating

Parameters updating via newly observed CM data is an important part for the real time RUL prediction. To estimate the PDF of the RUL based on the in situ CM data Y_0:k, Si et al. [34] proposed a general process for the estimation of the unknown parameters through the combination of Bayesian updating and EM algorithm. Inspired by this algorithm, we use both Bayesian updating and EM algorithm to estimate the unknown parameters for the linear WPME. However, the result of the EM algorithm can be directly derived and an interesting result is obtained.

For the combination algorithm, the estimate of the unknown parameters includes the following steps:

Step 1

Determine the prior information for the unknown parameters.

First, to update the PDF of the RUL, the unknown parameters need to be estimated by the history degradation data, lifetime data or accelerated degradation data. The initial parameters are estimated by the MLE method presented in the above section. Additionally, in order to simplify the updating process, it is assumed that the prior distribution of λ₀ follows $N({\mu }_0{,\sigma }_0^2)$. Consequently, such prior distribution falls into the conjugate family of sampling distribution p(Y_0:k|λ), and the posterior estimate of λ conditional on Y_0:k is still normal, that is, $\lambda |{\mathbf{\text{Y}}}_{0:k}~N({\mu }_{\lambda ,k},{\sigma }_{\lambda ,k}^2)$.

Step 2

Update the posterior distribution of the random parameter.

Given that $\lambda ~N({\mu }_0,{\sigma }_0^2)$, based on the Bayesian theorem, the posterior distribution of λ can be obtained as follows:

$$\begin{array}{l}p\left(\lambda |{\mathit{Y}}_{0:k}\right)\propto p({\mathit{Y}}_{0:k}|\lambda )p(\lambda )\hfill \\ \propto exp\left[-\frac{1}{2{\sigma }_B^2}{\left(\mathrm{\Delta }{\mathit{Y}}_k-\lambda \mathrm{\Delta }{\mathit{t}}_k\right)}^{\prime }\Psi {(\text{f})}^{-1}\left(\mathrm{\Delta }{\mathit{Y}}_k-\lambda \mathrm{\Delta }{\mathit{t}}_k\right)\right]exp\left(-\frac{{(\lambda -{\mu }_0)}^2}{2{\sigma }_0^2}\right)\hfill \\ \propto exp\left[-\frac{1}{2}{\lambda }^2\left(\mathrm{\Delta }{{\mathit{t}}^{\prime }}_k\Psi {(\text{f})}^{-1}\mathrm{\Delta }{\mathit{t}}_k/{\sigma }_B^2+1/{\sigma }_0^2\right)+\lambda \left(\mathrm{\Delta }{{\mathit{Y}}^{\prime }}_k\Psi {(\text{f})}^{-1}\mathrm{\Delta }{\mathit{t}}_k/{\sigma }_B^2+{\mu }_0/{\sigma }_0^2\right)\right]\hfill \\ \propto exp\left(-\frac{{(\lambda -{\mu }_{\lambda ,k})}^2}{2{\sigma }_{\lambda ,k}^2}\right)\hfill \end{array}$$

(23)

where Δy_k = {y₁−y₀,y₂−y₁,…,y_k−y_k−1}, and Δt_k = {t₁−t₀,t₂−t₁,…,t_k−t_k−1}. Due to the property of the normal distribution of λ | Y_0:k, we have:

$$\lambda |{\mathbf{\text{Y}}}_{0:k}\sim N({\mu }_{\lambda ,k},{\sigma }_{\lambda ,k}^2)$$

(24)

with:

$$\begin{array}{l}{\mu }_{\lambda ,k}=\frac{\mathrm{\Delta }{{\mathit{y}}^{\prime }}_k\Psi {(\varphi )}^{-1}\mathrm{\Delta }{\mathit{t}}_k\cdot {\sigma }_0^2+{\mu }_0{\sigma }_B^2}{\mathrm{\Delta }{{\mathit{t}}^{\prime }}_k\Psi {(\varphi )}^{-1}\mathrm{\Delta }{\mathit{t}}_k\cdot {\sigma }_0^2+{\sigma }_B^2}\hfill \\ {\sigma }_{\lambda ,k}^2=\frac{{\sigma }_0^2{\sigma }_B^2}{\mathrm{\Delta }{{\mathit{t}}^{\prime }}_k\Psi {(\varphi )}^{-1}\mathrm{\Delta }{\mathit{t}}_k\cdot {\sigma }_0^2+{\sigma }_B^2}\hfill \end{array}$$

(25)

Step 3

Estimate the unknown parameters via the EM algorithm.

For simplicity, the linear WPME is referred to as Model 0, and the basic linear Wiener process as Model 1. The unknown parameters of Model 0 are ${\sigma }_B^2$, ${\sigma }_{\mathit{\epsilon }}^2$, and the parameters μ₀ and ${\sigma }_0^2$ in prior distribution p(λ), denoted by $\mathrm{\Phi }=\{{\mu }_0,{\sigma }_0^2,{\sigma }_B^2,{\sigma }_{\mathit{\epsilon }}^2\}$. Accordingly, Model 1 includes three unknown parameters, i.e., ${\mathrm{\Phi }}^{\prime }=\{{{\mu }^{\prime }}_0{,{\sigma }^{\prime }}_0^2,{{\sigma }^{\prime }}_B^2\}$. For the real time RUL prediction, once new observation y_k is available, the unknown parameters can be calculated by the MLE method as follows:

$${\widehat{\mathrm{\Phi }}}_k=arg\underset{\mathrm{\Phi }}{max}lnp({\mathbf{\text{Y}}}_{0:k}|\mathrm{\Phi })$$

(26)

where p(Y_0:k|Φ) is the joint PDF of the observed degradation data Y_0:k.

However, it is difficult to maximize the log-likelihood function due to the random effect and unobservability of λ. Generally, the EM algorithm provides a possible way to solve this problem [54]. The fundamental principle of the EM algorithm is to manipulate the relationship between p(Y_0:k|Φ) and p(Y_0:k,λ|Φ) via the Bayesian theorem so that the estimating of Φ can be achieved by two steps:E-step and M-step:

• E-step:

  $$\text{Calculate}L\left(\mathrm{\Phi }|{\widehat{\mathrm{\Phi }}}_k^{(i)}\right)=E_{\lambda |Y_{0:k}.{\widehat{\mathrm{\Phi }}}_k^{(i)}}\left\{lnp({\mathbf{\text{Y}}}_{0:k},\lambda |\mathrm{\Phi })\right\}$$

  (27)
  where ${\widehat{\mathrm{\Phi }}}_k^{(i)}$ denotes the estimated parameters in the i-th step of the EM algorithm conditional on Y_0:k.

• M-step:

  $${\widehat{\mathrm{\Phi }}}_k^{(i+1)}=arg\underset{\mathrm{\Phi }}{max}L\left(\mathrm{\Phi }|{\widehat{\mathrm{\Phi }}}_k^{(i)}\right)$$

  (28)

Then, the E-step and M-step are iterated multiple times until a criterion of convergence is satisfied. The commonly used criterion of convergence is that the difference between ${\widehat{\mathrm{\Phi }}}_k^{(i)}$ and ${\widehat{\mathrm{\Phi }}}_k^{(i+1)}$ falls below a pre-defined threshold. For more details of the convergence properties of the EM algorithm, see [54].

For Model 0, we first evaluate the complete log-likelihood function as follows:

$$\begin{array}{ll}lnp({\mathbf{Y}}_{0:k},\lambda |\mathrm{\Phi })\hfill & =lnp({\mathbf{Y}}_{0:k}|\lambda ,\mathrm{\Phi })+lnp(\lambda |\mathrm{\Phi })\hfill \\ \hfill & =-\frac{k+1}{2}ln2\pi -\frac{k}{2}ln{\sigma }_B^2-\frac{1}{2}ln\left|\Psi (\varphi )\right|-\frac{1}{2{\sigma }_B^2}{\left(\mathrm{\Delta }{\text{y}}_k-\lambda \mathrm{\Delta }{\text{t}}_k\right)}^{\prime }\Psi {(\varphi )}^{-1}\left(\mathrm{\Delta }{\mathit{y}}_k-\lambda \mathrm{\Delta }{\mathit{t}}_k\right)-\frac{1}{2}ln{\sigma }_{0,k}^2-\frac{1}{2{\sigma }_{0,k}^2}{(\lambda -{\mu }_{0,k})}^2\hfill \end{array}$$

(29)

Given ${\widehat{\mathrm{\Phi }}}_k^{(i)}=\{{\widehat{\mu }}_{0,k}^{(i)},{\widehat{\sigma }}_{0,k}^{2(i)},{\widehat{\sigma }}_{B,k}^{2(i)},{\widehat{\varphi }}_k^{(i)}\}$ as the estimate in the i-th step based on Y_0:k. The expectation of ln p(Y_0:k,λ|Φ) can be calculated as follows:

$$\begin{array}{ll}{L}_k\left(\mathrm{\Phi }|{{\widehat{\mathrm{\Phi }}}_k}^{(i)}\right)\hfill & =E_{\lambda |{\mathit{y}}_{0:k}.{{\widehat{\mathrm{\Phi }}}_k}^{(i)}}\left\{lnp({\mathit{y}}_{0:k},\lambda |\mathrm{\Phi })\right\}\hfill \\ \hfill & =-\frac{k+1}{2}ln2\pi -\frac{1}{2}ln\left|\Psi (\varphi )\right|-\frac{k}{2}ln{\sigma }_{B,k}^{(j)2}-\frac{1}{2}ln{\sigma }_{0,k}^{(j)2}-\frac{{{\mu }_{\lambda ,k}}^2+{\sigma }_{\lambda ,k}^2-2{\mu }_{\lambda ,k}{\mu }_{0,k}^{(j)}+{\mu }_{0,k}^{(j)2}}{2{\sigma }_{0,k}^{(j)2}}\hfill \\ \hfill & -\frac{1}{2{\sigma }_{B,k}^{(j)2}}\left[{(\mathrm{\Delta }{\mathit{y}}_k-\mathrm{\Delta }{\mathit{t}}_k)}^{\prime }\Psi {({\widehat{\varphi }}_k^{(i)})}^{-1}(\mathrm{\Delta }{\mathit{y}}_k-\mathrm{\Delta }{\mathit{t}}_k)+\mathrm{\Delta }{{\mathit{t}}^{\prime }}_k\Psi {({\widehat{\varphi }}_k^{(i)})}^{-1}\mathrm{\Delta }{\mathit{t}}_k{\sigma }_{\lambda ,k}^2\right]\hfill \end{array}$$

(30)

After deriving $L_k\left(\Phi |{\widehat{\mathrm{\Phi }}}_k^{(i)}\right)$, the results of the estimated parameters in the (i + 1)-th step can be summarized in the following theorems:

Theorem 2.

(1)

${\widehat{\mathrm{\Phi }}}_k^{(i+1)}$, by maximizing $L_k\left(\mathrm{\Phi }|{\widehat{\mathrm{\Phi }}}_k^{(i)}\right)$, is given by:

$$\begin{array}{l}{\sigma }_{B,k}^{2(i+1)}=\frac{1}{k}{(\mathrm{\Delta }{\mathit{y}}_k-\mathrm{\Delta }{\mathit{t}}_k)}^{\prime }\Psi {({\widehat{\varphi }}_k^{(i+1)})}^{-1}(\mathrm{\Delta }{\mathit{y}}_k-\mathrm{\Delta }{\mathit{t}}_k)+\mathrm{\Delta }{{\mathit{t}}^{\prime }}_k\Psi {({\widehat{\varphi }}_k^{(i+1)})}^{-1}\mathrm{\Delta }{\mathit{t}}_k{\sigma }_{\lambda ,k}^2\hfill \\ {\mu }_{0,k}^{(i+1)}={\mu }_{\lambda ,k}\hfill \\ {\sigma }_{0,k}^{2(i+1)}={\sigma }_{\lambda ,k}^2\hfill \\ {\widehat{\varphi }}_k^{(i+1)}=arg\underset{\varphi }{max}{L}_k\left(\varphi |{\widehat{\mathrm{\Phi }}}_k^{(i)},{\sigma }_{B,k}^{2(i+1)},{\mu }_{0,k}^{(i+1)},{\sigma }_{0,k}^{2(i+1)}\right)\hfill \end{array}$$

(31)

where:

$$\begin{array}{l}{L}_{\text{k}}\left(\varphi |{\widehat{\mathrm{\Phi }}}_k^{(i)},{\sigma }_{B,k}^{2(i+1)},{\mu }_{0,k}^{(i+1)},{\sigma }_{0,k}^{2(i+1)}\right)\hfill \\ =C-\frac{1}{2}ln\left|\Psi ({\widehat{\varphi }}_k^{(i)})\right|-\frac{k}{2}ln\left[\frac{1}{k}{(\mathrm{\Delta }{\mathit{\text{y}}}_k-{\mu }_{\lambda ,k}\mathrm{\Delta }{\text{t}}_k)}^{\prime }\Psi {({\widehat{\varphi }}_k^{(i)})}^{-1}(\mathrm{\Delta }{\mathit{y}}_k-{\mu }_{\lambda ,k}\mathrm{\Delta }{\mathit{t}}_k)+\mathrm{\Delta }{{\mathit{t}}^{\prime }}_k\Psi {({\widehat{\varphi }}_k^{(i)})}^{-1}\mathrm{\Delta }{\mathit{t}}_k{\sigma }_{\lambda ,k}^2\right]\hfill \end{array}$$

and C denotes a generic constant that may change from line to line throughout the paper.

(2)

${\widehat{\mathrm{\Phi }}}_k^{(i+1)}$ is uniquely determined and located at the maximum.

The proof is given in Appendix. From Theorem 2, it can be observed that the M-step can be solved analytically and each iteration of the EM algorithm can be performed with only a single computation. To further simplify the computation, we propose the following Theorem.

Theorem 3.

(1)

The iteration result of the combination algorithm is:

$$\begin{array}{l}{\sigma }_{B,k}^2=\frac{1}{k}{(\mathrm{\Delta }{\mathit{y}}_k-\mathrm{\Delta }{\mathit{\text{t}}}_k)}^{\prime }\Psi {({\widehat{\varphi }}_k)}^{-1}(\mathrm{\Delta }{\mathit{y}}_k-\mathrm{\Delta }{\mathit{t}}_k)\hfill \\ {\mu }_{0,k}=\frac{\mathrm{\Delta }{{\mathit{y}}^{\prime }}_k\Psi {(\varphi )}^{-1}\mathrm{\Delta }{\mathit{t}}_k}{\mathrm{\Delta }{{\mathit{t}}^{\prime }}_k\Psi {(\varphi )}^{-1}\mathrm{\Delta }{\mathit{t}}_k}\hfill \\ {\sigma }_{0,k}^2=0\hfill \\ {\widehat{\varphi }}_k=arg\underset{\varphi }{max}{L}_k\left(\varphi |{\sigma }_{B,k}^2,{\mu }_{0,k},{\sigma }_{0,k}^2\right)\hfill \end{array}$$

(32)

where:

$$L_k\left(\varphi |{\sigma }_{B,k}^2,{\mu }_{0,k},{\sigma }_{0,k}^2\right)=C+-\frac{1}{2}ln\left|\Psi ({\varphi }_k)\right|-\frac{k}{2}ln\left[\frac{1}{k}{(\mathrm{\Delta }{\mathit{y}}_k-\mathrm{\Delta }{\mathit{t}}_k)}^{\prime }\Psi {({\widehat{\varphi }}_k)}^{-1}(\mathrm{\Delta }{\mathit{y}}_k-\mathrm{\Delta }{\mathit{t}}_k)\right]$$

(2)

The result of the EM algorithm is equal to the estimate via the MLE method by assuming a fixed λ.

The proof is given in Appendix.

Based on Theorem 3, we can obtain similar results for Model 1 as follows:

Corollary 1

(1)

The result of the combination algorithm for Model 1 is:

$$\begin{array}{l}{\sigma }_{B,k}^2=\frac{1}{k}\sum _{i=1}^k{(\mathrm{\Delta }{y}_i-{\mu }_{0,k}\mathrm{\Delta }{t}_i)}^2\hfill \\ {\mu }_{0,k}=y_k/t_k\hfill \\ {\sigma }_{0,k}^2=0\hfill \end{array}$$

(33)

(2)

It is equal to the estimate via the MLE method by assuming a fixed λ.

The proof of Corollary 1 can be easily obtained by setting ${\sigma }_{\mathit{\epsilon }}^2=0$ in the proof of Theorem 3.

The reason why the EM algorithm obtains the same results with the MLE method can be explained as follows. The EM algorithm is used to find the maximum likelihood estimate of $\mathrm{\Phi }=\{{\mu }_0{,\sigma }_0^2{,\sigma }_B^2,{\sigma }_{\mathit{\epsilon }}^2\}$, where μ₀ and ${\sigma }_0^2$ are the mean and variance of λ represent the heterogeneity among different items. However, λ is a fixed value for a specific item, therefore, the maximum likelihood estimate of ${\sigma }_0^2$ by the EM algorithm should be zero. And thus the EM algorithm turns into the traditional MLE method.

The second result of Theorem 3 indicates that the EM algorithm can overcome the impact of improper prior information on the RUL prediction. However, it completely gets rid of the impact of prior information and only depends on the in situ observed data. Therefore, how to reasonably integrate the prior information and in situ information is an important issue in the real time RUL prediction. The following subsection attempts to address this issue.

4.2. A Heuristic Parameter Updating Method

Note that the convergent result of the combination algorithm equals to maximum likelihood estimate. Consequently, if no iteration is performed, it is the traditional Bayesian updating method. As the number of iterations increases, the estimation depends more on the in situ information, i.e., the individual characteristic. To integrate the prior information and in situ information, the following two basic guidelines should be followed. The first is that the less confident on the prior information, the more iterations. The second is that the more in situ degradation data are observed, the more iteration times are performed.

Based on the above two guidelines, we present a heuristic parameter updating algorithm as shown in Figure 1.

The heuristic algorithm includes three main steps as follows:

Step 1

Calculate the priori information by the estimation method presented in Section 3, i.e., $\mathrm{\Phi }=\{{\mu }_0{,\sigma }_0^2{,\sigma }_B^2,{\sigma }_{\mathit{\epsilon }}^2\}$;

Step 2

According to the confidence of the priori information, determine the length of iteration interval;

Step 3

Suppose that there are k degradation data observed. If (k −1)/L_it ≤ i < k/L_it, the EM algorithm presented in Section 4.1 is performed i times of the iteration, where L_it is the iteration interval length determined in Step 2. Then, the RUL is predicted based on the updated parameters.

As the drift parameter λ is assumed as a random parameter to represent the heterogeneity among different items, and ${\{\sigma }_B^2,{\sigma }_{\mathit{\epsilon }}^2\}$ are assumed to be fixed which are common to all items, we only update the drift parameter by the heuristic algorithm in this paper. As more degradation data are observed, the iteration times of the EM algorithm increase and thus the estimation depends more on the in situ information. This is consistent with the second basic guideline presented above. Additionally, it can be observed that the key issue of the heuristic method is to determine the iteration interval length via the prior information. In this paper, the iteration interval length is determined by the historical experience or expert's information. The selection of L_it should base on the following principle. The more confidence there is on the prior information, the smaller L_it is.

5. Experimental Studies

In this section, we provide several numerical examples to analyse the effect without considering the restriction that w − x_k > 0, to compare the performance of the fixed MLE method with the existing method in the literature, and to verify the interesting consequence derived in this paper. Then, a practical case study of the lithium-ion batteries is illustrated to demonstrate the application of the developed method for the real time RUL prediction.

5.1. Numerical Examples

First, we provide a numerical example to show the effect without considering the restriction that w − x_k > 0. For illustrative purposes, the parameters in the degradation process are assumed as w = 10, y_k = 9 and Φ = {1,0.01,0.09,2} . Then, the corresponding PDF and CDF of the lifetime are illustrated in Figure 2, where ND represents the traditional normal distribution. Obviously, the PDF and CDF by the traditional normal distribution are smaller than those by the TND.

The reason why the traditional normal distribution obtains smaller PDF and CDF can be explained as follows. When the traditional normal distribution is used, the PDF and CDF of the RUL can be calculated by the law of total probability as follows:

$$\begin{array}{c}\begin{array}{ll}{f}_{T_{\mathit{\epsilon }}}(t)\hfill & ={\int }_{-\infty }^{+\infty }{f}_{T|\mathit{\epsilon }}(t)f(w_k)\text{d}{w}_k\hfill \\ \hfill & ={\int }_{-\infty }^0{f}_{T|\mathit{\epsilon }}(t)f(w_k)\text{d}{w}_k+{\int }_0^{+\infty }{f}_{T|\mathit{\epsilon }}(t)f(w_k)\text{d}{w}_k\hfill \\ \hfill & =f_{T_{\mathit{\epsilon }1}}(t)+f_{T_{\mathit{\epsilon }2}}(t)\hfill \end{array}\\ =f_{T_{\mathit{\epsilon }1}}(t)+f_{T_{\mathit{\epsilon }2}}(t)\end{array}$$

(34)

$$\begin{array}{ll}{F}_{T_{\mathit{\epsilon }}}(t)\hfill & ={\int }_{-\infty }^{+\infty }{F}_{T|\mathit{\epsilon }}(t)f(w_k)\text{d}{w}_k\hfill \\ \hfill & ={\int }_{-\infty }^0{F}_{T|\mathit{\epsilon }}(t)f(w_k)\text{d}{w}_k+{\int }_0^{+\infty }{F}_{T|\mathit{\epsilon }}(t)f(w_k)\text{d}{w}_k\hfill \\ \hfill & =F_{T_{\mathit{\epsilon }1}}(t)+F_{T_{\mathit{\epsilon }2}}(t)\hfill \end{array}$$

(35)

where w_k = w − x_k. It can be observed that f_Tε1(t) and F_Tε1(t) are considered when the traditional normal distribution is used to model w − x_k. However, f_Tε1(t) and F_Tε1(t) will be negative since the integral is performed on condition that w − x_k < 0. This is why the method via the traditional normal distribution derives the smaller PDF and CDF. The underestimated PDF could lead to delayed-maintenance, and then result in failure of the system. However, the TND is used to satisfy the condition that w − x_k < 0. Then, f_Tε1(t) and F_Tε1(t) are not considered for the RUL estimation. And thus the TND could provide a probability distribution with total probability one.

Then, five degradation paths are generated to demonstrate the validity of the proposed offline MLE method. For illustrative purposes, the parameters for the degradation process are assumed as Φ = {3×10⁻³,4×10⁻⁸, 1×10⁻⁵,2×10⁻⁵} . The degradation paths are simulated by the Euler approximation method with measurement frequency of 2.5 h, as shown in Figure 3.

For simplicity, the MLE method proposed in this paper is referred to as M_A, the method presented by Reference [31] as M_B. Then, the estimation results of the parameters and the log-likelihood function value are shown in

Table 1. Comparisons of the two maximum likelihood estimation (MLE) methods with simulated degradation data.

Methods	μλ	${\sigma }_{\lambda }^2$	${\sigma }_B^2$	${\sigma }_{\mathit{\epsilon }}^2$	Log-LF
Real value	3 × 10−3	4 × 10−8	1 × 10−5	2 × 10−5	-
MA	2.8735 × 10−3	1.6847 × 10−7	1.0467 × 10−5	1.8566 × 10−5	416.48
MB	2.8750 × 10−3	−6.9337 × 10−8	1.3679 × 10−5	1.4693 × 10−5	414.50
MB (with restriction)	2.8747 × 10−3	0	1.2749 × 10−5	1.5704 × 10−5	414.39

for comparison. From Table 1, it can be observed that the estimate of ${\sigma }_{\lambda }^2$ by M_B is negative. This demonstrates the better performance of the presented model. If the ${\sigma }_{\lambda }^2$ is restricted as ${\sigma }_{\lambda }^2>0$ in the process of maximizing the log-likelihood function, the estimation results are given in the last line of Table 1. Compared with M_B, our method produces closer estimations of the unknown parameters to the real results and has better fit in terms of the log-LF.

To demonstrate the interesting consequence derived in this paper, we use the third simulated degradation paths to illustrate the iteration process of the combination algorithm of the Bayesian updating and EM algorithm through estimating the unknown parameters of a specific item. Using the MLE method in Reference [40], the estimated results are Φ̂ |_MLE = {3.3279×10⁻³,0,1.3245×10⁻⁶,2.3885×10⁻⁵} . The real and random prior parameters are respectively used for the parameters updating, where the random prior parameters are chosen as Φ = {1×10⁻³,9×10⁻⁸,3×10⁻⁵,1×10⁻⁵}. The iterative process of 100 times is illustrated in Figure 4. It is shown that the iteration based on the right prior parameters converges faster than that under the random prior parameters. However, they all converge to the results by using the MLE directly.

5.2. A Practical Case Study

In this experiment, we use the data collected from the National Aeronautics and Space Administration (NASA) Ames Prognostics Center of Excellence to demonstrate the effectiveness of our algorithm [55]. The degradation data includes four Li-ion batteries run through three different operational profiles (charge, discharge and impedance) at room temperature. The use of the lithium-ion battery is a process of repeatedly charging and discharging. Based on the analysis of the available performance measures of the lithium-ion batteries presented by Jin et al. [12], the capacity can be used to characterize the long-term degradation process induced by the charge-discharge operational cycle. The cycle life of the battery is defined as the number of times a battery can be recharged before its capacity has faded beyond acceptable limits (20%∼30% of the rated capacity). The degradation data from NASA included two time scales, i.e., calendar time and cycle time [12]. The degradation data based on the two time scales are shown in Figure 5.

It can be observed that some rest time exists in the experiment, which leads to the recovery of the battery. This phenomenon, called relaxation effect, increases the available capacity for the next cycle [12]. Due to the uncertainty of the rest time, the RUL could not be accurately predicted. Therefore, the relaxation effect should be extracted to study the RUL prediction of the cycle time. Jin et al. [12] used an exponential function of the calendar time presented in Reference [14] to extract the relaxation effect. For more details about how to extract the relaxation effect, see Reference [14]. In this paper, we also extract the relaxation effect by the method presented in Reference [14]. The transformed data by extracting the relaxation effect is illustrated in Figure 6, where the data is collected until failure. When the full charge capacity is reduced to below 70% of its rated value (1.4 A h), it is considered as the end of life and the battery needs to be replaced.

In the following, we compare our method with Bayesian algorithm [43]. For simplicity, the method proposed in this paper is referred to as M₁, the method based on the Bayesian updating algorithm as M₂, and the method based on the MLE method as M₃. To evaluate the proposed method, B0006 is chosen to compare these methods at different CM times and the data of other items are used to evaluate the prior parameters, the results of which are Φ = {−6.30×10⁻³,3.92×10⁻⁶,4.51×10⁻⁵,2.55×10⁻⁵} . For comparison, the inappropriate prior parameters are chosen as μ₀ = −2.10 × 10⁻³ and ${\sigma }_0^2=9.8\times {10}^{-7}$.

We set iteration interval length L_it = 10 for the RUL prediction with right prior information and L_it = 5 for that with inappropriate prior information. The corresponding PDFs of RULs under M₁ and M₂, and the actual RULs at some different CM points are shown in Figure 7. From Figure 7, it can be observed that, for the case with the right prior information, the range of the PDFs of the RULs based on both methods could cover the actual RUL. However, the uncertainty in the estimated RULs of M₁ is less than that by M₂. For the case with inappropriate prior information, the PDFs based on M₁ could cover the actual RULs. However, the PDFs based on M₂ could not cover the actual RULs and thus the maintenance is delayed if such predictive result is applied in maintenance schedule. This indicates that the proposed method could obtain desirable results regardless of the quality of the prior information.

To further test the goodness of fit, we use the mean square error (MSE) to compare the methods in terms of the total prediction error of the RUL. The MSE at each observation point is calculated as follows:

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

(34)

where T is the real failure time for a specific data path. Figure 8 presents the results of the MSEs at some CM points. This criterion demonstrates that our method is better than the traditional Bayesian method. For the case with inappropriate prior information, the accuracy of the proposed method improves significantly. The reason is that, after more iteration of the EM algorithm, the estimation is more adaptive to the characteristic of the individual item under study. Moreover, when the prior information is appropriate, M₁ is a little better than M₂. This implies that the Bayesian method is useful for the case with appropriate prior information.

Since the combination of the Bayesian method and the EM algorithm is actually the MLE method for a fixed degradation rate [40], we further compare our method with the MLE method. The corresponding PDFs of RULs and the actual RULs at some different CM points are shown in Figure 9. From Figure 9, we observe that all the PDFs could cover the actual RUL. However, the MLE method has a wider distribution, and has low probability for the time close to the actual RUL values.

Additionally, we give the results of the MSEs at some CM points as shown in Figure 10. It can be observed that the proposed method obtains better accuracy than the MLE method and has certain robustness over the selection of the prior information. Moreover, the MLE method cannot apply the prior information for the RUL prediction. This further demonstrates the superiority of the proposed method.

Overall, the numerical examples and the practical study demonstrate that our developed method can work well and efficiently. Moreover, we verify that reasonably integrating the prior information and in situ information into degradation modeling can improve the accuracy of the RUL prediction.

6. Conclusions

This paper proposes a novel RUL prediction algorithm for lithium-ion batteries based on the WPME. First, some issues regarding the RUL prediction for the WPME have been studied. In order to ensure that w − x_k > 0, we use the TND to model the estimated state w − x_k and obtain an exact and closed-form RUL distribution by considering the ME and the distribution of the estimated drift parameter simultaneously. A new likelihood function is established for the offline parameters estimation, and thus the requirement of a fixed drift parameter for a specific item is satisfied. For the online parameter updating, we infer that the combination of the Bayesian updating algorithm and the EM algorithm derives the same results with that by the MLE method with a fixed drift parameter, and verify the interesting consequence through numerical examples. Based on this interesting result, we propose a heuristic parameter updating algorithm. Finally, the usefulness of the proposed method is demonstrated by a real-world degradation data of lithium-ion batteries from NASA. Compared with the existing approach, the proposed method can generate better results in predicting the RUL and has application potential.

We primarily discuss the issues associated with estimating the RUL for the batteries with linear degradation characteristic. However, in many cases, lithium-ion batteries exhibit nonlinear degradation trend, especially when they approach the end of life. Therefore, using Wiener process for the RUL estimation of lithium-ion batteries with nonlinear degradation trend is necessary in future research. Additionally, the relaxation effect could lead to the recovery of the battery, i.e., increases the available capacity for the next cycle. This phenomenon gives rise to a new studying issue that how to set the rest time for increasing the utilization of the batteries. Moreover, the heuristic parameter updating algorithm proposed in this paper is still preliminary, the issue regarding how to reasonably incorporate the prior information and in situ information together should be further exploited.