An Adaptive Prediction Model for the Remaining Life of an Li-Ion Battery Based on the Fusion of the Two-Phase Wiener Process and an Extreme Learning Machine

Abstract

Lithium-ion batteries (LiBs) are the most important part of electric vehicle (EV) systems. Because there are two different degradation rates during LiB degradation, there are many two-phase models for LiBs. However, most of these methods do not consider the randomness of the changing point in the two-phase model and cannot update the change time in real time. Therefore, this paper proposes a method based on the combination of the two-phase Wiener model and an extreme learning machine (ELM). The two-phase Wiener model is used to derive the mathematical expression of the remaining useful life (RUL), and the ELM is implemented to adaptively detect the changing point. Based on the Poisson distribution, the distribution of the changing time is derived as a gamma distribution. To evaluate the theoretical results and practicality of the proposed method, we perform both numerical and practical simulations. The results of the simulations show that due to the precise and adaptive detection of changing points, the proposed method produces a more accurate RUL prediction than existing methods. The error of our method for detecting the changing point is about 4% and the mean prediction error of RUL in the second phase is improved from 4.39 cycles to 1.61 cycles.

1. Introduction

Due to their ability to decrease environmental pollution, electric vehicles (EVs) have become a research hotspot [1]. An EV system consists of a differential, a mechanical transmission system, an electric motor, a power converter, a battery management system (BMS), and a battery pack [2]. In EV systems, batteries are the most important component [3,4]. Although there are many types of batteries that can be used in EVs, lithium-ion batteries (LiBs) are the most widely used due to their high energy density and long cycle life [2]. When using a battery, its state of health (SOH) will decrease over time. However, due to current engineering practice, the SOH of LiB typically decreases faster than intended due to random effects in real projects. When the SOH of the battery is too low for EVs to use, EVs may stop functioning, which may lead to property damage and casualties. Therefore, we must know the remaining useful life (RUL) of lithium-ion batteries in advance to avoid potential dangers and losses.

The methods used to obtain the RUL of a lithium-ion battery can be divided into three categories: Physical models; data-driven models; and hybrid models [5]. To build a physical model for LiB, we must know the internal electrochemical process of LiB. However, due to the complex structure of LiB, it is difficult to observe its internal electrochemical process [6]. Therefore, it is better to use a data-driven model to obtain the RUL of LiB. There are two types of data-driven models: Machine learning techniques and stochastic techniques.

In the literature on machine learning techniques, the number of artificial neural networks appears to be the most popular due to their strong ability to fit, classify, and predict LiB RULs. For example, Hossain et al. [7] proposed a method that uses a time-delay neural network (TDNN) to estimate the state of charge (SOC) of LiB. Furthermore, Wu et al. [8] provided a method that combines a feed forward neural network (FFNN) with importance sampling to predict the RUL of LiB. Additionally, Ren et al. [9] and Chemali et al. [10] used the deep learning approach to extract features of LiB and then predicted the RUL of LiBs based on these features. Because the weights of the neural network in references [7,8,9,10] are adjusted via continuous iteration, if many data are used, the time complexity will be sufficiently high, leading to a decrease in the accuracy of the RUL prediction. To decrease the time complexity, the extreme learning machine (ELM) can be used to predict the RUL of the LiB. ELM is a single hidden layer FFNN whose weights are obtained by random assignment and calculation; thus, we can determine its weights faster than the neural network described above [11]. Razavi-Far et al. [12] used ELM to predict the RUL of LiB using one-step and multistep prediction. Moreover, Hossain et al. [13] presented an improved ELM with the gravitational search algorithm (GSA) to estimate the SOC of LiB and experimentally showed that their method performed better than the back-propagation neural network (BPNN) and the radial basis function neural network (RBFNN). Additionally, Ma et al. [14] proposed a multilayer ELM whose initial weights were obtained by the heuristic Kalman algorithm (HKA) to predict the RUL of the LiB. However, with the addition of a hidden layer, their method will be more difficult to train than the neural network proposed in [14]. Therefore, Ma et al. [15] presented a model called the broad learning extreme learning machine (BL-ELM) to improve the ELM precision.

Although machine learning techniques are widely used, they tend to accumulate errors, and longer prediction times tend to produce more accumulated errors. Another drawback of machine learning techniques is that their RUL predictions are always expected RULs, which cannot describe the uncertainty within the degradation trajectory.

With their strong ability to describe the uncertainty within the degradation trajectory, stochastic techniques have been widely investigated and applied to many degradation systems [16]. There are three common stochastic processes (the gamma process (GP), inverse Gaussian process (IGP), and Wiener process (WP)), which are typically used in stochastic techniques. Peng et al. [17] and Chen et al. [18] presented a method based on IGP to describe the degradation of a device. Furthermore, Wu et al. [19] used the GP to predict the RUL of LiB. However, the IGP and GP can only depict a monotonous degradation trajectory [20]. Therefore, the IGP and GP will not be useful if any non-monotonous portion exists. However, regardless of whether it is monotonous, the WP performance will not change, and thus, it is better to use the WP to predict the RUL of LiB. Li et al. [21] and Tang et al. [22] presented a linear Wiener process model (LWPM) to predict the RUL of LiB and provided a method of updating parameters online to describe the differences between individuals. However, based on the survey reported by Burgess et al. [23], the capacity degradation of batteries can be divided into two phases: A slow degradation process, followed by a faster degradation process. Therefore, it is not suitable to describe the degradation of LiB using LWPM. Zhang et al. [20], Wang et al. [24,25], and Kong et al. [26] used the two-phase Wiener process model (TPWPM) to describe a device with two different degradation rates and determined the RUL of the TPWPM. The TPWPM can mitigate the drawbacks of the LWPM; however, the changing point detection of the TPWPM is not as good in these methods. Zhang [20] used the maximum likelihood estimation (MLE) to estimate the changing time without any prior knowledge; thus, the changing time detection of Zhang’s method had a low accuracy. Additionally, the distribution of the changing time in TPWPM was not explained clearly in [20]. Wang [24,25] and Kong [26] determined the changing time by Bayesian theory and the expectation-maximization algorithm, respectively; however, they did not consider the randomness of the changing time. Therefore, the methods of Wang and Kong could not adaptively eliminate the difference between the different units.

