Online Semiparametric Identification of Lithium-Ion Batteries Using the Wavelet-Based Partially Linear Battery Model

Abstract

Battery model identification is very important for reliable battery management as well as for battery system design process. The common problem in identifying battery models is how to determine the most appropriate mathematical model structure and parameterized coefficients based on the measured terminal voltage and current. This paper proposes a novel semiparametric approach using the wavelet-based partially linear battery model (PLBM) and a recursive penalized wavelet estimator for online battery model identification. Three main contributions are presented. First, the semiparametric PLBM is proposed to simulate the battery dynamics. Compared with conventional electrical models of a battery, the proposed PLBM is equipped with a semiparametric partially linear structure, which includes a parametric part (involving the linear equivalent circuit parameters) and a nonparametric part [involving the open-circuit voltage (OCV)]. Thus, even with little prior knowledge about the OCV, the PLBM can be identified using a semiparametric identification framework. Second, we model the nonparametric part of the PLBM using the truncated wavelet multiresolution analysis (MRA) expansion, which leads to a parsimonious model structure that is highly desirable for model identification; using this model, the PLBM could be represented in a linear-in-parameter manner. Finally, to exploit the sparsity of the wavelet MRA representation and allow for online implementation, a penalized wavelet estimator that uses a modified online cyclic coordinate descent algorithm is proposed to identify the PLBM in a recursive fashion. The simulation and experimental results demonstrate that the proposed PLBM with the corresponding identification algorithm can accurately simulate the dynamic behavior of a lithium-ion battery in the Federal Urban Driving Schedule tests.

1. Introduction

Currently, lithium-ion batteries play a significant role in the energy storage devices used in electric vehicles (EVs) because of their high energy density, low self-discharge rate, lack of a memory effect, high operating voltage, and long life cycle. The primary concerns in the design of EVs are how to maintain optimum battery performance and extend the battery’s expected life. To achieve these goals, the battery must have a well-designed battery management system (BMS). The BMS must track the dynamic behaviors of the battery for reliable and efficient operation, which requires battery models that can accurately determine battery behavior under various operating conditions. Another important task for the battery model is to provide precise battery data for simulation of an EV energy system.

Due to the complex chemical and physical processes that occur in batteries, developing an accurate battery model is challenging. Researchers worldwide have developed a wide variety of battery models for different purposes. Almost all the existing battery models can be classified into the following two types: (1) parametric models, e.g., electrochemical models [1,2] and electrical models [3,4,5,6]; and (2) nonparametric models, e.g., artificial neural network (ANN) models [7,8,9] and abstract mathematical models [10]. Parametric battery models have been developed in terms of the equivalent electric-circuit parameters (for electrical models) or electrochemical parameters (for electrochemical models). These parameters have clear physical meaning and can be used to investigate the working status of batteries. For example, electrical model parameters, such as the open-circuit voltage (OCV) and equivalent internal resistance, can provide useful information that can be used to estimate the state of charge (SoC) and state of health (SoH) based on the battery’s intrinsic characteristics [11]. Battery parameters typically rely on many factors, such as the SoC, SoH, internal structure, service time, and operating conditions. Some of these factors are difficult or almost impossible to measure, which may lead to misspecification of the models and cause severe modeling biases. In contrast, the nonparametric modeling approach is useful for exploring hidden structures and reducing modeling biases. For example, ANN-based nonparametric battery models were developed in [7] and [9]; in these reports, variables that are hard to measure were replaced by the weights of the neurons. However, the nonparametric battery models lack direct physical interpretation, which limits their fields of application.

To overcome the above-mentioned disadvantages, this paper presents a semiparametric battery model based on the Thevenin equivalent circuit model, which is a compromise between the nonparametric and parametric battery models. Since the introductory work of Engle in 1986 [12], the semiparametric model has been widely studied and applied in many fields, such as applied statistics [13], signal processing [14,22], and dynamic system modeling [23,24,25]. However, to our knowledge, this model has not yet been investigated for use in battery modeling. The proposed semiparametric battery model consists of two major components: a linear parametric component that involves the equivalent circuit parameters (such as ohmic internal resistance, polarization resistance, and polarization capacitance) and a nonparametric component that involves the OCV (because accurate modeling of OCV is a difficult task). In the literature, the semiparametric model described by such a form is also termed the “partially linear model” because of its relationship to the standard linear model [14]. Compared with the standard linear model, the partially linear model could be more adaptive when it is believed that the response variable has a linear relationship with some covariates but a nonlinear or unknown relationship with others. The proposed semiparametric battery model, hereafter referred to as the partially linear battery model (PLBM), features the flexibility of the nonparametric component while maintaining the direct physical interpretation of the parametric component.

Because the battery parameters cannot be measured directly during the run time, electrical models are often estimated through model-based optimization techniques, such as the use of recursive least square (RLS) estimators [5,15], sliding-mode observers [16], the adaptive control approach [17], the subspace identification method [18,19], and the extended Kalman filter [20,21]. For the PLBM, an efficient estimation procedure must be developed for both the linear parameters and nonparametric component. Several semiparametric estimators [14,22,23,24,25] based on basis function expansion methods have been proposed to identify the partially linear model. The current study is concerned with modeling the nonparametric component in the PLBM using the wavelet multiresolution analysis (MRA) expansion due to the wavelet’s good time—frequency localization property. Another advantage of the modeling method is that the nonparametric component can be parsimoniously represented by large timescale wavelet basis functions (or a truncated wavelet MRA expansion) because of the slowly varying behavior of the OCV. Furthermore, the wavelet MRA expansion typically leads to a sparse representation of the original signal [26], such as the OCV; this sparsity can be exploited to improve the estimation performance by employing an l₁-norm penalized wavelet estimator. However, most existing penalized wavelet estimators of the partially linear model are processed in batch form using iterative algorithms [14,23,27], which are not suitable for online implementation purposes. Recent advances in sparse linear model identification have shown that the cyclic coordinate descent (CCD) algorithm provides an efficient means of solving the penalized least squares (LS) problem [28,29] and can be implemented online in a recursive fashion [30,31]. This feature motivates us to extend this algorithm to the wavelet-based partially linear model and develop a recursive penalized wavelet estimator based on a modified online CCD algorithm. The performance of the proposed semiparametric identification approach for lithium-ion batteries is finally investigated by performing simulations and experiments.