In this paper, we primarily focus on the changing point detection of TPWPMs. We combine the ELM with the TPWPM and propose a method called ELM-TPWPM to monitor the changing point adaptively in the degradation trajectory of LiB. First, the distribution of the changing time is inferred by the Poisson distribution. Then, we obtain the health index (HI) of LiB based on the one-step prediction by ELM. Therefore, we can determine the degradation phase using the HI. Then, we can update the RUL prediction of the TPWPM adaptively based on the identified degradation phase.

The remainder of this paper is organized as follows. Section 2 describes the problem. The mathematical expression of RUL for TPWPM is derived in Section 3. The parameter estimation of the TPWPM is proposed in Section 4. The experiment is presented in Section 5. Section 6 provides the conclusion to the study.

2. Problem Description

Based on the analysis of Burgess [23], there are two different degradation rates in the degradation of batteries: A slow degradation process, followed by a faster degradation process. For example, Figure 1 shows the degradation data of the LiB used in [27,28,29].

The ordinate of Figure 1 indicates the lost capacity of LiB and the abscissa is the time of the experiment. The LiB data in Figure 1 were obtained from the University of Maryland, and the four LiBs were subjected to the same experimental conditions. As shown in Figure 1, the degradation rate in the period of 0 to 70 cycles is much slower than the degradation rate after 70 cycles, which shows the specific characteristic of two-phase degradation. Moreover, many researchers have already started the TPWPM of LiB [20,24,25,26], thus, it is better to describe the degradation of LiB by a TPWPM than an LWPM. The TPWPM consists of two LWPMs with different degradation rates. As Wang et al. [30] discussed, the mathematical expression of TPWPM is defined as follows:

$$X\left(t\right)=\{\begin{array}{c}{x}_0+{\mu }_{1}t+{\sigma }_{1}B\left(t\right) 0<t<\tau \\ x_{\tau }+{\mu }_{2}t+{\sigma }_{2}B\left(t\right) \tau <t\end{array}$$

(1)

where ${\mu }_1\mathrm{and}{\mu }_2$ are drift coefficients in the first and second phases, respectively; ${\sigma }_1\mathrm{and}{\sigma }_2$ are diffusion coefficients in the first phase and second phase, respectively; $B\left(t\right)$ is the standard Brownian motion; $\tau$ is the time of the changing point; $x_0$ is the initial state in the first phase; and $x_{\tau }$ is the state at time $\tau$. Then, with the concept of the first passage time (FPT), we can obtain the mathematical expression of the LiB lifetime as follows:

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

(2)

where $T$ is the lifetime of LiB and $\omega$ is the failure threshold of LiB. To describe the uncertainty in the degradation of LiB, let $f_T\left(t\right)$ represent the probability density function (PDF) of $T$.

Based on the discussion in [22], the lifetime of the LWPM obeys an inverse Gaussian distribution. The TPWPM is composed of two LWPMs in two different periods of time; thus, the lifetime of the TPWPM also obeys an inverse Gaussian distribution in two different periods of time. If all parameters of the TPWPM are known in advance, the PDF of the TPWPM is as follows:

$$f_T\left(t\right)=\{\begin{array}{c} \frac{\omega -x_0}{\sqrt{2\pi t^3{\sigma }_1^2}}exp\left(-\frac{{\left(\omega -x_0-{\mu }_{1}t\right)}^2}{2{\sigma }_1^{2}t}\right) 0<t<\tau \\ \frac{\omega -x_{\tau }}{\sqrt{2\pi {\left(t-\tau \right)}^3{\sigma }_2^2}}exp\left(-\frac{{\left(\omega -x_{\tau }-{\mu }_2\left(t-\tau \right)\right)}^2}{2{\sigma }_2^2\left(t-\tau \right)}\right) \tau <t\end{array}$$

(3)

Equation (3) shows that the changing time $\tau$ influences the lifetime in both the first and second degradation phase. Therefore, if we want an accurate lifetime prediction, it is important to have an accurate estimation of $\tau$. The degradation data shown in Figure 1 demonstrate that different LiBs have different changing times. Therefore, before the arrival of $\tau$, we should regard it as a random variable, and when $\tau$ arrives, we should detect it adaptively so that the proposed TPWPM can be more widely applicable in the field of LiB.

From the description above, the following issues will be researched in this paper:

(1)

How to obtain a precise estimation of changing time $\tau$.

(2)

How to derive a suitable distribution for $\tau$.

3. RUL of TPWPM

As we discussed above, we must know the value of all parameters in Equation (3) to determine the lifetime of LiB using Equation (3). It is impossible to know the true value of $\tau$ in advance; thus, we cannot know the true value of $x_{\tau }$ in advance. Therefore, to obtain $f_T\left(t\right)$, we must obtain the PDF of $x_{\tau }$ first. Let $g_{\tau }\left(x_{\tau }\right)$ denote the PDF of $x_{\tau }$. To obtain the $g_{\tau }\left(x_{\tau }\right)$, Lemma 1 is as follows:

Lemma 1.

If$X\left(t\right)={\mu }_{ }t+{\sigma }_{ }B\left(t\right)$denotes a linear Brownian motion with the initial value$x_0=0$, its transition density function with an absorbing boundary$\xi$has the following form [31]:

$$g\left(x,t\right)=\frac{1}{\sqrt{2\pi t{\sigma }_B^2}}\left\{exp\left(-\frac{{\left(x-\mu t\right)}^2}{2{\sigma }_B^{2}t}\right)-exp\left(\frac{2\mu \xi }{{\sigma }_B^2}\right)exp\left(-\frac{{\left(x-\mu t-2\xi \right)}^2}{2{\sigma }_B^{2}t}\right)\right\}$$

(4)

Based on (4), we can obtain the mathematical expression of $g_{\tau }\left(x_{\tau }\right)$ as follows:

$$\begin{array}{cc}\hfill g_{\tau }\left(x_{\tau }\right)=\frac{1}{\sqrt{2\pi \tau {\sigma }_1^2}}& \{exp\left(-\frac{{\left(x_{\tau }-{\mu }_1\tau \right)}^2}{2{\sigma }_1^2\tau }\right)\hfill \\ & \hfill -exp\left(\frac{2{\mu }_1\omega }{{\sigma }_1^2}\right)exp\left(-\frac{{\left(x_{\tau }-{\mu }_1\tau -2\omega \right)}^2}{2{\sigma }_1^2\tau }\right)\}\end{array}$$

(5)

Then, we must determine the other parameters in Equation (3). Based on Figure 1, although the four LiBs are the same type of battery and are tested under the same experimental conditions, their degradation trajectories are different, which indicates that their ${\mu }_1$, ${\mu }_2$, and $\tau$ vary. Therefore, to describe the uncertainty of the parameters, we typically think of ${\mu }_1$, ${\mu }_2$, and $\tau$ as random variables. Based on [20,21,22], ${\mu }_1\mathrm{and}{\mu }_2$ are assumed to be normally distributed; thus, we assume that ${\mu }_1~N\left({\mu }_{1p},{\sigma }_{1p}^2\right) \mathrm{and} {\mu }_2~N\left({\mu }_{2p},{\sigma }_{2p}^2\right)$, where ${\mu }_{1p} \mathrm{and} {\sigma }_{1p}^2$ are the mean value and variance of ${\mu }_1$, respectively, and ${\mu }_{2p} \mathrm{and} {\sigma }_{2p}^2$ are the mean value and variance of ${\mu }_2$, respectively. Then, we use the Poisson distribution to infer the PDF of $\tau$ and provide Theorem 1 as follows:

Theorem 1.

When there aren measurements of degradation before the changing point, and n obeys a Poisson distribution, the time of the changing point obeys a gamma distribution.

Due to space limitations, the derivation of Theorem 1 is shown in Appendix A.

Typically, the RUL of a LiB at time $t_k$ has more engineering significance than the lifetime shown in Equation (2). Based on the concept of Equation (2), the RUL of a LiB at time $t_k$ has the following form:

$$L_k=inf\left\{l_k:X\left(t_k+l_k\right)\ge \omega |X\left(t_k\right)\le \omega \right\}$$

(6)

where $L_k$ denotes the RUL of an LiB at time $t_k$. Similar to $f_T\left(t\right)$, we let $f_L\left(l_k\right)$ represent the PDF of $L_k$.

Based on Equations (3), (5) and (6), and considering the distribution of ${\mu }_1$ and ${\mu }_2,$ the conditional PDF of the RUL of the TPWPM at time $t_k$ is as follows [20]:

When $0<t_k<\tau$,

$$f_L\left(l_k\mid \tau \right)=\{\begin{array}{cc}\frac{\omega -x_k}{\sqrt{2\pi \left({\sigma }_{1p}^2{l}_k^2+{\sigma }_1^2{l}_k\right)l_k^2}}exp\left(-\frac{{\left(\omega -x_k-{\mu }_{1p}{l}_k\right)}^2}{2\left({\sigma }_{1p}^2{l}_k^2+{\sigma }_1^2{l}_k\right)}\right)0<l_k+t_k\le \tau & \\ A-B\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\tau <l_k+t_k\end{array}$$

(7)

where

$$\begin{gathered} A=\frac{1}{\sqrt{2\pi {\left(l_k^{}-\tau +t_k\right)}^2\left({\sigma }_a^2+{\sigma }_b^2\right)}}exp\left(-\frac{{\left({\mu }_a-{\mu }_b\right)}^2}{2\left({\sigma }_a^2+{\sigma }_b^2\right)}\right) \\ \left\{\frac{{\mu }_b{\sigma }_a^2+{\mu }_a{\sigma }_b^2}{{\sigma }_a^2+{\sigma }_b^2}\mathsf{\Phi }\left(\frac{{\mu }_b{\sigma }_a^2+{\mu }_a{\sigma }_b^2}{\sqrt{\left({\sigma }_a^2+{\sigma }_b^2\right){\sigma }_a^2{\sigma }_b^2}}\right)+\sqrt{\frac{{\sigma }_a^2{\sigma }_b^2}{\left({\sigma }_a^2+{\sigma }_b^2\right)}}\mathsf{\Phi }\left(\frac{{\mu }_b{\sigma }_a^2+{\mu }_a{\sigma }_b^2}{\sqrt{\left({\sigma }_a^2+{\sigma }_b^2\right){\sigma }_a^2{\sigma }_b^2}}\right)\right\} \\ B=exp\left(\frac{2{\mu }_{1p}\left(\omega -x_k\right)}{{\sigma }_1^2}+\frac{2\left({\left(\omega -x_k\right)}^2{\sigma }_{1p}^4\tau +{\left(\omega -x_k\right)}^2{\sigma }_{1p}^2{\sigma }_1^2\right)}{\left({\sigma }_1^2+\left(\tau -t_k\right){\sigma }_{1p}^2\right){\sigma }_1^4}\right) \\ \frac{1}{\sqrt{2\pi {\left(l_k^{}-\tau +t_k\right)}^2\left({\sigma }_a^2+{\sigma }_b^2\right)}}exp\left(-\frac{{\left({\mu }_a-{\mu }_c\right)}^2}{2\left({\sigma }_a^2+{\sigma }_b^2\right)}\right) \\ \left\{\frac{{\mu }_c{\sigma }_a^2+{\mu }_a{\sigma }_b^2}{{\sigma }_a^2+{\sigma }_b^2}\mathsf{\Phi }\left(\frac{{\mu }_c{\sigma }_a^2+{\mu }_a{\sigma }_b^2}{\sqrt{\left({\sigma }_a^2+{\sigma }_b^2\right){\sigma }_a^2{\sigma }_b^2}}\right)+\sqrt{\frac{{\sigma }_a^2{\sigma }_b^2}{\left({\sigma }_a^2+{\sigma }_b^2\right)}}\mathsf{\Phi }\left(\frac{{\mu }_c{\sigma }_a^2+{\mu }_a{\sigma }_b^2}{\sqrt{\left({\sigma }_a^2+{\sigma }_b^2\right){\sigma }_a^2{\sigma }_b^2}}\right)\right\} \\ {\mu }_a={\mu }_{2p}\left(l_k-\tau +t_k\right),{\mu }_b=\omega -x_k-{\mu }_{1p}\left(\tau -t_k\right) \\ {\mu }_c=-\omega +x_k-{\mu }_{1p}\left(\tau -t_k\right)-\frac{{\sigma }_{1p}^2\left(\tau -t_k\right)}{{\sigma }_1^2} \\ {\sigma }_a^2={\sigma }_2^2\left(l_k-\tau +t_k\right)+{\sigma }_{2p}^2{\left(l_k-\tau +t_k\right)}^2,{\sigma }_b^2={\sigma }_1^2\left(\tau -t_k\right)+{\sigma }_{1p}^2{\left(\tau -t_k\right)}^2 \end{gathered}$$

when $\tau <t_k$,

$$f_L(l_k|\tau )=\frac{\omega -x_k}{\sqrt{2\pi \left({\sigma }_{2p}^2{l}_k^2+{\sigma }_2^2{l}_k\right)l_k^2}}exp\left(-\frac{{\left(\omega -x_k-{\mu }_{2p}{l}_k\right)}^2}{2\left({\sigma }_{2p}^2{l}_k^2+{\sigma }_2^2{l}_k\right)}\right)$$

(8)

Then, considering the distribution of $\tau$, we can determine the PDF of the RUL in TPWPM by Bayesian theory. The PDF is calculated by the following formula:

$$f_L\left(l_k\right)={\int }_{t_k}^{\infty }{f}_L\left(l_k|\tau \right)p\left(\tau \right)d\tau$$

(9)

where $p\left(\tau \right)$ is the PDF of $\tau$, which obeys the gamma distribution. The integral in Equation (9) is difficult to calculate; however, we can obtain the numerical solution of $f_L\left(l_k\right)$ using certain mathematical methods.

4. Parameter Estimation of TPWPM

In this section, we estimate the parameter of TPWPM by combining offline estimation with online estimation. Specifically, the ELM is used during online estimation to monitor the arrival of changing points in real time, which can decrease the error caused by the offline estimation.

4.1. Offline Method of Parameter Estimation

Let the parameter of the $i$th LiB be ${\theta }_i=\left({\mu }_{i,1},{\mu }_{i,2},{\sigma }_{i,1}^2,{\sigma }_{i,2}^2,{\tau }_i\right)$, where ${\mu }_{i,1}\mathrm{and}{\mu }_{i,2}$ are the drift coefficients of the first phase and second phase, respectively; ${\sigma }_{i,1}^2\mathrm{and}{\sigma }_{i,2}^2$ are the diffusion coefficients of the first phase and second phase, respectively; and ${\tau }_i$ is the changing time of the $i$th LiB. We assume that we have acquired $n$ sets of LiB degradation data, which are obtained under the same conditions. The degradation data of the $i$th LiB are $X_i=\left(x_{i,1},x_{i,2}\dots \dots x_{i,k_i}\right)$, where $k_i$ is the total number of degradation data of the $i$th LiB. Therefore, the time of the degradation data can be described as $T_i=\left(t_{i,1},t_{i,2}\dots \dots t_{i,k_i}\right)$. To analyse and dispose of the data more conveniently, the following hypotheses are introduced:

Hypothesis 1 (H1).

The time in $T_i$has the same interval.

Hypothesis 2 (H2).

${\tau }_i$belongs to$T_i$.

Then, we use the maximum likelihood estimation (MLE) to estimate ${\theta }_i$. Because ${\tau }_i$ is the changing time, $X_{i,1:{\tau }_i}=\left\{x_{i,1},x_{i,2}\dots \dots x_{i,{\tau }_i}\right\}$ and $X_{i,{\tau }_i+1:k_i}=\left\{x_{i,{\tau }_i+1},x_{i,{\tau }_i+2}\dots \dots x_{i,k_i}\right\}$ are the data in the first and second phase, respectively. Their measurement times are $T_{i,1:{\tau }_i}=\left(t_{i,1},t_{i,2}\dots \dots t_{i,{\tau }_i}\right)$ and $T_{i,{\tau }_i+1:k_i}=\left(t_{i,{\tau }_i+1},t_{i,{\tau }_i+2}\dots \dots t_{i,k_i}\right)$, respectively. Based on Equation (1), $X\left(t\right)$ is normally distributed. Therefore, the likelihood function of the $i$th LiB under the condition of $X_{i,1:k_i}$ is as follows:

$$\begin{array}{c}\mathit{ln}L\left({\theta }_i|X_{i,1:k_i}\right)=\hfill \\ -\frac{k}{2}\mathit{ln}\left(2\pi \right)-\frac{1}{2}\mathit{ln}\left|{\phi }_{i,1}\right|-\frac{1}{2}\mathit{ln}\left|{\phi }_{i,2}\right|-\hfill \\ \frac{1}{2}{\left(X_{i,1:{\tau }_i}-x_{i,0}-{\mu }_{i,1}{T}_{i,1:{\tau }_i}\right)}^T{\phi }_{i,1}^{-1}\left(X_{i,1:{\tau }_i}-x_{i,0}-{\mu }_{i,1}{T}_{i,1:{\tau }_i}\right)\hfill \\ \hspace{1em}\hspace{1em}-\frac{1}{2}{\left(X_{i,{\tau }_i+1:k_i}-x_{i,{\tau }_i}-{\mu }_{i,2}{T}_{i,{\tau }_i+1:k_i}\right)}^T{\phi }_{i,2}^{-1}\left(X_{i,{\tau }_i+1:k_i}-x_{i,{\tau }_i}-{\mu }_{i,2}{T}_{i,{\tau }_i+1:k_i}\right)\hfill \end{array}$$

(10)

where ${\phi }_{i,1}$ and ${\phi }_{i,2}$ are the covariance matrices of $X_{i,1:{\tau }_i}$ and $X_{i,{\tau }_i+1:k_i}$, respectively, and

$${\phi }_{i,2}={\sigma }_{i,2}^2{\left[\begin{array}{cccc}{t}_{i,{\tau }_i+1}& t_{i,{\tau }_i+1}& \dots & t_{i,{\tau }_i+1}\\ t_{i,{\tau }_i+1}& t_{i,{\tau }_i+2}& \dots & t_{i,{\tau }_i+2}\\ ⋮& ⋮& \ddots & ⋮\\ t_{i,{\tau }_i+1}& t_{i,{\tau }_i+2}& \dots & t_{i,k_i}\end{array}\right]}_{\left(k_i-{\tau }_i\right)\ast \left(k_i-{\tau }_i\right)},{\phi }_{i,1}={\sigma }_{i,1}^2{\left[\begin{array}{cccc}{t}_{i,1}& t_{i,1}& \dots & t_{i,1}\\ t_{i,1}& t_{i,2}& \dots & t_{i,2}\\ ⋮& ⋮& \ddots & ⋮\\ t_{i,1}& t_{i,2}& \dots & t_{i,{\tau }_i}\end{array}\right]}_{{\tau }_i\ast {\tau }_i}$$

where $x_{i,0}$ is the initial value of the first phase of the $i$th LiB, and $x_{i,{\tau }_i}$ is the degradation state of the $i$th LiB at time ${\tau }_i$. When $\mathit{ln}L\left({\theta }_i|X_{i,1:k_i}\right)$ in Equation (10) is the largest, ${\theta }_i$ is the parameter we want. Therefore, we can obtain the largest ${\theta }_i$ by the following formula:

$${\theta }_i=\arg max\left(\ln L\left({\mu }_{i,1},{\mu }_{i,2},{\sigma }_{i,1}^2,{\sigma }_{i,2}^2,{\tau }_i|X_{1:k_i}\right)\right)$$

(11)

We can obtain the maximization of Equation (11) by the ‘fmincon’ function in MATLAB. For n LiBs, there are $\Theta =\{{\theta }_1,{\theta }_2\dots \dots {\theta }_n\}.$ Therefore, we can obtain the distribution parameters of ${\mu }_1,{\mu }_2,$ and $\tau$ by $\Theta .$.

Remark 1.

For the$i$th LiB, its degradation data from the initial state to failure have been collected so that${\mu }_{i,1}, {\mu }_{i,2}, and {\tau }_i$are constant rather than random variables.

4.2. Online Method of Parameter Estimation

In this part, we propose a method to monitor the working condition of LiB in real time so that we can decrease the estimation error by updating ${\mu }_1,{\mu }_2$, and $\tau$ in time.

4.2.1. Monitoring the Changing Point by ELM

The structure of ELM is shown in Figure 2. $\left\{x_1,x_2\dots \dots x_j\right\}$ and $\left\{h_1,h_2\dots \dots h_k\right\}$ are the input layer and hidden layer of the ELM, respectively, and the weights between them are random values. $\left\{{\beta }_1,{\beta }_2\dots \dots {\beta }_l\right\}$ represents the weights between the hidden layer and output layer and can be obtained mathematically.

We use the ELM to build the health index (HI) of a LiB to determine the changing point. First, we use the LiB degradation data in the first phase to train the ELM so that the degradation feature of LiB in the first phase is learned by the ELM. Then, the trained ELM is used to perform one-step prediction for the LiB that is working. The HI is constructed with the following formula:

$$H{I}_{1:m}=\left(h{i}_1,h{i}_2\dots \dots h{i}_m\right)=|Y_{1:m}-X_{1:m}|,$$

(12)

where $Y_{1:m}=\left(y_1,y_2\dots \dots y_m\right)$ is the output of the trained ELM and $X_{1:m}=\left(x_1,x_2\dots \dots x_m\right)$ is the real degradation data of the LiB. Because the trained ELM has learned the degradation feature of LiB in the first phase, if the input is the data from the second phase, the output will deviate significantly from normal, and the HI will be larger. Therefore, we construct the evaluation index of HI to adaptively detect the changing point. The form of the evaluation index is as follows:

$$\begin{array}{c}mea{n}_{HI}{}_{1:m}=\frac{{{\displaystyle \sum }}_{j=1}^{m}h{i}_j}{m}\\ sigm{a}_{HI}{}_{1:m}=\sqrt{\frac{{{\displaystyle \sum }}_{j=1}^m{\left(h{i}_j-mea{n}_{HI}{}_{1:m}\right)}^2}{m-1}}\end{array}$$

(13)

In the degradation model of Equation (1), $X\left(t\right)$ is normally distributed; thus, we assume that HI is also normally distributed. Therefore, we use the ‘$3\sigma$ criterion’ in a normal distribution to detect the changing point. When we obtain the new data $x_{m+1}$, its HI is $h{i}_{m+1}$. If $h{i}_{m+1}\notin [mea{n}_{HI}{}_{1:m}-3sigm{a}_{HI}{}_{1:m},mea{n}_{HI}{}_{1:m}+3sigm{a}_{HI}{}_{1:m}]$, $x_{m+1}$ will be classified as suspicious data. Let the number of inputs of ELM be n; when all of the data from $x_{m+1}$. to $x_{m+\mathrm{n}}$. are suspicious data, we judge that LiB has entered the second degradation phase starting from data $x_{m+1}$. Therefore, we mark the changing time as $\tau =m+1$. Because the one-step prediction of ELM is precise, we can always detect the changing point accurately.

4.2.2. Updating the Drift Coefficient by Bayesian Theory

We know from Section 2 that the ${\mu }_1\mathrm{and}{\mu }_2$ in TPWPM are random variables with distributions of ${\mu }_1~N\left({\mu }_{1p},{\sigma }_{1p}^2\right)\mathrm{and}{\mu }_2~N\left({\mu }_{2p},{\sigma }_{2p}^2\right), \mathrm{respectively}$. However, for a specific LiB, ${\mu }_1\mathrm{and}{\mu }_2$ are a certain number. Therefore, for a working LiB, we use Bayesian theory to update its ${\mu }_1\mathrm{and}{\mu }_2$ adaptively to determine the estimation of ${\mu }_1\mathrm{and}{\mu }_2$ as accurately as possible.

When the current time $t_m$ is below $\tau$, we update the distribution parameters of ${\mu }_1$. When $t_m>\tau$, we update the distribution parameters of ${\mu }_2$. The updating formula is as follows:

When $t_m<\tau$:

$$\begin{array}{c}{\mu }_{1p,m}=\frac{X_{1:m}^T{\phi }_{1,1:m}^{-1}{T}_{1:m}{\sigma }_{1p,0}^2+{\mu }_{1p,0}}{T_{1:m}^T{\phi }_{1,1:m}^{-1}{T}_{1:m}{\sigma }_{1p,0}^2+1}\\ {\sigma }_{1p,m}^2=\frac{{\sigma }_{1p,0_{}}^2}{T_{1:m}^T{\phi }_{1,1:m}^{-1}{T}_{1:m}{\sigma }_{1p,0}^2+1}\end{array}$$

(14)

where ${\mu }_{1p,m}\mathrm{and}{\sigma }_{1p,m}^2$ are the posterior estimations of ${\mu }_{1p}\mathrm{and}{\sigma }_{1p}^2$, respectively; ${\mu }_{1p,0}\mathrm{and}{\sigma }_{1p,0}^2$ are obtained from offline estimation and are the a priori estimations of ${\mu }_{1p}\mathrm{and}{\sigma }_{1p}^2$, respectively; $T_{1:m}=\left(t_1,t_2\dots \dots t_m\right)$ is the time of degradation data $X_{1:m}$; and ${\phi }_{1,1:m}={\sigma }_1^2{\left[\begin{array}{cccc}{t}_1& t_1& \dots & t_1\\ t_1& t_2& \dots & t_2\\ ⋮& ⋮& \ddots & ⋮\\ t_1& t_2& \dots & t_m\end{array}\right]}_{m\ast m}$ is the covariance matrix of $X_{1:m}$.

When $t_m>\tau :$

$$\begin{array}{c}{\mu }_{2p,m}=\frac{X_{\tau :m}^T{\phi }_{2,\tau :m}^{-1}{T}_{\tau :m}{\sigma }_{2p,0}^2+{\mu }_{2p,0}}{T_{\tau :m}^T{\phi }_{2,\tau :m}^{-1}{T}_{\tau :m}{\sigma }_{2p,0}^2+1}\\ {\sigma }_{2p,m}^2=\frac{{\sigma }_{2p,0_{}}^2}{T_{\tau :m}^T{\phi }_{2,\tau :m}^{-1}{T}_{\tau :m}{\sigma }_{2p,0}^2+1}\end{array}$$

(15)

where ${\mu }_{2p,m}\mathrm{and}{\sigma }_{2p,m}^2$ are the posterior estimations of ${\mu }_{2p}\mathrm{and}{\sigma }_{2p}^2,$ respectively; ${\mu }_{2p,0}\mathrm{and}{\sigma }_{2p,0}^2$ are obtained from offline estimation and are the a priori estimations of ${\mu }_{2p}\mathrm{and}{\sigma }_{2p}^2$, respectively; $T_{\tau :m}=\left(\tau ,\tau +1\dots t_m\right)$ is the time of degradation data $X_{\tau :m}$; and ${\phi }_{2,\tau :m}={\sigma }_2^2\left[\begin{array}{cccc}\tau & \tau & \dots & \tau \\ \tau & \tau +1& \dots & \tau +1\\ ⋮& ⋮& \ddots & ⋮\\ \tau & \tau +1& \dots & t_m\end{array}\right]$ is the covariance matrix of $X_{\tau :m}$. The proofs of Equations (14) and (15) are shown in Appendix B. Based on Equations (14) and (15) and the method proposed in Section 4.2.1, we can describe accurate and adaptive RUL prediction of a working LiB.

According to the method shown above, we provide the procedure of our algorithm in

Table 1. Implementation procedures of remaining useful life (RUL) estimation for LiB.

Algorithm Procedure
Step 1.	Use the offline estimation method in Section 4.1 to obtain the priori estimation of parameters.
Step 2.	Train the ELM using the LiB degradation data from the first phase.
Step 3.	Online estimation begins; obtain new degradation data of LiB.
Step 4.	Construct the health index (HI) using the trained ELM from Step 2 and the new data from Step 3.
Step 5.	Use the evaluation index in Equation (13) and the HI from Step 4 to judge whether a changing point is coming or not.
Step 6.	If a changing point appears, update the parameters in the second phase with Equation (15); otherwise, update the parameters in the first phase with Equation (14).
Step 7.	If a changing point appears, calculate the RUL of LiB with Equation (8); otherwise, calculate the RUL of LiB with Equation (7).
Step 8.	If the current LiB degradation data reaches the failure threshold, the algorithm ends; otherwise, loop the algorithm from Step 3.

.

5. Case Study

To verify the accuracy of the proposed algorithm, we provide two examples: (1) A numerical example with data created by a computer and (2) a practical example with data from [27,28,29].

5.1. Numerical Example

Because we cannot know the real parameters of realistic LiB degradation data, we used a computer to stimulate the degradation data of the TPWPM using set parameters. Then, we used these data to verify the online parameter estimation ability of the proposed method. These parameters were set as follows: ${\mu }_{1p}=0.005; {\mu }_{2p}=0.02; {\sigma }_{1p}^{ }=9\times {10}^{-4};{\sigma }_{2p}^{ }=4.5\times {10}^{-3}; {\sigma }_1^2=1\times {10}^{-4};{\sigma }_2^2=3\times {10}^{-4}; \alpha =140; \mathrm{and}\beta =1$, where $\alpha \mathrm{and} \beta$ are the parameters of changing time $\tau$, and $\tau ~Ga\left(\alpha ,1/\beta \right)$. We obtained 100 sets of data, and there were 300 data points in each set. The data we obtained are shown in Figure 3.

First, we used the ELM, which is shown in Section 4.2.1, to detect the changing point. The input dimensions and hidden nodes are important for ELM: Too many input dimensions will cause long time delays for changing point detection, and too many hidden nodes may cause overfitting. Therefore, to determine the number of input dimensions and hidden nodes in the ELM, we performed experiments with different numbers of input dimensions and hidden nodes, and used the relative error of changing time estimation (RECTE) to evaluate them. The formula of RECTE is as follows:

$$RECT{E}_{1:q}=\frac{{\sum }_{i=1}^q\left|\frac{{\tau }_{rea{l}_i}-{\tau }_{detectio{n}_i}}{{\tau }_{rea{l}_i}}\right|}{q}$$

(16)

where ${\tau }_{real}$ and ${\tau }_{detection}$ are the real changing time and the changing time detected by the ELM, respectively, and $q$ is the number of changing times detected by the ELM. For $\tau ~Ga\left(140,1\right)$, $p\left(\tau \le 100\right)\approx 0.01\%$, which indicates that the probability of the degradation data from time 1 to 100 belonging to the first phase is $99.99\%$. Therefore, we used the data from time 1 to 100 to train the ELM, and the number of trained sets was 70. To reduce the randomness of the experiment and prove the robustness of our method, we performed 100 experiments for the same input dimensions and hidden nodes. Then, we took the mean RECTE of 100 trials as the final RECTE. The final RECTE is shown in Figure 4.

Figure 4 shows that when we choose five hidden nodes and five input dimensions, the RECTE is the smallest. However, the RECTE of three input dimensions and four hidden nodes is near the RECTE of five input dimensions and five hidden nodes; thus, the fewer input dimensions there are, the shorter the delay of changing point detection. Therefore, it is better to choose three input dimensions and four hidden nodes for the ELM structure.

Then, to prove the robustness of changing point detection, we created 1000 sets of data using the parameters shown above and used the ELM with three input dimensions and four hidden nodes to detect the changing point. We created another 100 sets of data to train the ELM, and 1000 sets of data were used as test data. We used the detection deviation between the detection time and real time of the changing point to evaluate the accuracy of the ELM. The probability of each detection deviation is shown in Figure 5. Then, we compared the proposed method with Zhang’s method [20], and the results shown in Figure 5 indicate that the proposed method for detecting the changing point is more efficient and accurate. The mean relative error of our method for detecting the changing point is about 4%. Compared with Figure 5b, the result of our method shows a strong robustness.

To update other parameters in real time, we created a set of degradation data, as shown in Figure 6. The real parameters of the data shown in Figure 6 are defined as follows: ${\mu }_{1,sample}=0.0059$; ${\mu }_{2,sample}=0.0202$; ${\sigma }_{1,sample}^2=1\times {10}^{-4}$; ${\sigma }_{2,sample}^2=3\times {10}^{-4}$; and ${\tau }_{sample}=128$. Then, we used (14) and (15) to update the parameters, and the initial values were calculated using 20 samples created in $N\left({\mu }_{1p},{\sigma }_{1p}^2\right)$ and $N\left({\mu }_{2p},{\sigma }_{2p}^2\right)$: ${\mu }_{sample,1p,0}=0.0053; {\mu }_{sample,2p,0}=0.0206; {\sigma }_{sample,1p,0}^{ }=1.1\times {10}^{-3}; \mathrm{and} {\sigma }_{sample,2p,0}^{ }=5.5\times {10}^{-3}$. Then, we used the data in Figure 6 to estimate the parameters in real time. The process of parameter updating is shown in Figure 7.

Figure 7 shows that as increasingly more degradation data are considered during parameter updating, ${\mu }_{1p} \mathrm{and} {\mu }_{2p}$ move closer to ${\mu }_{1,sample} \mathrm{and} {\mu }_{2,sample},$ respectively. ${\sigma }_{1p}^{ } \mathrm{and} {\sigma }_{2p}^{ }$ also decrease, which indicates that ${\mu }_{1p} \mathrm{and} {\mu }_{2p}$ become increasingly certain. The changing time we obtained is 126, and its deviation with ${\tau }_{sample}$ is 2. Therefore, the proposed method of parameter online estimation is efficient and precise.

5.2. Practical Example

We used the model shown in Equation (1) to predict the RUL of real LiBs and used the method described in Section 4 to estimate the parameter of these real LiBs. Based on the algorithm procedure shown in Table 1, we used realistic LiB degradation data from the Center for Advanced Life Cycle Engineering (CALCE) of Maryland University, College Park, MD, USA [27,28,29] (CS2-35, CS2-36, CS2-37, and CS2-38) to perform a practical experiment. The cell chemistry of these LiBs is as follows: LiCoO₂ cathode, and EDS results also showed trace elements of manganese. Each of these LiBs has a capacity of 1350 mAh and they underwent the same charging profile, which was a standard constant current/constant voltage protocol with a constant current rate of 0.5 C until the voltage reached 4.2 V and then 4.2 V was sustained until the charging current dropped to below 0.05 A. The discharge cut off voltage for these batteries was 2.7 V. Each LiB was cycled multiple times under the above conditions, and each cycle was about 30 h. All these cells were randomly numbered and named accordingly. The realistic LiB degradation data are shown in Figure 1. We eliminated the trend items of these data, and the results are shown in Figure 8.

Figure 8 shows that the random items of the LiB are not monotonous; thus, we cannot use the inverse Gaussian process [17,18] or the gamma process [19] to describe the degradation of LiB. Then, based on the offline parameter estimation method shown in Section 4.1, we used the data in CS2-35, CS2-36, CS2-37, and CS2-38 to estimate the offline parameters of CS2-35, CS2-36, CS2-37, and CS2-38, respectively, which are shown in

Table 2. Offline parameter estimates of CS2-35, CS2-36, CS2-37, and CS2-38.

	${\mathit{\mu }}_1$	${\mathit{\mu }}_2$	${\mathit{\sigma }}_1^2$	${\mathit{\sigma }}_2^2$	$\mathit{\tau }$
CS2-35	$3.75\times {10}^{-3}$	$2.28\times {10}^{-2}$	$1.61\times {10}^{-4}$	$3.86\times {10}^{-4}$	62
CS2-36	$4.97\times {10}^{-3}$	$2.32\times {10}^{-2}$	$1.67\times {10}^{-4}$	$2.65\times {10}^{-4}$	70
CS2-37	$4.16\times {10}^{-3}$	$2.31\times {10}^{-2}$	$1.31\times {10}^{-4}$	$3.77\times {10}^{-4}$	76
CS2-38	$3.95\times {10}^{-3}$	$2.16\times {10}^{-2}$.	$1.46\times {10}^{-4}$	$2.03\times {10}^{-4}$	78

.

To demonstrate the efficiency of the proposed method for $\tau$ detection, we regarded the parameters of CS2-36, CS2-37, and CS2-38 as prior estimations, and the parameters of CS2-35 were estimated by the proposed online method. Then, we used the RECTE shown in Equation (16) to evaluate the performance of the ELM with different structures. We used the history degradation data of CS2-36, CS2-37, and CS2-38 to obtain the RECTE of ELM with different structures. The degradation data of CS2-36 in the first phase were used to train the ELM, and the degradation data of CS2-37 and CS2-38 were used to construct the HI and detect the changing point. The RECTE of different ELM are shown in Figure 9.

Figure 9 shows that the RECTE of ELM with four input dimensions and six hidden nodes is the smallest; therefore, we chose this type of ELM to detect the changing point. Then, we began the real-time estimation using the degradation data of CS2-35. The degradation data of CS2-36, CS2-37, and CS2-38 in the first phase were used to train the ELM with four input dimensions and six hidden nodes. With the trained ELM, we used Equation (12) to construct the HI of CS2-35, and the HI at different times is shown in Figure 10.

With the obtained HI, we could determine the position of $\tau$ by the evaluation index shown in Equation (13), and we found that $\tau =62 \mathrm{cycles}$, which is the same as the offline estimation of $\tau$. Therefore, the estimation of $\tau$ shows that the proposed method is accurate. Based on $\tau$, we could update the distribution parameters of the drift coefficients ${\mu }_1~N\left({\mu }_{1,p},{\sigma }_{1,p}^2\right)$ and ${\mu }_2~N\left({\mu }_{2,p},{\sigma }_{2,p}^2\right)$ with Equations (14) and (15), and these results are shown in Figure 11.

${\mu }_{1,p}$ and ${\mu }_{2,p}$ in Figure 11 change with the degradation data, which shows that the proposed method of parameter updating has a fast response to the arrival of new degradation data. Conversely, when more degradation data are considered in the process of parameter updating, ${\sigma }_{1,p}^2$ and ${\sigma }_{2,p}^2$ decrease, which shows that the more degradation data we obtain, the more accurate the estimation of ${\mu }_1$ and ${\mu }_2$ will be. Then, based on the results of Figure 10 and Figure 11, we used Equations (7)–(9) to predict the RUL of CS2-35. We assumed that when the lost capacity of CS2-35 arrives at 0.75 Ah, the RUL of CS2-35 becomes 0. Therefore, according to the degradation data shown in Figure 1, the real lifetime of CS2-35 is 85 cycles. To demonstrate the accuracy of the proposed method, the methods of Razavi et al. [12], Li et al. [21], and Zhang et al. [20] were compared, and the results of the RUL prediction are shown in Figure 12.

Figure 12a shows that the proposed method for detecting the changing point can provide the real-time performance. As soon as the ELM detects the arrival of a changing point, the proposed model enters the second degradation phase so that the RUL prediction becomes much more accurate. The changing time of CS2-35 is much earlier than the others shown in Table 2, and when we predict the RUL of CS2-35 before the arrival of the changing time, we can only estimate the changing time of CS2-35 based on the changing time of the others. Therefore, it is inevitable that there is a relatively large error of the RUL in the first phase.

Figure 12b shows the RUL at different times obtained by different methods. To demonstrate the improved accuracy of RUL prediction due to the proposed changing point detection method, we used the mean error of RUL prediction (MERULP) to evaluate the accuracy of these methods in the second degradation phase. The calculation of MERULP is as follows:

$$MERUL{P}_{62:q}=\frac{{\sum }_{i=62}^q\left|RU{L}_{rea{l}_i}-RU{L}_{predic{t}_i}\right|}{q},$$

(17)

where $RU{L}_{rea{l}_i}$ is the $i\mathrm{th}$ real RUL in Figure 12b and $RU{L}_{predic{t}_i}$ is the $i\mathrm{th}$ RUL prediction value of different methods. There is a large error when we use Li’s method to obtain the RUL prediction, which shows that the linear Wiener model cannot describe the degradation path of LiB. Although the RUL prediction obtained by Razavi’s method appears better than the result of Li’s method, there is a drawback that the RUL prediction of Razavi’s method cannot be obtained early in the degradation path. This phenomenon can be attributed to the degradation speed being too slow for the ELM early in the degradation path. Therefore, the ELM is likely to create a negative degradation speed so that the degradation data can never arrive at the failure threshold. Compared to the two-phase Wiener model of Zhang, the proposed method yields a better performance in the second degradation phase due to the proposed changing point detection method.

To describe the error in the RUL prediction visually, we show the difference between the RUL prediction and the real RUL at different times in Figure 12c, which demonstrates that the RUL prediction of the proposed method is markedly better than those of the other methods in the second degradation phase. Based on Table 2, the changing time of CS2-35 is 62 cycles; thus, we started the calculation of $MERULP$ at 62 cycles. The results of $MERULP$ are shown in

Table 3. The mean error of RUL prediction ($MERULP$) of different methods.

	Zhang’s Method	Li’s Method	Razavi’s Method	Our Method
$MERULP$	4.39 cycles	20.7 cycles	5.17 cycles	1.61 cycles

. Table 3 demonstrates that, due to the proposed changing point detection method, the accuracy of the RUL prediction in the second phase improves considerably.

Overall, based on the results of the numerical and practical experiments, both the theoretical and practical results of the proposed method have been explained explicitly. The data created by the computer and the practical data were used to prove the validity of our method, which means the robustness of our method is good. The high accuracy of the proposed method for RUL prediction and changing point detection has been demonstrated, highlighting the ability of the proposed method to be applied in practical engineering.

6. Conclusions

This paper primarily considers the changing point detection and RUL prediction of the two-phase Wiener process model (TPWPM) for LiBs. We used ELM to construct the HI of a LiB and detect the changing point of the LiB by the HI. Then, based on the concept of FPT, we obtained the PDF of the RUL of the LiB and determined its formula. Finally, we investigated numerical and practical cases to determine the theoretical and practical results of the proposed method. Based on a comparison with other methods, the higher accuracy of the proposed method has been demonstrated. The error of our method for detecting the changing point is about 4% and the mean prediction error of RUL in the second phase is improved from 4.39 cycles to 1.61 cycles. However, there are several directions for future research:

(1)

The accuracy of RUL prediction in the first degradation phase must be improved.

(2)

The accuracy of changing point detection could be improved.