2. Modeling of the Lithium-Ion Battery

The modeling of a battery is aimed at investigating the working status of the battery and finding a suitable mathematical expression that can accurately describe its input-to-output dynamics. Accurate electrical battery modeling is difficult to achieve because the electrochemical process is complex and a variety of operating conditions are used. Therefore, many complex electrical modeling methods have been developed to compensate for the modeling errors and uncertainties but at the costs of high computational burden, system complexity, and resource requirements. In this section, we present a simple, flexible, and comprehensive semiparametric battery model based on straightforward circuit analysis of the Thevenin equivalent circuit model.

2.1. Battery Equivalent Circuit Equations

In general, the electrical model is described by an equivalent circuit based on a combination of basic elements, such as voltage sources, resistors, and capacitors, to approximate the electrochemical processes and input-to-output dynamics of the battery. By combining the basic elements in different ways, the electrical models can be further divided into four classes [3]: the Rint model, radio control model, Thevenin equivalent circuit model, and Partnership for a New Generation of Vehicles model. Among these models, the Thevenin equivalent circuit model is widely employed to model the lithium-ion battery because it has a simple model structure and can capture the battery’s dynamic response well within an acceptable error [17,21]. Therefore, we select this model as the basic element of our proposed model.

The equivalent circuit model is shown in Figure 1 and consists of three parts: (1) the equivalent ohmic internal resistor R_s; (2) the resistor-capacitor (RC) parallel network C_p // R_p (where R_p is the equivalent polarization resistance and C_p is the equivalent polarization capacitance), which is used to simulate transient responses of the battery during charging–discharging transients; and (3) the OCV v_oc(h(t)), which is a nonlinear function of SoC h(t). The equivalent circuit model considers the current as the model control input and the terminal voltage as the measured output.

Figure 1. Equivalent circuit model of the lithium-ion battery.

Figure 1. Equivalent circuit model of the lithium-ion battery.

Theoretically, the equivalent circuit parameters R_s, R_p, and C_p in this model are functions of the temperature and SoC. In fact, because the variations in the SoC and temperature with respect to time are both very small, these parameters can be assumed as quasi-stationary, i.e., time invariant over a short observation time window. Further, for real-time application, parameter identification can only be performed over a finite time sliding window or exponentially decaying time window of the most recent measurement. Thus, within some error tolerance, these parameters can be assumed as constant during the identification. In other words, online identification can capture these parameters faster than the temperature or SoC variations.

The analysis of the equivalent circuit shown in Figure 1 is straightforward. Using Kirchhoff’s law, the dynamics of the circuit can be expressed as follows:

$$v_b\left(t\right)=v_{oc}\left(h\right(t\left)\right)-R_s{i}_b\left(t\right)-v_c\left(t\right)$$

(1)

$$\frac{\mathrm{d}{v}_c\left(t\right)}{\mathrm{d}t}=-\frac{1}{C_p{R}_p}{v}_c\left(t\right)+\frac{1}{C_p}{i}_b\left(t\right)$$

(2)

where v_b(t) and i_b(t) are the terminal voltage and current, respectively; and v_c(t) is the voltage across the RC network, which cannot be measured directly.

2.2. PLBM

The OCV of a battery is a nonlinear function of its SoC. Theoretically, if the nonlinear relationship OCV–SoC is obtained prior to the identification procedure, then the model parameters can be directly identified using the LS estimators. Traditionally, this relationship is experimentally established in the offline mode by measuring the OCV at each SoC set point and assuming that it remains the same during the lifetime of the battery. However, in practice, this relationship is not exactly the same for every battery, even if these batteries are part of the same production batch. Moreover, this relationship depends on the number of charge–discharge cycles of the batteries. Therefore, for any online identification algorithms, the a priori OCV–SoC data may lead to an inaccurate identification result. Thus, we should treat v_oc(h(t)) as an unknown function that belongs to some functional space, which must be identified online, as well as the other model parameters. In this case, v_oc(h(t)) can be simply represented as v_oc(t).

As can be seen from Equations (1) and (2), the equivalent circuit equations involve the OCV v_oc(t) and the polarization voltage v_c(t), both of which are unknown. However, the unknown term v_c(t) can be removed from these equations based on the following three assumptions [17]: (1) The variation in SoC with respect to time is very small because the consumed (or regained) energy is very small compared with the total useful capacity; (2) The battery temperature is monitored in real-time and controlled at a predetermined level by a well-designed BMS module. Thus, the temperature variation can be ignored for normal operating conditions; (3) The usage history of the battery has a long-term effect on battery behavior, which can be ignored during on-line identification procedure. These assumptions also have been applied in [21]. With the unknown term v_c(t) removed, the dynamics of the battery model can be approximated as follows:

$$\frac{\mathrm{d}{v}_b\left(t\right)}{\mathrm{d}t}=-\frac{v_b\left(t\right)}{C_p{R}_p}-R_s\frac{\text{d}{i}_b\left(t\right)}{\text{d}t}-\frac{R_p+R_s}{C_p{R}_p}{i}_b\left(t\right)+\frac{v_{oc}\left(t\right)}{C_p{R}_p}.$$

(3)

The Laplace transformation of (3) can be expressed as follows:

$$V_b\left(s\right)=-(R_s+\frac{R_p}{1+C_p{R}_{p}s})I_b\left(s\right)+\frac{v_{oc}\left(s\right)}{1+C_p{R}_{p}s}.$$

(4)

Equation (4) is a mathematical model of a continuous system in the s-domain, which can be mapped to its discrete form as follows:

$$v_b\left(k\right)={\theta }_1{v}_b(k-1)+{\theta }_2{i}_b\left(k\right)+{\theta }_3{i}_b(k-1)u\left(k\right)$$

(5)

where v_b(k) and i_b(k) are the samples of v_b(t) and i_b(t), respectively, at time point t_k = kT_c and T_c is the sampling period. θ₁, θ₂, and θ₃ are defined as:

$$\begin{gathered} {\theta }_1=C_p{R}_p/(C_p{R}_p+T_c) \\ {\theta }_2=-R_s-R_p{T}_c/(C_p{R}_p+T_c) \\ {\theta }_3=C_p{R}_p{R}_s/(C_p{R}_p+T_c) \end{gathered}$$

(6)

and u(k) is defined as:

$$u\left(k\right)=(1-{\theta }_1)v_{oc}\left(k\right)$$

(7)

by letting v_oc(k) denote the virtual sample of v_oc(t) at time point t_k = kT_c , which is assumed to be an unknown function because v_oc(k) is unknown.

We should note that v_b(k) and i_b(k) are acquired only via observation, and during the observation, these parameters often suffer additional random disturbance. Let y(k) and x(k) denote the observed values of v_b(k) and i_b(k), respectively. From (5), the discrete observation equation can be obtained as follows:

$$y\left(k\right)={\mathit{\phi }}_k^T\mathit{\theta }+u\left(k\right)+\epsilon \left(k\right)$$

(8)

where the superscript T denotes the transpose of a vector or matrix, ${\mathit{\phi }}_k={[y(k-1),x\left(k\right),x(k-1)]}^T$ is the measurement vector, $\mathit{\theta }={[{\theta }_1,{\theta }_2,{\theta }_3]}^T$ is the unknown linear parameter vector, and ε(k) is assumed to be an additive zero-mean white noise error that contains not only the measurement errors of v_b(k) and i_b(k) but also the interior noise of the battery system.

The battery model (8) is nonlinear because the OCV exhibits a nonlinear relationship with the charging–discharging current. In general, nonlinear system identification and modeling techniques can be classified into three categories: parametric, nonparametric, and, more recently, semiparametric methods. Parametric methods assume that the model structure based on the physical modeling principles is known but that the model parameters are unknown. Parametric methods can provide a direct physical interpretation of the system’s behavior. Conversely, nonparametric methods are necessary when the physical-based model structure is unknown. Nonparametric methods assume a functional relationship between the input and output data. In this case, these methods rely on functional approximators, such as the neural network, the kernel function, or abstract mathematical regression models. However, these methods focus on the input–output representation ability and neglect some important properties of the highly successful parametric modeling, such as the simple model structure involved and clear physical meaning. As an attractive alternative, semiparametric methods offer a compromise between nonparametric and parametric methods, thus reducing the high risk of model misspecification related to the fully parametric methods while simultaneously avoiding some of the disadvantages of fully nonparametric methods. From (8), we can see that the battery model is equipped with a semiparametric, partially linear structure, i.e., the right-hand side of (8) consists of two major components: a linear parametric component ${\mathit{\phi }}_{\mathrm{k}}^{\mathrm{T}}\mathit{\theta }$ involving the parameter vector θ and a nonparametric component u(k) involving the unknown function v_oc(k). Thus, we can naturally describe the dynamics of the battery using the semiparametric partially linear model; this description is useful not only because our understanding of the OCV does not provide us with a specific function or a parametric model of OCV but also because it provides a better understanding of the model structure and the parameters.

3. Online Identification of the Wavelet-Based PLBM

The goal of our study is to estimate both the unknown linear parameter vector θ and the nonparametric component u(k) from the input–output data. One attractive approach in addressing this problem is to expand the nonparametric component as a linear combination of the parameterized basis functions. The problem is then reduced to a linear estimation problem because the output is linear not only with the linear components but also with the basis functions. The identification method developed in this section is based on the wavelet MRA expansion of the nonparametric component u(k) of the PLBM.

3.1. Wavelet MRA

To make this paper self-contained and improve its readability, we first provide a brief overview of the basic concepts of wavelet MRA theory, which are essential for this paper (the readers are referred to [32] and [33] for a detailed treatment of this subject). The wavelet MRA expansion has proven to be a powerful tool for parsimoniously representing arbitrary signals because of its inherent property of localization in both time and frequency; wavelet representations have been successfully used in dynamic system identification and modeling [32,33,34,35].

The essential idea of the wavelet MRA is to decompose a finite energy signal into a subspace, which is spanned using the scaling and wavelet basis functions at different resolutions. For any signal $f\left(t\right)\in L^2\left(R\right)$, the MRA expansion with respect to the mother wavelet ψ(t) and scaling function ϕ(t) can be expressed as [34]:

$$f\left(t\right)=\sum _{m=-\infty }^{\infty }{c}_{j_0,m}{\varphi }_{j_0,m}\left(t_k\right)+\sum _{j=j_0}^{\infty }\sum _{m=-\infty }^{\infty }{d}_{j,m}{\psi }_{j,m}\left(t\right),j_0,j,m\in Z$$

(9)

where the scaling basis function ${\varphi }_{j_0,m}\left(t\right)=2^{j_0/2}\varphi (2^{j_0}t-m)$ is the scaled and translated version of ϕ(t); the wavelet basis function ${\psi }_{j,m}\left(t\right)=2^{j/2}\varphi (2^{j}t-m)$ is the scaled and translated version of ψ(t); c_j0,m and d_j,m are respectively termed the scaling and wavelet coefficients at scale levels j₀; and j; and j₀; is an arbitrary integer that represents the coarsest resolution level (defined by the user). The wavelet MRA representation has localization property in both the time domain (via translations) and the frequency domain (via scale). With this property, wavelets can capture global (low frequency) and local (high frequency) characteristics of any finite energy signal easily.

Theoretically, the resolution level (or timescale) can span from −∞ to ∞. In practice, the observed signals are sampled with finite length. Thus, the MRA expansion can be performed only over a finite number of resolution levels. Without loss of generality, we assume that f(t) is defined over the time window [0,T]. Let f(k) denote the sample of f(t) at time point t_k = kT_c, k = 1, 2,···, N where N = T / T_c is the length of the sampled signal and assumed to be a power of two for simplicity. Then, the MRA expansion of f(k) can be expressed as:

$$f\left(k\right)=\sum _{m=0}^{2^{j_0}-1}{c}_{j_0,m}{\varphi }_{j_0,m}\left(t_k\right)+\sum _{j=j_0}^{j_{\max }-1}\sum _{m=0}^{2^j-1}{d}_{j,m}{\psi }_{j,m}\left(t_k\right),k=1,···,N$$

(10)

where j₀ ≥ 0 is the coarsest resolution level (or the largest timescale) and j_max = log₂ (N) is the finest resolution level (or the smallest timescale) where the sampled signal is originally represented.

From the viewpoint of signal representation, the scaling coefficients carry the approximation information of the original signal at the coarsest scale j₀, and the wavelet coefficients at scale j carry the information that is different between two approximations to the original signal at scales j and j − 1. When these smaller timescale components at scale j ≥ J make little contribution to the original signal, the original signal f(k) can be approximated by a truncated MRA expansion from scale j₀ up to J as:

$$f^J\left(k\right)=\sum _{m=0}^{2^{j_0}-1}{c}_{j_0,m}{\varphi }_{j_0,m}\left(t_k\right)+\sum _{j=j_0}^J\sum _{m=0}^{2^j-1}{d}_{j,m}{\psi }_{j,m}\left(t_k\right),J<j_{\max },k=1,···,N$$

(11)

where all of the wavelet coefficients at scale J ≤ j < j_max are assumed to be zero. The higher the upper resolution level J, the more accurate the approximation. The truncated MRA representation can describe not only low-frequency fluctuations but also some rapid localized changes. The advantage of the truncated MRA expansion is that it provides a parsimonious approximation representation of the original signal and is particularly suitable for original signals that do not vary significantly during a short time interval.

3.2. Truncated Wavelet MRA Expansion of the Nonparametric Component

As mentioned earlier, the OCV of a battery has a complex nonlinear relationship with its SoC. Obtaining this relationship in real time is difficult because the dynamics of the battery are actually the result of a series of complex electrochemical processes affected by the operating conditions. The variation in SoC with respect to time is very small because the consumed (or regained) energy is very small compared with the total useful capacity. Thus, as a monotonic function with respect to SoC, the OCV is considered to be slowly varying within the battery’s usable capacity region [17]. Based on this consideration, we assume that the OCV is the superposition of physical and chemical phenomena that occur at larger timescales. In other words, the OCV can be modeled as a signal that belongs to a subspace spanned by large-scale wavelets. Moreover, we assume that (1 − θ₁), another part involved in the nonparametric component, is time invariant during the identification process. Thus, an appropriate model for the nonparametric component u(k) in (8) can be provided using a truncated MRA expansion as:

$$u^J\left(k\right)=\sum _{m=0}^{2^{j_0}-1}{c}_{j_0,m}^u{\varphi }_{j_0,m}\left(t_k\right)+\sum _{j=j_0}^J\sum _{m=0}^{2^j-1}{d}_{j,m}^u{\psi }_{j,m}\left(t_k\right),k=1,···,N$$

(12)

where the truncation parameter J is defined as the finest timescale that characterizes the complexity of the OCV. The truncated MRA expansion leads to a parsimonious model structure that is highly desirable for identification of the wavelet-based, partially linear model discussed in the next section.

The MRA expansion of a signal is not unique and depends on the type of wavelet used. In this study, we use the Haar wavelet because of its simple structure and good localization properties [32]. Furthermore, the Haar scaling and wavelet basis functions are completely supported and provided explicitly, and these features enable the operations of the MRA expansion of (12) to be more convenient.

The Haar scaling function ϕ^haar(t) and mother wavelet ψ^haar(t) can be described as:

$${\varphi }^{Haar}\left(t\right)=\{\begin{array}{cc}1,\hfill & 0<t\le 1\\ 0,\hfill & otherwise\end{array}$$

(13)

and:

$${\psi }^{Haar}\left(t\right)=\{\begin{array}{cc}1,\hfill & 0<t\le 1/2\\ -1,\hfill & 1/2<t\le 1\\ 0,\hfill & otherwise\end{array}$$

(14)

For notational convenience, we denote ϕ^haar(t) and ψ^haar(t) as g₁(t) and g₂(t); that is, g₁(t) and g₂(t) are the scaling and wavelet basis functions at the coarsest scale of zero, respectively. All other Haar wavelet basis functions derived from g₂(t) are denoted as:

$$g_l\left(t\right)=2^{j/2}{g}_2(2^{j}t-m),\hspace{1em}\hspace{1em}\hspace{1em}l=2^j+m+1$$

(15)

where m = 0,1,···,2^j − 1 for each scale j = 1,2,···,J. We define g(t) = [g₁(t), g₂(t),···,g_L(t)]^T, where L = 2^J+1. For example, Figure 2 shows the Haar basis functions of g(t) when J = 3. The discrete representation of g(t) is defined as:

$${\mathit{g}}_k={\left[g_1\right(\frac{k}{N}),g_2(\frac{k}{N}),···,g_L(\frac{k}{N}\left)\right]}^T,\hspace{1em}\hspace{1em}k=1,···,N$$

(16)

Figure 2. Haar basis functions.

Figure 2. Haar basis functions.

Using the above notations and from (12), the truncated MRA expansion of u(k) from the initial scale j₀ = 0 up to J can be compactly expressed as:

$$u^J\left(k\right)=\sum _{l=1}^L{\eta }_l{g}_l\left(\frac{k}{N}\right)={\mathit{g}}_k^T\mathit{\eta },\hspace{1em}\hspace{1em}k=1,···,N$$

(17)

where the expansion coefficient vector is defined as η = [η₁,η₂,···,η_L]^T with η_l as the expansion coefficient that corresponds to the associated Haar wavelet basis function g_l(t).

3.3. Recursive Penalized Wavelet Estimator for Online PLBM Identification

By substituting u(k) in (8) with the truncated MRA representation u^J(k), the observation equation can be rewritten as:

$$y\left(k\right)={\mathit{\phi }}_k^T\mathit{\theta }+{\mathit{g}}_k^T\mathit{\eta }+{\epsilon }_k={\mathit{z}}_k^T\mathit{\beta }+\epsilon \left(k\right),\hspace{1em}\hspace{1em}k=1,···,N$$

(18)

where ${\mathit{z}}_k^T=[{\mathit{\phi }}_k^T,{\mathit{g}}_k^T]$ and ${\mathit{\beta }}^T=[{\mathit{\theta }}^T,{\mathit{\eta }}^T]$ are both P = 3 + L -dimensional vectors. Comparing (18) with (8), it is seen that the semiparametric identification problem has been transformed into a linear-in-parameter identification problem, which can be solved by the ordinary LS (or RLS) estimator as follows:

$$\hat{\mathit{\beta }}=\underset{\beta }{\mathrm{argmin}}\sum _{k=1}^N{\left(y\right(k)-{\mathit{z}}_k^T\mathit{\beta })}^2$$

(19)

The wavelet MRA expansion typically leads to a sparse representation of the original signal [26], meaning that most of the wavelet coefficients are zero or close to zero, whereas the important features of the signal are captured by a small number of large coefficients. However, the ordinary LS and RLS estimators cannot account for the underlying sparsity. For the wavelet-based partially linear model (18), the sparsity can be exploited to improve the estimation performance by the penalized wavelet estimator, similar to that in [14] and [23]:

$$\hat{\mathit{\beta }}=\underset{\beta }{\mathrm{argmin}}{J}_N(\mathit{\beta })$$

(20)

where the cost function is:

$$J_N(\mathit{\beta })=\frac{1}{2}\sum _{k=1}^N{\left(y\right(k)-{\mathit{z}}_k^T\mathit{\beta })}^2+\lambda \sum _{i\ge 5}^P|{\beta }_i|$$

(21)

where λ is a given penalty factor; and β_i is the ith entry of β. In (21), the first summation term measures how well the candidate solution fits the observed data in the LS sense, whereas the second summation term is a regularizer that considers an l₁-norm penalty on the wavelet coefficients. λ controls the relative weight of the two terms, which means that the larger λ is, the more the wavelet coefficients will shrink to zero. The penalty term in (21) penalizes only the wavelet coefficients of the nonparametric part of the model and not the scaling coefficients. The penalized wavelet estimator (20) can be regarded as an extension of the wavelet shrinkage estimators, which are typically processed in batch form by iterative soft-thresholding algorithms [14,22,23]. However, the batch estimators suffer from high computational complexity and increased memory requirements as time progresses and are thus not appropriate for online implementation. In contrast to batch estimators, their recursive counterparts offer computational and memory savings and enable tracking of the slowly time-varying system. This feature motivated us to develop a recursive penalized wavelet estimator for PLBM identification.

Now, we consider the recursive solution of (20). Neglecting the constant terms, the cost function (21) at the kth recursive step can be equivalently expressed as follows:

$$J^{\left(k\right)}(\mathit{\beta })=\frac{1}{2}{\mathit{\beta }}^{\mathit{T}}{\mathbf{R}}^{\left(k\right)}\mathit{\beta }-{\mathit{\beta }}^{\mathit{T}}{\mathbf{r}}^{\left(k\right)}+\lambda \sum _{i\ge 5}^P|{\beta }_i|,\hspace{1em}\hspace{1em}k=1,···,N$$

(22)

where the matrix R^(k) and vector r^(k) can be defined recursively as follows:

$${\mathbf{R}}^{\left(k\right)}={\mathbf{R}}^{(k-1)}+{\mathbf{z}}_k{\mathbf{z}}_k^T$$

(23)

and:

$${\mathbf{r}}^{\left(k\right)}={\mathbf{r}}^{(k-1)}+y\left(k\right){\mathbf{z}}_k$$

(24)

Unfortunately, the cost function (22) is nondifferentiable and is thus not available in closed-form minimization. One attractive approach to solve such an optimization problem is to run an online CCD algorithm due to its speed and numerical stability [30,31]. The CCD algorithm separately minimizes the cost function (22) for each entry of β and can admit a closed-form solution.

Let ${\widehat{\beta }}^{(k)}={[{\widehat{\beta }}_1^{\left(k\right)},{\widehat{\beta }}_2^{\left(k\right)},···,{\widehat{\beta }}_p^{\left(k\right)}]}^T$ denote the solution of the CCD algorithm at the kth recursive step. The pth entry of ${\widehat{\beta }}^{(k)}$ is cyclically updated as follows:

$${\widehat{\beta }}_p^{\left(k\right)}=\underset{\beta }{\mathrm{argmin}}{J}^{\left(k\right)}({\widehat{\beta }}_1^{\left(k\right)},···,{\widehat{\beta }}_{p-1}^{\left(k\right)},\beta ,{\widehat{\beta }}_{p+1}^{(k-1)},···,{\widehat{\beta }}_p^{(k-1)})$$

(25)

In the above expression, only the pth entry of β is currently optimized. We recall that the penalty in (22) is applied only on the wavelet coefficients (β₅:β_P) and not on the scaling coefficient (β₄) or linear parameters (β₁:β₃). Thus, (25) can be rewritten as:

$${\widehat{\beta }}_p^{\left(k\right)}=\{\begin{array}{ll}\underset{\beta }{\mathrm{argmin}}[\frac{1}{2}{R}^{\left(k\right)}(p,p){\beta }^2-w_p^{\left(k\right)}\beta ],& 1\le p\le 4\\ \underset{\beta }{\mathrm{argmin}}[\frac{1}{2}{R}^{\left(k\right)}(p,p){\beta }^2-w_p^{\left(k\right)}\beta +\lambda |\beta |],& 5\le p\le P\end{array}$$

(26)

where R^(k)(p,p) is the (p,p)th entry of the matrix R^(k) , and $w_p^{\left(k\right)}$ is defined as:

$$w_p^{\left(k\right)}=r^{\left(k\right)}\left(p\right)-\sum _{q=1}^{p-1}{R}^{\left(k\right)}(p,q){\widehat{\beta }}_q^{\left(k\right)}-\sum _{q-p+1}^P{R}^{\left(k\right)}(p,q){\widehat{\beta }}_q^{(k-1)}$$

(27)

where r^(k)(p) is the pth entry of r^(k). Being a scalar minimization problem, the solution of (26) can be easily obtained in closed form as:

where sgn(a) = a/|a| and [a]₊ = max (a,0). Because the Haar wavelet basis functions are compactly supported, local data processing can be employed to reduce the computational burden [35], i.e., the wavelet coefficients ${\widehat{\beta }}_p^{\left(k\right)},p=5,6,\mathrm{...},P,$ are updated only when at least one of the incoming inputs falls into the support of the associated wavelet functions. The modified online CCD algorithm for the recursive penalized wavelet estimator is summarized as follows:

1: Initialize ${\widehat{\beta }}^{\left(0\right)}={\mathbf{0}}_P$, ${\mathbf{R}}^{\left(0\right)}={\mathbf{0}}_{P\times P}$ and ${\mathbf{r}}^{\left(0\right)}={\mathbf{0}}_P$.

2: For k = 1,2,...,N do

3: Update R^(k) and vector r^(k) using (23) and (24), respectively.

4: For p = 1,···,4 do

5: Calculate $w_p^{\left(k\right)}$ using (27).

6: Update ${\widehat{\beta }}_p^{\left(k\right)}$ using (28.a).

7: end for

8: For p = 5,···,P do

9: If g_{p − 3}(k/N) = 0

10:  Do not update ${\widehat{\beta }}_p^{\left(k\right)}$.

11: else

12:   Calculate $w_p^{\left(k\right)}$ using (27).

13:   Update ${\widehat{\beta }}_p^{\left(k\right)}$ according to the thresholding rule (28.b).

14:  end if

15: end for

16: end for

After acquiring ${\widehat{\mathit{\beta }}}^{\left(k\right)}$, i.e., ${\widehat{\mathit{\theta }}}^{\left(k\right)}$ and ${\widehat{\mathit{\eta }}}^{\left(k\right)}$ (the recursive estimations of θ and η), the nonparametric component u(k) can be approximately calculated via (17) as follows:

$$u\left(k\right)\approx u^J \left(k\right)={\mathit{g}}_k^T{\widehat{\mathit{\eta }}}^{\left(k\right)},\hspace{1em}\hspace{1em}\hspace{1em}k=1,···,N$$

(29)

Also, according to (6) and (7), the trajectory of the estimated battery equivalent circuit parameters and OCV can be determined using u(k) and ${\widehat{\mathit{\theta }}}^{\left(k\right)}$.

4. Simulation and Experimental Results

This section discusses the numerical simulations and experiments performed to verify the accuracy and effectiveness of the proposed identification approach in accordance with the Federal Urban Driving Schedule (FUDS) test. The FUDS test is a standard driving test used by vehicle designers to evaluate performance under urban driving conditions. In all simulations and experiments, the nonparametric components of the PLBM are expanded using the truncated Haar wavelet MRA expansion with the coarsest resolution level j₀ = 0 and finest resolution level J = 5. The original input-output data in the observation window [0,T] are normalized into the unit interval [0,1] for convenient implementation [33].

4.1. Simulations

We first verify the accuracy of the proposed identification method by performing simulation studies. For all simulations, the input {x(k)} is acquired from one FUDS testing cycle as model excitation, whereas the output {y(k)} is generated according to (5) with an initial value y(0) = 0. The real linear parameters θ₁, θ₂ and θ₃ are generated according to (6) with R_s = 50 mΩ, R_p = 30 mΩ, and C_p = 2500 F. To illustrate the robustness of the proposed identification method, y(k) is corrupted by an additive white Gaussian noise with zero-mean and variance σ². The signal-to-noise ratio (SNR) is defined as 1/σ². The OCV is assumed to be monotonically varying throughout the entire simulation process and is modeled as a smooth polynomial for simplicity:

v_oc (k) = −0.5t_k (1.5 − t_k) + 41.5,  t_k = k/N,  k = 1,2,···,N

(30)

The trajectory of the estimated OCV, R_s, R_p, and C_p are shown individually in Figure 3 in comparison with their real values. The time axis is rescaled according to the inverse time mapping [0,1] → [0,T]. The penalty threshold adopted here is the traditional universal threshold value $\lambda =\sigma \sqrt{2\log \left(N\right)}$ derived from [26] and used in [14] and [23] for the wavelet semiparametric estimation. The SNR is set at 23 dB. The estimated OCV follows the true OCV variations quite well, and the other estimated parameters all converge to their real values with different convergence speeds. This result demonstrates that the proposed estimator can accurately identify the PLBM and time-invariant linear parameters.

Figure 3. Real and estimated results of the simulation: (a) OCV; (b) R_s; (c) R_p; (d) C_p.

Figure 3. Real and estimated results of the simulation: (a) OCV; (b) R_s; (c) R_p; (d) C_p.

To illustrate the flexibility of the proposed semiparametric approach, a model bias is first introduced to the OCV model as:

$$\text{Δ}=d(0.5+1.5t_k^2)$$

(31)

where d is defined as the bias coefficient. Thus, the OCV can be remodeled as:

$${\tilde{v}}_{oc}\left(k\right)=0.5t_k(1.5-t_k)+41.5+\text{Δ}$$

(32)

Assuming the OCV model is known as (30), and from (5) and (6), the parametric model of a battery can be expressed by the following linear regression equation:

[v_b(k) − v_oc(k)] = θ₁[v_b(k − 1) − v_oc(k)] + θ₂i_b(k) + θ₃i_b(k − 1)

(33)

The above parametric model can be identified by standard RLS estimator, and the identification results were compared with the proposed semiparametric approach as shown in Figure 4. As expected, the proposed semiparametric approach can accurately estimate the equivalent circuit parameters since it does not rely on any prior knowledge of OCV model.

Figure 4. Estimated results of equivalent circuit parameters with the semiparametric and parametric approaches: (a) R_s; (b) R_p; (c) C_p.

Figure 4. Estimated results of equivalent circuit parameters with the semiparametric and parametric approaches: (a) R_s; (b) R_p; (c) C_p.

It can be seen from Figure 4a that the parametric approach can also estimate the ohmic internal resistance R_s accurately. This is because the fact that the ohmic internal resistance represents the terminal voltage variation caused by a current variation in the battery and the slowly varying OCV bias has very little effect on the terminal voltage variation. However, as can be seen from Figures 4b,c the parametric approach cannot estimate the equivalent polarization resistance and capacitance accurately, and the larger the model bias (i.e., the larger d), the larger estimation error will be. This result indicates that the polarization parameters are much more sensitive to OCV model bias in comparison with the ohmic internal resistance.

We next compare the performance of the proposed recursive penalized wavelet estimator, the RLS wavelet estimator according to (19), and the adaptive control approach [17] with respect to the mean-squared error (MSE) of the OCV, which is defined as:

(34)

where E{} denotes the expectation operator and the lower bound of the summation is set as N₀ = 100 to ensure a fair comparison. To evaluate the MSE, we conducted 500 replicate Monte Carlo simulations for each SNR. Moreover, to clearly illustrate the sensitivity of each approach to observed noise, the input-output data are not pre-filtered. For the adaptive control approach, the OCV was identified based on the discrete version of Equations (11) and (12) in [17]. The results are shown in Figure 5 and illustrate the sensitivity of the three identification approaches to the observed noise. As expected, the performance of all three approaches increases with the SNR. In the low-SNR regime, the wavelet estimators clearly outperform the adaptive control approach due to the fact that the wavelet estimators identify OCV only using the large scale expansion coefficients as shown in (17), thus the high frequency noise ( spanned by small scale expansion coefficients) can be effectively reduced from the identified OCV. However, the MSE performance of the adaptive control approach increases significantly with increasing SNR. In the high-SNR regime, the MSEs of the adaptive control approach and the proposed recursive penalized wavelet estimator are almost indistinguishable. Additionally, it can be seen that the proposed penalized wavelet estimator outperforms the RLS wavelet estimator because of the power of the l₁-norm penalty.

Figure 5. MSEs of OCV plotted against SNR.

Figure 5. MSEs of OCV plotted against SNR.

From the above simulation results, it can be seen that the recursive penalized wavelet estimator works well for identifying the equivalent circuit parameters and OCV. The PLBM, as a semiparametric model, combines the simple linear structure and direct physical interpretation ability of the parametric model with the flexibility of the nonparametric model. This is done by modeling the unknown term u(k) in (8) as a nonparametric component, which exhibits high fitting flexibility compared to parametric approaches [4,15], i.e., the nonparametric approach can be used to accurately identify the OCV without making any prior assumptions on the functional form of OCV (except it is a slowly varying signal). Moreover, unlike the fully nonparametric approaches [7,9], the proposed semiparametric approach assumes the availability of equivalent circuit model that has explicit physical meaning, i.e., it can directly extracts the equivalent circuit parameters and OCV from the identification results.

It should be pointed out that the wavelet-based PLBM (18) is a generic battery model, which is capable of provide more rich and useful description of battery dynamics in comparison with the baseline model (8). For example, if only the Haar scale basis function at scale zero is used to expand the nonparametric component u(k), the wavelet-based PLBM will be reduced to the baseline model. Since the adaptive identification approaches proposed in [17,21] are based on the baseline model, a better performance of these adaptive identification approaches could be anticipated by employing the wavelet-based PLBM, which is an interesting topic for future study.

4.2. Experiments

The objective of the experiments is to validate whether the identified results properly represent the real input-to-output behavior of the battery. The experimental studies were conducted on a battery module composed of nine LiMn₂O₄ lithium-ion cells connected in series Each healthy cell has a nominal voltage of 3.6 V and a nominal capacity of 15 A·h. The experimental setup comprised a Digatron EV battery testing system EVT300-500, a battery management module, a controller area network communication unit, and a host computer running the built-in BMS software BTS-600, as shown in Figure 6. The EVT300-500 system includes a programmable charger/electronic load (which can charge or discharge the battery module with a maximum voltage of 500 V and maximum current of 300 A) and a real-time data acquisition unit that can collect test data, including terminal voltage, outflow current, temperature, and accumulated ampere-hours and watt-hours. The battery management module can also measure the terminal voltage and temperature of each cell in the battery module. The BTS-600 software can deal with the collected data and generate control signals for the EVT300-500 system based on the designed program. To reduce the influence of temperature fluctuation, the experiments were performed under a predefined temperature of 28 °C using a temperature-controlled chamber.

To predict the voltage response of the PLBM, the linear parameters and nonparametric component must first be identified experimentally. After full charging (SoC = 100%) followed by a rest period, the battery module was discharged to 90% SoC at a constant current of 5 A. The terminal voltage and current collected for one FUDS driving cycle are shown in Figure 7a,b, respectively. The sampling period was a half second, and the penalty threshold was empirically chosen as λ = 1.5. The results are shown in Figure 8.

Figures 8a–d show that, given appropriate initial parameter selection, all estimated parameters converged after 400 s. The OCV reduced slowly from its initial value of 34.036 V to the final value of 33.951 V because it decreased as a monotonic function of the SoC. The convergence rate of R_s was much faster than those of the other parameters under a steady value of 25.6 mΩ, which implied that the ohmic internal resistance R_s was less sensitive to the charge–discharge behavior of the lithium-ion battery than the other parameters. The polarized resistance R_p and capacitance C_p converged to 14.3 mΩ and 3900 F, respectively, with some fluctuations.

Based on the previously obtained identification results, we calculated the predicted value of the voltage response of the PLBM using a one-step forward predictor. The voltage response predicted by the model was compared with the measured voltage response as shown in Figure 9a. The relative error analysis result is also shown in Figure 9b. It can be seen that the voltage predicted by the model accurately tracked the measured value despite the large initial value error, and the relative voltage error converged into ± 0.25% of the error band only.

Figure 6. Schematic diagram of the battery testing system.

Figure 6. Schematic diagram of the battery testing system.

Figure 7. Measured input–output data in a FUDS cycle: (a) Terminal voltage; (b) Current.

Figure 7. Measured input–output data in a FUDS cycle: (a) Terminal voltage; (b) Current.

Figure 8. Identification results for a FUDS cycle: (a) OCV; (b) R_s; (c) R_p; (d) C_p.

Figure 8. Identification results for a FUDS cycle: (a) OCV; (b) R_s; (c) R_p; (d) C_p.

Figure 9. Validation results of the voltage responses in a FUDS cycle: (a) Measured and predicted battery voltage responses; (b) Relative error rate.

Figure 9. Validation results of the voltage responses in a FUDS cycle: (a) Measured and predicted battery voltage responses; (b) Relative error rate.

5. Conclusions

To characterize the dynamic behavior of lithium-ion batteries for EVs, this paper has proposed a novel semiparametric identification approach based on the wavelet-based PLBM model. For simplicity, the PLBM developed in this study was based on the first-order equivalent circuit model. However, the analysis results can be easily extended to higher order models. The essential features and advantages of the proposed identification scheme may be summarized as follows:

(1) As a semiparametric model, the PLBM features the flexibility of the nonparametric component while maintaining the direct physical interpretation of the parametric component. Due to its flexible structure, the PLBM can be identified even with poor preliminary knowledge of the OCV.

(2) Because the wavelet method provides a powerful tool for nonstationary signal analysis, we advocate use of the wavelet MRA expansion to model the nonparametric component of the PLBM.

(3) The identification algorithm is very simple because it requires only elementary computations, and simple recursive routines are available for updating both R^(k) and r^(k), similar to the RLS algorithm.

(4) To exploit the sparsity in the MRA representation and allow for online implementation, a recursive penalized wavelet estimator is proposed, which uses a modified CCD algorithm with l₁-norm penalty on the wavelet coefficients to identify the wavelet-based PLBM.

(5) Because the Haar wavelet basis functions are compactly supported, only a small fraction of the wavelet coefficients must be updated for every new measurement data obtained, resulting in low-complexity online computation.