Multi-Level Model Reduction and Data-Driven Identification of the Lithium-Ion Battery

Abstract

The lithium-ion battery is a complicated non-linear system with multi electrochemical processes including mass and charge conservations as well as electrochemical kinetics. The calculation process of the electrochemical model depends on an in-depth understanding of the physicochemical characteristics and parameters, which can be costly and time-consuming. We investigated the electrochemical modeling, reduction, and identification methods of the lithium-ion battery from the electrode-level to the system-level. A reduced 9th order linear model was proposed using electrode-level physicochemical modeling and the cell-level mathematical reduction method. The data-driven predictor-based subspace identification algorithm was presented for the estimation of lithium-ion battery model in the system-level. The effectiveness of the proposed modeling and identification methods was validated in an experimental study based on LiFePO₄ cells. The accuracy and dynamic characteristics of the identified model were found to be much more likely related to the operating State of Charge (SOC) range. Experimental results showed that the proposed methods perform well with high precision and good robustness in the SOC range of 90% to 10%, and the tracking error increases significantly within higher (100–90%) or lower (10–0%) SOC ranges. Moreover, to achieve an optimal balance between high-precision and low complexity, statistical analysis revealed that the 6th, 3rd, and 5th order battery model is the optimal choice in the SOC range of 90% to 100%, 90% to 10%, and 10% to 0%, respectively.

1. Introduction

The lithium-ion battery is currently one of the most promising energy storage devices due to its high energy density, high power density, and long cycle life performance, and has been widely used in many fields, such as electric vehicles, energy storage, and industrial electronics [1]. In order to ensure the safe, reliable, and efficient operation of lithium-ion batteries, researchers need to understand the electrochemical and thermodynamic characteristics of the battery, which is usually costly and time-consuming through experimental and model-based approaches [2]. In addition, the performance of the battery can be strongly affected by operation conditions. These facts increase the difficulties in building the battery model and reducing computational complexity in applications with high real-time requirements. From an application-oriented perspective, the battery model needs to be simple and easy to calculate, while still being complex enough to provide valid and accurate results. This creates a problem where the models are often either too simplified and approximated to provide high-precision analysis or too complicated for low-complexity processing. However, little research to date has focused on the model reduction and optimization with the goal of achieving a balance between high-precision and low-complexity for real-time applications.

Multi-type models of the lithium-ion battery have been developed to study different issues, such as thermal behavior, inner polarization, microstructure, and capacity fading [2]. Moreover, many model-based techniques have been extensively applied to application-oriented researches, in particular with regard to evaluating the overall performance of lithium-ion batteries. At present, battery models can be divided into three categories: the equivalent-circuit models (ECM), the empirical models, and the electrochemical models.

The ECMs are the most used battery models and they mainly focus on capturing the input-to-output behavior of batteries using different lumped resistor–capacitor (RC) circuits. The parameters of the ECMs are usually obtained through system-identification approaches based on circuit theory. Many kinds of ECM topologies have been proposed, including the series RC circuit model [3,4], the parallel RC circuit model [5], the runtime-based circuit model [6,7], and the impendence-based circuit model [8,9]. The structures of ECMs are usually quite simple, and they have been widely applied in real-time applications. However, the parameters of ECMs have little connection with the electrochemical parameters of the battery. As a result, little physical insight into the physicochemical processes of battery can be inferred from the ECMs.

The empirical models belong to the category of data-based modeling approaches. They build the battery models based entirely on the experimental data instead of physicochemical theory. Many kinds of empirical models have been proposed using different techniques, such as statistical analysis [10,11], support vector [12], big data [13,14], and deep neural networks [15,16]. The empirical battery models mainly deal with the battery properties that related to capacity [17], State-of-Charge [18] and State-of-Health [19,20]. When using empirical-based approaches, researchers do not need to have a deep understanding of the electrochemical processes, which makes the empirical models easily implemented in practice. However, a common disadvantage of the empirical models is that experimental data are always needed; they are applicable only for the modeled experiment under a particular operating condition. Therefore, empirical models are not able to predict beyond a particular operating range.

The electrochemical models are also known as the first-principle models. The electrochemical model was originally developed by Doyle et al. based on the theories of porous electrodes and concentrated solutions [21]. Electrochemical models are usually coupled with many physicochemical processes including mass and charge conservations as well as the electrochemical kinetics. In addition, electrochemical models usually consist of a series of partial differential equations, (PDEs), which make them computationally expensive. Recently, many reduced-order electrochemical models [22,23,24,25] have been proposed for different purpose, including State-of-Charge estimation [25], aging prediction [26,27], optimal charging control [28,29], and the optimization of battery design [30]. Compared to the empirical models and ECMs, the electrochemical models can achieve higher accuracy with clear physical meanings. However, the calculation of the electrochemical model is a time-consuming process, and there are a large number of electrochemical parameters that need to be predefined, which lead to the electrochemical models being rarely used in real-time applications.

The goal of this research was two-fold. First, we aimed at proposing a reduced order model of the lithium-ion battery using electrode-level electrochemical modeling and cell-level mathematical reduction methods. The second goal was to present an easily implemented data-driven predictor-based subspace identification algorithm for the estimation of the battery model in the system-level. In particular, aimed to show that the accuracy and dynamic characteristics of the identified LiFePO₄ battery model were related to the operating State of Charge (SOC) range; statistical analysis revealed that 6th, 3rd, and 5th order battery model is the optimal choice in the SOC range of 90% to 100%, 90% to 10%, and 10% to 0%, respectively. Models presented in this paper can be used for the model-based estimation of the lithium-ion battery in real-time applications.

This paper is organized as follows. Section 2 introduces the electrode-level electrochemical model of the lithium-ion battery. Section 3 surveys the cell-level model reduction methods. Section 4 outlines the methods and experiments for the system-level model identification. In Section 5 the model validation tests as well as the statistical analysis of the test results are presented. Section 6 and Section 7 contain the discussion and conclusions, respectively.

2. Electrode-Level Modeling of Lithium-Ion Battery

In this section, the electrochemical fundamentals of lithium-ion batteries are introduced, and the electrode-level electrochemical modeling process is then studied. Symbols employed are listed below in Nomenclature. It should be noted that this research is based on the ideal assumption that we do not consider the side reactions, and we mainly focus on the main reactions within the lithium-ion battery, i.e., the lithium insertion/extraction reactions.

2.1. Governing Equations

Figure 1 shows a schematic diagram of the lithium-ion battery with three main domains: a negative electrode (width δ_n), a separator (width δ_sep), and a positive electrode (width δ_p). We can treat the lithium-ion battery as a one-dimensional model on electrode-level from the negative electrode (x = 0) to the positive electrode (x = L). During discharge, lithium ions diffuse to the surface of the active material particles in the negative electrode and migrate to the positive electrode through the electrolyte via diffusion and ionic conduction. Meanwhile, electrons are released through the external circuit to supply power to the load. The exact opposite process happens during charge.

To start, we summarize the four PDEs governing the charge and discharge dynamics of the lithium-ion battery in Equations (1) to (8). The governing equations were first developed by Doyle [21] and have been widely applied in many research studies [22,23,24,25]. They mainly describe the charge and species conservation both in the solid and electrolyte phases.

Conservation of species in solid phase:

$$\frac{\partial c_s}{\partial t}=\frac{D_s}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial c_s}{\partial r}\right),r\in (0,R_s),$$

(1)

with the boundary condition:

$${\frac{\partial c_s}{\partial r}|}_{r=0}=0,\begin{array}{c}{D_s\frac{\partial c_s}{\partial r}|}_{r=R_s}=-\frac{j^{\mathrm{Li}}}{a_{s}F}=-\frac{j^{\mathrm{Li}}{R}_s}{3{\epsilon }_{s}F}\end{array}.$$

(2)

Conservation of species in electrolyte phase:

$${\epsilon }_e\frac{\partial c_e}{\partial t}=D_e^{eff}\frac{{\partial }^2{c}_e}{\partial x^2}+\frac{(1-t_{+}^0)}{F}{j}^{\mathrm{Li}},x\in (0,L),$$

(3)

with the boundary condition:

$${\frac{\partial c_e}{\partial x}|}_{x=0}={\frac{\partial c_e}{\partial x}|}_{x=L}=0.$$

(4)

Conservation of charge in solid phase:

$${\delta }^{\mathrm{eff}}\frac{{\partial }^2{\varphi }_{\mathrm{s}}}{\partial x^2}-j^{\mathrm{Li}}=0,x\in (0,L).$$

(5)

with the boundary condition:

$$\{\begin{array}{l}{\frac{\partial {\varphi }_s}{\partial x}|}_{x={\delta }_n}={\frac{\partial {\varphi }_s}{\partial x}|}_{x={\delta }_n+{\delta }_{sep}}=0,\\ -{\delta }_n^{eff}{\frac{\partial {\varphi }_s}{\partial x}|}_{x=0}={\delta }_p^{eff}{\frac{\partial {\varphi }_s}{\partial x}|}_{x=L}=\frac{I}{A}.\end{array}$$

(6)

Conservation of charge in electrolyte phase:

$${\kappa }^{eff}\frac{{\partial }^2{\varphi }_e}{\partial x^2}+{\kappa }_d^{eff}\frac{{\partial }^2\ln c_e}{\partial x^2}+j^{\mathrm{Li}}=0,x\in (0,L).$$

(7)

with the boundary condition:

$$\{\begin{array}{l}{\frac{\partial {\varphi }_e}{\partial x}|}_{x=0}={\frac{\partial {\varphi }_e}{\partial x}|}_{x=L}=0,\\ {\left({\kappa }^{eff}\frac{\partial {\varphi }_e}{\partial x}+{\kappa }_d^{eff}\frac{\partial \ln c_e}{\partial x}\right)|}_{x={\delta }_n}={\left({\kappa }^{eff}\frac{\partial {\varphi }_e}{\partial x}+{\kappa }_d^{eff}\frac{\partial \ln c_e}{\partial x^2}\right)|}_{x={\delta }_n+{\delta }_{sep}}.\end{array}$$

(8)

The four PDEs are coupled by the Butler–Volmer equation [31], which describes the reaction kinetics at the solid/electrolyte interface:

$$j^{\mathrm{Li}}=a_s{i}_0\left[\exp \left(\frac{{\alpha }_{a}F}{RT}\eta \right)-\exp \left(\frac{-{\alpha }_{c}F}{RT}\eta \right)\right],\mathrm{for}x\in (0,L),$$

(9)

where symbols and their meanings are listed below in Nomenclature.

The overpotential η is defined as:

$$\eta ={\varphi }_s-{\varphi }_e-U,\mathrm{for}x\in (0,L),$$

(10)

where U is the thermodynamic equilibrium potential. Based on Equation (10), we obtain:

$$\{\begin{array}{l}{\varphi }_s(0,t)={\eta |}_{x=0}+{{\varphi }_e|}_{x=0}+{U|}_{x=0}\triangleq {\eta }_n+{\varphi }_{e,n}+U_n,\\ {\varphi }_s(L,t)={\eta |}_{x=L}+{{\varphi }_e|}_{x=L}+{U|}_{x=L}\triangleq {\eta }_p+{\varphi }_{e,p}+U_p.\end{array}$$

(11)

The voltage across the cell terminals is determined by:

$$V(t)={\varphi }_s(L,t)-{\varphi }_s(0,t)-R_{f}I(t)/A,$$

(12)

where R_f is the contact resistance.

According to Equations (11) and (12), we can rewrite the terminal voltage in Equation (12) as follows:

$$V=(U_p-U_n)+({\eta }_p-{\eta }_n)+({\varphi }_{e,p}-{\varphi }_{e,n})-R_{f}I(t)/A,$$

(13)

where U_p and U_n are the equilibrium potential of positive and negative electrode, respectively.

2.2. Simplification of Model Parameters

In this subsection, we mainly focus on the simplification of the electrochemical parameters by neglecting the non-uniform distribution of c_s, c_e, and j^Li. According to Equations (5) and (6), we obtain:

$$\{\begin{array}{l}{\delta }_n{j}_{n,avg}^{\mathrm{Li}}={\displaystyle {\int }_0^{{\delta }_n}\left({\delta }^{eff}{\partial }^2{\varphi }_{\mathrm{s}}/\partial x^2\right)}dx=I/A,\\ {\delta }_p{j}_{p,avg}^{\mathrm{Li}}=-{\displaystyle {\int }_{{\delta }_n+{\delta }_{sep}}^L\left({\delta }^{eff}{\partial }^2{\varphi }_{\mathrm{s}}/\partial x^2\right)}dx=-I/A,\end{array}$$

(14)

where $j_{n,avg}^{\mathrm{Li}}$ and $j_{p,avg}^{\mathrm{Li}}$ are the volume-averaged reaction current density in the negative and positive solid phase, respectively. Furthermore, we replace the real values of the reaction current density with the volume-average ones:

$$\{\begin{array}{l}{j}_n^{\mathrm{Li}}=j_{n,avg}^{\mathrm{Li}}=I/\left(A{\delta }_n\right),\\ j_p^{\mathrm{Li}}=j_{p,avg}^{\mathrm{Li}}=-I/\left(A{\delta }_p\right).\end{array}$$

(15)

Suppose that c_s,e is evenly distributed along the x-axes, then ${\partial }^2\ln c_e/\partial x^2=0$, and Equation (7) can be simplified as follows:

$$\{\begin{array}{l}{\mathsf{\kappa }}^{eff}{\partial }^2{\varphi }_{e,n}/\partial x_n^2+j_n^{\mathrm{Li}}=0,\\ {\mathsf{\kappa }}^{eff}{\partial }^2{\varphi }_{e,p}/\partial x_p^2+j_p^{\mathrm{Li}}=0.\end{array}$$

(16)

According to Equations (15) and (16), and the boundary conditions in Equation (8), we can obtain:

$$\{\begin{array}{l}{\varphi }_{e,n}(x_n)={\varphi }_{e,n}(0)-x_n^{2}I/\left(2A{\mathsf{\kappa }}^{eff}{\delta }_n\right),\\ {\varphi }_{e,p}(x_p)={\varphi }_{e,p}(0)+x_p^{2}I/\left(2A{\mathsf{\kappa }}^{eff}{\delta }_p\right)-x_{p}I/A{\mathsf{\kappa }}^{eff}.\end{array}$$

(17)

In Equation (17), $x_n=x$, $x_p=x-{\delta }_{sp}$. ϕ_e,n and ϕ_e,p are the electrical potential in the negative and positive electrolyte phase, respectively.

Since the electric potential is equal at both ends of the separator, it yields the following:

$${\varphi }_{e,p}(0)={\varphi }_{e,n}({\delta }_{sp})={\varphi }_{e,n}(0)-I{\delta }_{sp}^2/\left(2A{\mathsf{\kappa }}^{eff}{\delta }_n\right)\approx {\varphi }_{e,n}(0)-I{\delta }_n/\left(2A{\mathsf{\kappa }}^{eff}\right)-I{\delta }_{sep}/\left(A{\mathsf{\kappa }}^{eff}\right).$$

(18)

Based on Equations (17) and (18), we can approximate that:

$${\varphi }_{e,p}-{\varphi }_{e,n}={{\varphi }_{e,p}(x_p)|}_{x_p=L-{\delta }_{sp}}-{{\varphi }_{e,n}(x_n)|}_{x_n=0}=-I\left({\delta }_n+2{\delta }_{sep}+{\delta }_p\right)/\left(2A{\mathsf{\kappa }}^{eff}\right).$$

(19)

Furthermore, according to Equations (9) and (14), we can obtain:

$$\{\begin{array}{l}{j}_n^{\mathrm{Li}}=\frac{I}{A{\delta }_n}=a_{s,n}{i}_{0,n}\left[\exp \left(\frac{{\alpha }_{a}F}{RT}{\eta }_n\right)-\exp \left(\frac{-{\alpha }_{c}F}{RT}{\eta }_n\right)\right],\\ j_p^{\mathrm{Li}}=\frac{-I}{A{\delta }_p}=a_{s,p}{i}_{0,p}\left[\exp \left(\frac{{\alpha }_{a}F}{RT}{\eta }_p\right)-\exp \left(\frac{-{\alpha }_{c}F}{RT}{\eta }_p\right)\right],\end{array}$$

(20)

where α_a and α_c are the anodic and cathodic transfer coefficient, respectively. Assuming that α_a = α_c = α = 0.5, η_p and η_n can be approximated as:

$$\{\begin{array}{l}{\eta }_n=RT\ln \left({\xi }_n+\sqrt{{\xi }_n^2+1}\right)/\left(\alpha F\right),\\ {\eta }_p=RT\ln \left({\xi }_p+\sqrt{{\xi }_p^2+1}\right)/\left(\alpha F\right),\end{array}$$

(21)

where ${\xi }_p=-I/\left(2A{\delta }_p{a}_{s,p}{i}_{0,p}\right)$, ${\xi }_n=I/\left(2A{\delta }_n{a}_{s,n}{i}_{0,n}\right)$.

Finally, according to Equations (19) and (21), the terminal voltage in Equation (13) can be calculated by:

$$V(t)=\left(U_p(c_{s,p})-U_n(c_{s,n})\right)+\frac{RT}{\alpha F}\ln \left(\frac{{\xi }_p+\sqrt{{\xi }_p^2+1}}{{\xi }_n+\sqrt{{\xi }_n^2+1}}\right)-I\left({\delta }_n+2{\delta }_{sep}+{\delta }_p\right)/\left(2A{\mathsf{\kappa }}^{eff}\right)-R_{f}I/A.$$

(22)

In Equation (22), U_p and U_n can be considered as functions of the li-ion concentration in the solid phase.

3. Cell-level Model Reduction

The coupled governing partial differential equations (PDEs) in Equations (1) to (8) are complex and difficult to calculate, which makes their use in control-oriented applications impractical. In this section, our goal is to find a low-dimensional approximation for the nonlinear electrochemical model using cell-level mathematical theories including discretization and linearization methods.

3.1. Model Discretization

The cell terminal voltage in Equation (22) mainly describes the steady-state behavior of the cell; it can be written as a function of current I, concentration of lithium-ion in solid phase c_s,p and c_s,n. In this subsection, we aim to build the transient-state model of the cell in the state-space form. In the state-space model structure, the li-ion concentration in both the solid and electrolyte phase are modeled as the state vectors, which are highly correlated with the transient characteristics of the lithium-ion cell. Therefore, the emphasis of this subsection is on the transient modelling of the li-ion concentration.

Comparing Equation (7) with (16), we can find that Equation (16) fails to model the influence of c_e on electrical potential, which may result in the loss of model information and accuracy. Alexander et al. [32] give an approximate solution to Equation (7) that involves the effects of c_e, then Equation (19) can be approximated by:

$${\varphi }_{e,p}-{\varphi }_{e,n}=RT\beta \left(\ln c_{e,p}-\ln c_{e,n}\right)/F-I\left({\delta }_n+2{\delta }_{sep}+{\delta }_p\right)/2A{\kappa }^{eff},$$

(23)

where R is the gas constant, T is the temperature in Kelvin, F is the Faraday’s constant, and β is a constant. Consequently, the terminal voltage of the cell in Equation (22) can be rewritten as:

$$V(t)=\left(U_p(c_{s,p})-U_n(c_{s,n})\right)+\frac{RT}{\alpha F}\ln \left(\frac{{\xi }_p+\sqrt{{\xi }_p^2+1}}{{\xi }_n+\sqrt{{\xi }_n^2+1}}\right)+\frac{RT\beta \left(\ln c_{e,p}-\ln c_{e,n}\right)}{F}-\frac{I\left({\delta }_n+2{\delta }_{sep}+{\delta }_p\right)}{2A{\mathsf{\kappa }}^{eff}}-\frac{R_{f}I}{A}.$$

(24)

Equation (24) gives the output equation of the state-space model of the lithium-ion cell. The next goal is to establish the equation of state with respect to state vectors including c_s,p, c_s,n, c_e,p, and c_e,n.

Volume integration of Equation (1) yields:

$$\frac{\partial c_{s,avg}}{\partial t}=\frac{1}{V_s}{\displaystyle \int \frac{\partial c_s}{\partial t}}d{V}_s=\frac{1}{V_s}{\displaystyle \int \left[\frac{D_s}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial c_s}{\partial r}\right)\right]}d{V}_s,$$

(25)

where $V_s=4\pi R_s^3/3$, $d{V}_s=4\pi r^{2}dr$, subscript avg means volume-average. Further, substituting Equations (2) and (15) into Equation (25) yields:

$$\{\begin{array}{l}\partial c_{s,n,avg}/\partial t=-I/\left({\delta }_{n}A{\epsilon }_{s,n}F\right),\\ \partial c_{s,p,avg}/\partial t=I/\left({\delta }_{p}A{\epsilon }_{s,p}F\right),\end{array}$$

(26)

where c_s,n,avg and c_s,p,avg are the volume-averaged li-ion concentration in negative and positive solid phase, respectively.

However, input/output data are sampled in discrete form in real applications. We can rewrite Equation (26) in the discrete form as:

$$\{\begin{array}{l}{c}_{s,n,avg}(k+1)=c_{s,n,avg}(k)-\Delta tI(k)/\left({\delta }_{n}A{\epsilon }_{s,n}F\right),\\ c_{s,p,avg}(k+1)=c_{s,p,avg}(k)+\Delta tI(k)/\left({\delta }_{p}A{\epsilon }_{s,p}F\right),\end{array}$$

(27)

where $\Delta t$ is the sampling period. It can be shown that the li-ion concentration decreases in the negative solid phase during discharge (i.e., I(k) > 0), while it increases in the positive solid phase, which is consistent with the electrochemical reaction process.

Substituting Equations (4) and (15) into the volume integration of Equation (3) yields:

$$\{\begin{array}{l}\frac{\partial c_{e,n,avg}}{\partial t}=\frac{D_{e,n}^{eff}}{{\delta }_n{\epsilon }_{e,n}}\left({{\nabla }_x{c}_{e,n}|}_{x={\delta }_n}\right)+\left(\frac{(1-t_{+}^0)}{{\delta }_n{\epsilon }_{e,n}FA}\right)I,\\ \frac{\partial c_{e,p,avg}}{\partial t}=-\frac{D_{e,p}^{eff}}{{\delta }_p{\epsilon }_{e,p}}\left({{\nabla }_x{c}_{e,p}|}_{x={\delta }_{sp}}\right)-\left(\frac{(1-t_{+}^0)}{{\delta }_p{\epsilon }_{e,p}FA}\right)I.\end{array}$$

(28)

Similarly, we rewrite Equation (28) in the discrete form, and the negative electrolyte phase yields:

$$c_{e,n,avg}(k+1)=c_{e,n,avg}(k)+\left(\frac{\Delta t{D}_{e,n}^{eff}}{{\delta }_n{\epsilon }_{e,n}}\right)\left[\left({{\nabla }_x{c}_{e,n}|}_{x={\delta }_n}\right)(k)\right]+(1-t_{+}^0)\Delta tI(k)/\left({\delta }_n{\epsilon }_{e,n}FA\right),$$

(29)

where $c_{e,n,avg}$ is the volume-averaged li-ion concentration in negative electrolyte phase, ${\nabla }_x{c}_{e,n}$ is the gradient of $c_{e,n}$ in the x-direction, and ${\nabla }_x$ denotes the vector differential operator, ${\nabla }_x{c}_{e,n}=\partial c_{e,n}/\partial x$.

The positive electrolyte phase yields:

$$c_{e,p,avg}(k+1)=c_{e,p,avg}(k)-\left(\frac{\Delta t{D}_{e,p}^{eff}}{{\delta }_p{\epsilon }_{e,p}}\right)\left[\left({{\nabla }_x{c}_{e,p}|}_{x={\delta }_{sp}}\right)(k)\right]-(1-t_{+}^0)\Delta tI(k)/\left({\delta }_p{\epsilon }_{e,p}FA\right),$$

(30)

where $c_{e,p,avg}$ is the volume-averaged li-ion concentration in positive electrolyte phase and ${\nabla }_x{c}_{e,p}$ is the gradient of $c_{e,p}$ in the x-direction, ${\nabla }_x{c}_{e,p}=\partial c_{e,p}/\partial x$.

Due to the hysteresis of lithium-ion migration in the electrolyte phase, the values of ${\nabla }_x{c}_{e,n}$ and ${\nabla }_x{c}_{e,p}$ at time k + 1 are not only affected by reaction current I(k), but also related to their own values the at time k. Therefore, ${\nabla }_x{c}_{e,n}$ and ${\nabla }_x{c}_{e,p}$ can be approximated by:

$$\{\begin{array}{l}\left({{\nabla }_x{c}_{e,n}(x)|}_{x={\delta }_n}\right)(k+1)\triangleq f_1\left(\left({{\nabla }_x{c}_{e,n}(x)|}_{x={\delta }_n}\right)(k),\Delta t\right)+f_2\left(I(k)\right),\\ \left({{\nabla }_x{c}_{e,p}(x)|}_{x={\delta }_{sp}}\right)(k+1)\triangleq f_3\left(\left({{\nabla }_x{c}_{e,p}(x)|}_{x={\delta }_{sp}}\right)(k),\Delta t\right)+f_4\left(I(k)\right).\end{array}$$

(31)

Finally, according to Equations (24), (27), and (29)–(31), we can rewrite the reduced order lithium-ion battery model in discretized state-space form as:

$$\{\begin{array}{l}{x}_{\mathrm{k}+1}=f(x_{\mathrm{k}},u_{\mathrm{k}}),\\ y_{\mathrm{out},\mathrm{k}}=g(x_{\mathrm{k}},u_{\mathrm{k}}),\end{array}$$

(32)

In Equation (32), $u_{\mathrm{k}}$ is the input matrix, $u_{\mathrm{k}}=I(k)$, $y_{\mathrm{out},\mathrm{k}}$ is the output matrix, $y_{\mathrm{out},k}=V(k)$, $x_k$ is the state vector:

$$\begin{array}{ll}{x}_k& ={\left[x_1(k),x_2(k),x_3(k),x_4(k),x_5(k),x_6(k),x_7(k),x_8(k),x_9(k)\right]}^{\mathrm{T}}\\ & =[c_{s,n,avg}(k),c_{s,p,avg}(k),c_{e,n,avg}(k),c_{e,p,avg}(k),\left({{\nabla }_x{c}_{e,n,avg}(x)|}_{x={\delta }_n}\right)(k),\\ & {\left({{\nabla }_x{c}_{e,p,avg}(x)|}_{x={\delta }_{sp}}\right)(k),\left({\varphi }_{e,p}-{\varphi }_{e,n}\right)(k),\left({\eta }_p-{\eta }_n\right)(k),\left(U_p-U_n\right)(k)]}^{\mathrm{T}},\end{array}$$

(33)

$g(x_{\mathrm{k}},u_{\mathrm{k}})$ is the output equation:

$$\begin{array}{l}g(x_k,u_k)=\left({\varphi }_{e,p}-{\varphi }_{e,n}\right)(k)+\left({\eta }_p-{\eta }_n\right)(k+1)+\left(U_p-U_n\right)(k)-R_{f}I(k)/A\\ \begin{array}{ccc}& & \end{array} =x_7(k)+x_8(k)+x_9(k)-R_{f}I(k)/A,\end{array}$$

(34)

$f(x_{\mathrm{k}},u_{\mathrm{k}})$ is the state equation:

$$\begin{array}{l}{x}_{k+1}=f(x_k,u_k)=\\ \{\begin{array}{l}{x}_1(k+1)=c_{s,n,avg}(k+1)=c_{s,n,avg}(k)-\Delta tI(k)/\left({\delta }_{n}A{\epsilon }_{s,n}F\right)\triangleq x_1(k)-C_{1}I(k),\\ x_2(k+1)=c_{s,p,avg}(k+1)=c_{s,p,avg}(k)+\Delta tI(k)/\left({\delta }_{p}A{\epsilon }_{s,p}F\right)\triangleq x_2(k)+C_{2}I(k),\\ x_3(k+1)=c_{e,n,avg}(k+1)=c_{e,n,avg}(k)+\Delta t{D}_{e,n}^{eff}\left[\left({{\nabla }_x{c}_{e,n,avg}(x)|}_{x={\delta }_n}\right)(k)\right]/\left({\delta }_n{\epsilon }_{e,n}\right)+(1-t_{+}^0)\Delta tI(k)/\left({\delta }_n{\epsilon }_{e,n}FA\right)\\ \triangleq x_3(k)+C_3{x}_5(k)+C_{4}I(k),\\ x_4(k+1)=c_{e,p,avg}(k+1)=c_{e,p,avg}(k)-\Delta t{D}_{e,p}^{eff}\left[\left({{\nabla }_x{c}_{e,p,avg}(x)|}_{x={\delta }_{sp}}\right)(k)\right]/\left({\delta }_p{\epsilon }_{e,p}\right)-(1-t_{+}^0)\Delta tI(k)/\left({\delta }_p{\epsilon }_{e,p}FA\right)\\ \triangleq x_4(k)+C_5{x}_6(k)-C_{6}I(k),\\ x_5(k+1)=\left({{\nabla }_x{c}_{e,n,avg}(x)|}_{x={\delta }_n}\right)(k+1)\triangleq f_1\left(x_5(k)\right)+f_2(I(k)),\\ x_6(k+1)=\left({{\nabla }_x{c}_{e,p,avg}(x)|}_{x={\delta }_{sp}}\right)(k+1)\triangleq f_3\left(x_6(k)\right)+f_4(I(k)),\\ x_7(k+1)=\left({\varphi }_{e,p}-{\varphi }_{e,n}\right)(k+1)=RT\beta \left(\ln c_{e,p,avg}(k)-\ln c_{e,n,avg}(k)\right)/F-I(k)\left({\delta }_n+2{\delta }_{sep}+{\delta }_p\right)/\left(2A{\mathsf{\kappa }}^{eff}\right)\\ \triangleq f_5\left(x_4(k)\right)-f_6\left(x_3(k)\right)-C_{7}I(k),\\ x_8(k+1)=\left({\eta }_p-{\eta }_n\right)(k+1)=RT\ln \left[\left({\xi }_p+\sqrt{{\xi }_p^2+1}\right)/\left({\xi }_n+\sqrt{{\xi }_n^2+1}\right)\right]/\left(\alpha F\right)\triangleq f_7\left(I(k)\right),\\ x_9(k+1)=\left(U_p-U_n\right)(k+1)=U_p\left(c_{s,p,avg}(k)\right)-U_n(c_{s,n,avg})\triangleq f_8\left(x_2(k)\right)-f_9\left(x_1(k)\right).\end{array}\end{array}$$

(35)

where, C₁, C₂, C₃, C₄, C₅, C₆, and C₇ are all constants.

3.2. Model Linearization

The proposed reduced-order model in Equation (32) is still extraordinarily complex due to the involved nonlinear equations. In this subsection, our goal is to linearize the reduced-order battery model in Equation (32). As such, we employ the small-signal analysis method to remove the nonlinearities of the battery model. The basic principle of the small-signal analysis method is to expand a nonlinear function into Taylor series, and retain only the constant and the first-order term of the series, which lead to a linear approximation of the nonlinear equation. More specifically, a nonlinear function y = f(x) can be approximated by:

$$y=f(x)\approx f(\overline{x})+{df/dx|}_{x=\overline{x}}(x-\overline{x}),$$

(36)

where $\overline{x}$ denotes the equilibrium working point and $(x-\overline{x})$ represents a small variation in x around $\overline{x}$ with $(x-\overline{x})\ll 1$.

Similarly, for a nonlinear system $y=f(x_1,x_2)$ which consists of two input variables x₁ and x₂, a linear approximation can be given by:

$$\{\begin{array}{l}y=f(x)\approx f({\overline{x}}_1,{\overline{x}}_2)+K_1(x_1-{\overline{x}}_1)+K_2(x_2-{\overline{x}}_2),\\ K_1={\partial f/\partial x_1|}_{x_1={\overline{x}}_1},\begin{array}{c}\end{array}{K}_2={\partial f/\partial x_2|}_{x_2={\overline{x}}_2}.\end{array}$$

(37)

A similar linearization process can be employed when the number of variables is larger than two.

The linearization of function $f_1\left(x_5(k)\right)+f_2(I(k))$ using small-signal analysis method yields:

$$x_5(k+1)=y=\overline{y}+K_1\left(x_5(k)-{\overline{x}}_5(k)\right)+K_2\left(I(k)-\overline{I}(k)\right),$$

(38)

When a lithium-ion cell operates at the equilibrium point, the reaction current j = I = 0. In this case, li-ions are stored in the solid electrode with very few remaining in the electrolyte; therefore, we can approximate that $c_{e,n,avg}=0$, so that $\overline{y}={\overline{x}}_5(k+1)={\overline{x}}_5(k)=0$. Consequently, Equation (38) can be simplified as:

$$x_5(k+1)=K_1{x}_5(k)+K_{2}I(k),$$

(39)

where $K_1={d{f}_1\left(x_5(k)\right)/d{x}_5(k)|}_{x_5(k)}$,$K_2={d{f}_2(I(k))/dI(k)|}_{I(k)=0}$.

Similar linearization results can be obtained for the nonlinear function with zero equilibrium variables (i.e., f₃, f₄, f₅, f₆, f₇).

For the functions with non-zero variables (i.e., f₈, f₉), the linearization process yields:

$$x_9(k+1)=\overline{y}+K_9\left(x_2(k)-{\overline{x}}_2(k)\right)+K_{10}\left(x_1(k)-{\overline{x}}_1(k)\right)=K_9{x}_2(k)+K_{10}{x}_1(k)+E(k),$$

(40)

with:

$$\{\begin{array}{l}{K}_9={d{f}_8\left(x_2(k)\right)/d{x}_2(k)|}_{x_2(k)={\overline{c}}_{s,p,avg}(k)},\\ K_{10}={-d{f}_9\left(x_1(k)\right)/d{x}_1(k)|}_{x_1(k)={\overline{c}}_{s,n,avg}(k)},\\ E(k)={\overline{U}}_p-{\overline{U}}_n-K_9{\overline{c}}_{s,p,avg}(k)-K_{10}{\overline{c}}_{s,n,avg}(k),\end{array}$$

(41)

where ${\overline{c}}_{s,p,avg}(k)$ and ${\overline{c}}_{s,n,avg}$ are the volume-averaged li-ion concentration at the equilibrium point, and they are non-zero constant. E(k) is a non-zero constant which denotes the equilibrium inner potential.

Finally, the battery model in Equation (32) can be linearized as a 9th order model:

$$\{\begin{array}{l}{x}_{\mathrm{k}+1}=A{x}_{\mathrm{k}}+B{u}_{\mathrm{k}},\\ y_{\mathrm{k}}=C{x}_{\mathrm{k}}+D{u}_{\mathrm{k}},\end{array}$$

(42)

where u_k is the input current sequence, y_k is the output voltage sequence, $y_k=V(k)-E(k)$, and A is the system matrix consists of constant:

$$A=\left[\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& C_3& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& C_5& 0& 0& 0\\ 0& 0& 0& 0& K_1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& K_3& 0& 0& 0\\ 0& 0& K_6& K_5& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ K_{10}& K_9& 0& 0& 0& 0& 0& 0& 0\end{array}\right],$$

(43)

B is the input matrix consists of constant:

$$B={\left[\begin{array}{ccccccccc}-C_1& C_2& C_4& -C_6& K_2& K_4& -C_7& K_8& 0\end{array}\right]}^{\mathrm{T}},$$

(44)

C is the observation matrix:

$$C=\left[\begin{array}{ccccccccc}0& 0& 0& 0& 0& 0& 1& 1& 1\end{array}\right],$$

(45)

D is the feed-through matrix:

$$D=\left[-R_f/A\right].$$

(46)

4. System-Level Model Identification

The linearized 9th order model in Equation (42) shows that there are numerous unknown parameters that need to be predetermined. Theoretically, the model parameters (i.e., the values of A, B, C, D) can be obtained from battery tests. However, the determination of these parameters is a non-trivial task. It is always the case that one must make empirical guesses about many parameters, which makes the model unsatisfactory. In this section, we aim to provide a complete and easily implementable parameter identification algorithm for the 9th order battery model in Equation (42) and to explore the relationship between model order and accuracy.

4.1. Subspace Identification Algorithm

Conventionally, a system is modeled by a transfer function, and is identified using optimization methods such as nonlinear least-squares. Subspace identification methods offer an alternative solution based on the state-space model structure. Subspace identification methods are non-iterative procedures and can achieve a globally optimal solution while avoid local minima problems [33], which makes the algorithm particularly suitable for the identification of complex systems like lithium-ion battery models.

Rewrite Equation (42) in the predictor form:

$$\{\begin{array}{l}{\widehat{x}}_{k+1}=\tilde{A}{\widehat{x}}_k+\tilde{B}{u}_k+K{y}_k,\\ y_k=C{\widehat{x}}_k+D{u}_k+e_k,\end{array}$$

(47)

where $\tilde{A}=A-KC$, $\tilde{B}=B-KD$. Further, we define the system order denoted by n, a past window denoted by p, a future window denoted by f, where n ≤ f ≤ p.

Based on the input and output sequence over a given time, we define the following stacked matrices:

$$\{\begin{array}{l}{\overline{y}}_{k-p,p}={[y_{k-p},y_{k-p+1},\cdots ,y_{k-1}]}^T,\\ {\overline{y}}_{k,f}={[y_k,y_{k+1},\cdots ,y_{k+f-1}]}^T,\\ Y=\left[\begin{array}{ccc}{y}_p,& \cdots ,& y_{N-1}\end{array}\right],\\ {\overline{Y}}_p=\left[\begin{array}{ccc}{\overline{y}}_{0,p},& \cdots ,& {\overline{y}}_{N-p,p}\end{array}\right].\end{array}$$

(48)

Similarly, we can define the vectors ${\overline{u}}_{k-p,p}$, ${\overline{u}}_{k,f}$, $U$, ${\overline{U}}_p$ in the same way.

Suppose that ${\tilde{A}}^j=0$ for all $j\ge p$, we rewrite the state equation of the system in Equation (47) as

$$\tilde{\Gamma }{\widehat{x}}_k=\widetilde{\Gamma L} {\overline{u}}_{k-p,p}+\widetilde{\Gamma K} {\overline{y}}_{k-p,p}.$$

(49)

$\widetilde{\Gamma L}$ and $\widetilde{\Gamma K}$ are defined as follows:

$$\widetilde{\Gamma L}=\left[\begin{array}{cccccc}{\tilde{\Xi }}^{(u_{k-p})}& {\tilde{\Xi }}^{(u_{k-p+1})}& \cdots & {\tilde{\Xi }}^{(u_{k-p+f-1})}& \cdots & {\tilde{\Xi }}^{(u_{k-1})}\\ 0& {\tilde{\Xi }}^{(u_{k-p})}& \cdots & {\tilde{\Xi }}^{(u_{k-p+f-2})}& \cdots & {\tilde{\Xi }}^{(u_{k-2})}\\ ⋮& \ddots & \ddots & ⋮& \ddots & ⋮\\ 0& \cdots & 0& {\tilde{\Xi }}^{(u_{k-p})}& \cdots & {\tilde{\Xi }}^{(u_{k-f})}\end{array}\right],$$

(50)

$$\widetilde{\Gamma K}=\left[\begin{array}{cccccc}{\tilde{\Xi }}^{(y_{k-p})}& {\tilde{\Xi }}^{(y_{k-p+1})}& \cdots & {\tilde{\Xi }}^{(y_{k-p+f-1})}& \cdots & {\tilde{\Xi }}^{(y_{k-1})}\\ 0& {\tilde{\Xi }}^{(y_{k-p})}& \cdots & {\tilde{\Xi }}^{(y_{k-p+f-2})}& \cdots & {\tilde{\Xi }}^{(y_{k-2})}\\ ⋮& \ddots & \ddots & ⋮& \ddots & ⋮\\ 0& \cdots & 0& {\tilde{\Xi }}^{(y_{k-p})}& \cdots & {\tilde{\Xi }}^{(y_{k-f})}\end{array}\right],$$

(51)

where ${\tilde{\Xi }}^{(u_{k-i})}=C{\tilde{A}}^{i-1}\tilde{B}$ for $i>0$, ${\tilde{\Xi }}^{(u_{k-i})}=D$ for $i=0$, and ${\tilde{\Xi }}^{(y_{k-i})}=C{\tilde{A}}^{i-1}K$ for $i>0$.

Similarly, with the assumption that ${\tilde{A}}^j=0$ for all $j\ge p$, the output equation of the system in Equation (47) can be approximated as:

$$y_{k|k-1,\cdots ,k-p}={\displaystyle \sum _{i=0}^p{\tilde{\Xi }}^{(u_{k-i})}}{u}_{k-i}+{\displaystyle \sum _{i=1}^p{\tilde{\Xi }}^{(y_{k-i})}}{y}_{k-i}+e_k=\tilde{\Xi }{\phi }_k+e_k.$$

(52)

where $\tilde{\Xi }\triangleq \left[{\tilde{\Xi }}^{(u_{k-p})}\cdots \begin{array}{cc}{\tilde{\Xi }}^{(u_k)}& {\tilde{\Xi }}^{(y_{k-p})}\end{array}\cdots {\tilde{\Xi }}^{(y_{k-1})}\right]$, and it is known as the Markov Matrix; ${\phi }_k$ is the input/output matrix ${\phi }_k=\left[\begin{array}{ccc}{\overline{u}}_{k-p,p}^T,& u_k^T,& {\overline{y}}_{k-p,p}^T\end{array}\right]$.

Finally, based on the above discrete and iterative processing of state-space systems, the detailed steps of predictor-based subspace identification (PBSID) algorithm [34,35] are summarized in

Table 1. Detailed scheme of the predictor-based subspace identification (PBSID) algorithm.

Step	Detailed Scheme of PBSID Algorithm
Step 1	First, sample the input and output sequence $u_k$ and $y_k$; Second, set the initial value of n, p, f, and the weight matrix W, require n ≤ f ≤ p; Third, construct the matrices Y, ${\overline{Y}}_p$, U, ${\overline{U}}_p$ given by Equation (48).
Step 2	First, construct the matrix $\Psi$, $\Psi ={\left[{\overline{U}}_p^T{U}^T{\overline{Y}}_{y-p,p}^T\right]}^T$; Second, solve the linear regression problem in Equation (52) using the least square method, then estimate the value of $\tilde{\Xi }$, $\tilde{\Xi }=Y{\Psi }^{†}$.
Step 3	First, build the matrices $\widetilde{\Gamma L}$ and $\widetilde{\Gamma K}$ given by Equations (50) and (51); Second, calculate $\begin{array}{c}\widehat{\Gamma X}\end{array}=W\left(\widetilde{\Gamma L} {\overline{U}}_p+\widetilde{\Gamma K} {\overline{Y}}_p\right)$ given by Equation (49), and estimate the state sequence through singular value decomposition (SVD) of $\begin{array}{c}\widehat{\Gamma X}\end{array}$: $$\begin{array}{c}\widehat{\Gamma X}\end{array}=\left[\begin{array}{cc}{U}_n& U_{\perp }\end{array}\right]\left[\begin{array}{cc}{\sum }_n& 0\\ 0& \sum \end{array}\right]\left[\begin{array}{c}{V}_n\\ V_{\perp }\end{array}\right];$$ Third, determine the system order n through detecting a big gap between the singular values, and calculate the state sequence: $\begin{array}{c}\widehat{X}\end{array}={\sum }_n^{1/2}{V}_n$.
Step 4	First, estimate the system matrix C and D through solving the following least squares problem: $$\underset{C,D}{\min }\left(\frac{1}{N}{\displaystyle {\sum }_{k=0}^{N-1}{\Vert y_k-C{x}_k-D{u}_k\Vert }_2^2}\right);$$ Second, estimate the system matrix A, B, and K in the same way: $$\underset{A,B,K}{\min }\left(\frac{1}{N}{\displaystyle {\sum }_{k=0}^{N-1}{\Vert {\widehat{x}}_{k+1}-\left(A{x}_k+B{u}_k+K{e}_k\right)\Vert }_2^2}\right).$$

.

4.2. Identification Experiments

The test bench in Figure 2 was used in order to identify the parameters of the battery model. It consists of a high precision battery test system for single cell test (BT2000, Arbin Instruments, College Station, TX, USA), a high voltage battery test system for the battery module test (BTS2000, Techpow Electric Co., Xiangyang, China), a thermal chamber for temperature control (HLT2005P, Hardy Technology Co., Chongqing, China), and a computer for user-machine interface and data analysis. In addition, a set of lithium iron phosphate cells were selected for testing and validation (LiFePO₄, 11Ah, Lishen Battery Co., Ltd., Tianjin, China). The detailed parameters of the LiFePO₄ cells are shown in

Table 2. Parameters of the LiFePO₄ cells for testing.

Capacity	Rated Voltage	Voltage Range	Charge Temperature	Discharge Temperature	Specific Energy Density
11 Ah	3.2 V	2.5 V–3.65 V	0–45 °C	−20 °C–60 °C	115 Wh/kg

.

Based on the 9th order battery model and the PBSID algorithm proposed above, the accuracy of the model will only be affected by the test method. In order to obtain a battery model with high accuracy, we expect that the test method can continuously excite most of the dynamic characteristics of the LiFePO₄ cells. In this research, the hybrid pulse tests at different SOC points were performed based on methods described in our previous paper [36]. The current profile of the battery identification test procedure consists of -6 A continuous discharge and hybrid pulse test, separated by 1800-s rest (Figure 3a). More specifically, -6 A continuous discharge was employed to make the battery be tested at fixed SOC points including 95%, 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, and 5%. In addition, the hybrid pulse test is a kind of pseudo-random charge and discharge sequence which was used to test the dynamic properties of the batteries. The hybrid pulse test has a wide frequency band, which can capture more detailed information of battery dynamics. After performing one set of identification test profile, we repeated it until the battery voltage reached the 2.5 V cut-off voltage. The corresponding variation in the battery voltage of the samples over time is shown in Figure 3b. The above tests were conducted at 25 °C.

4.3. Identification Results

In this subsection, the data sets of hybrid pulse tests at different SOC points were used to obtain the input and output sequence for the identification experiments. According to the 9th battery model in Equation (42), the collected current data were used as input sequence u_k, and the voltage data were used as output sequence y_k. Before identification, we set the model order n = 9, and the initial values f = p = 9 in the PBSID algorithm.

When analyzing the identification results, we used the Bode diagrams in the frequency domain instead of parametric displays in the time domain as the main analysis method. The main reason is that the identified A, B, C, and D matrices contain a large number of elements, which are not suitable for intuitive analysis and comparison. Figure 4 shows the Bode diagrams of the identified 9th battery model at different SOC points. The identification results suggest that the dynamic properties of the tested LiFePO₄ batteries vary across different SOC points and show high consistency in the SOC range of 10% to 90%. In contrast, we found significant differences in dynamic properties when the batteries were tested at 95% and 5% SOC points. Therefore, we can conclude that the dynamic characteristics of the LiFePO₄ batteries will be significantly different when the SOC is high or low.

In Figure 5, the singular values of matrix $\begin{array}{c}\widehat{\Gamma X}\end{array}$ are given using the PBSID-based model identification. As expected, the singular values of PBSID-based identification show clear gaps after the first three larger values in the SOC range of 10% to 90%. The fourth and subsequent values are less relevant during identification, since their contributions are small compared with the noise, which suggest the 3rd model may achieve a better balance between high-precision and low-complexity in a wide SOC range of 10% to 90%.

Furthermore, the singular values of matrix $\begin{array}{c}\widehat{\Gamma X}\end{array}$ at 5% and 95% SOC points are given in Figure 6. The results of the PBSID-based identification suggest that 5th order model is the preferred choice for a battery model at 5% SOC point (Figure 6a). In a similar way, we can infer that the 6th order model is a better choice for a battery model at 95% SOC point (Figure 6b).

5. Model Validation

To investigate the effectiveness and robustness of the proposed PBSID algorithm, two different model-based validation experiments were conducted in electric vehicle applications. The experiments were carried out on a new cell in the same group and included two different input current input sequences: the hybrid pulse test and the real-time UDDS driving cycle test. Our validation method is highly effective and less expensive than other options currently available.

The validation results were analyzed in time domains, which consists of the computation of tracking error and Variance-Accounted-For (VAF) on a data set that different from the data set used for determining the model. The VAF resembles the percentage of the output variation that is estimated by the model. The VAF is defined as:

$$\mathrm{VAF}\left(y_k,{\widehat{y}}_k\right)=\left(1-\frac{\mathrm{var}\left(y_k-{\widehat{y}}_k\right)}{\mathrm{var}\left(y_k\right)}\right)\times 100\%,$$

(53)

where $y_k$ is the measured output voltage, ${\widehat{y}}_k$ denotes the estimated output voltage obtained by the identified model, and var() denotes the variance of a quasi-stationary signal.

5.1. Experimental Validation

For the hybrid pulse validation experiment, the tracking error at each sampling point and the upper envelope bound was plotted to depict the model quality (Figure 7a). The corresponding SOC variation of the samples over time and the measured voltage output are also shown in Figure 7b,c, respectively. It can be observed that in the SOC range of 90% to 5%, the tracking error upper envelope varies from 0.88% to 2.99% with good robustness and high accuracy. In addition, within the entire range of SOC from 100% to 0%, the calculated value of the VAF is 96.11%, and the maximum tracking error is 6.40%. Importantly, the experiment results also show that the tracking error increases significantly within higher (100–90%) or lower (5–0%) SOC range.

The second experiment deals with the model validation in Urban Dynamometer Driving Schedule (UDDS)driving cycle. The tracking error at each sampling point was plotted to depict the model quality (Figure 8a). In addition, the current profile of the UDDS driving cycle test is shown in Figure 8b, and the corresponding measured voltage output response is also plotted in Figure 8c. It can be observed that the tracking errors remain less than 2.83% in the SOC range of 100% to 3% and increase significantly in the lower SOC range (SOC < 3%). In addition, within the entire range of SOC from 100% to 0%, the calculated value of the VAF is 94.12%, and the maximum tracking error is 5.77%. The experiment results also show that the identified 9th order battery model can achieve high-precision in the UDDS driving cycle. Therefore, we can conclude that the PBSID algorithm performs very well with high precision and good robustness under random high rate charging and discharging experimental conditions.

5.2. Statistical Analysis

The above tests show that the accuracy and dynamic characteristics of the model are related to the operating SOC range of the battery. In real-time applications, the higher the order of the battery model, the lower computational efficiency will be. To strike a balance between the model order and computation complexity, a series of model identification experiments were carried out based on the same data set of hybrid pulse test. During the identification experiments, the model order n was set to vary from 2 to 9 in the PBSID algorithm so as to obtain the nth order battery model. All the identified models were validated in the same UDDS driving cycle, and their performance was then evaluated by investigating the vector-plots of tracking errors (Figure 9). In addition, to address the question of whether the accuracy of the model depends on the operating SOC range, we used a boxplot graph to analyze the distribution of the tracking errors in the SOC range of 90% to 100% (Figure 10a), 50% to 90% (Figure 10b), 10% to 50% (Figure 10c), and 0 to 10% (Figure 10d), respectively. It can be observed that the distribution of the tracking errors vary depending on the model order and the operating SOC range. Meanwhile, the identification results demonstrate a significant variation of dynamic properties as well as model parameters across different SOC working points. The comparison plots in Figure 9 and Figure 10 reveal that the model order should be further reduced to achieve an optimal balance between high-precision and low-complexity. Therefore, it is important to choose and appropriate model order according to the operating SOC ranges in real-time applications.

As shown in Figure 10a, the boxplots are used to visually show the distribution of the tracking errors and skewness by displaying the data percentiles and averages in the SOC range of 90% to 100%. Each boxplot provides a visual summary of the tracking errors including five main percentiles (the 5%, 15%, 50%, 85%, and 95% percentile), outliers (the minimums, maximums, 1% and 99% percentile) the mean value. As expected, the 95% and 50% percentiles decrease until the 6th order model. After the 6th order model, the dispersion of the data set and signs of skewness show a high level of consistency, which suggests the 6th order battery model is the optimal choice in the SOC range of 90% to 100%.

In a similar way, we can conclude that the 3rd order battery model is the optimal choice in the wide SOC range of 10% to 90% (Figure 10b,c), and the 5th order battery model is optimal in the SOC range of 0% to 10% (Figure 10d). These results agree well with the findings of the singular values analysis in Section 4.3, which can achieve an optimal balance between high-precision and low-complexity.

6. Discussion

In this study, we showed that the lithium ion battery is a complicated non-linear electrochemical system with multi physicochemical processes. The simplified electrochemical modeling analyses suggest that the 9th order linear model structure can provide a comprehensive cell-level description of the main electrochemical reaction processes with little loss in accuracy. In addition, it has also been illustrated that the overall performance of 9th order battery model depends on not only electrolyte and electrode materials, but also on the operation condition and choice of massive physical parameters. Conversely, the calculation process of the model depends on an in-depth understanding of the thermodynamic characteristics and physical parameters of the battery, which can be costly and time-consuming with experimental methods. While not all the results were significant, the overall direction of the investigation showed trends that could be helpful to further simplifying the battery model from the electrode-level to the system-level.

Furthermore, the system-level experimental results of the LiFePO₄ battery model identification demonstrate a significant variation of dynamic properties as well as model parameters across different SOC working points. Meanwhile, the validation experiments indicated that it is important to choose the appropriate model order according to the operating SOC ranges. The boxplots of tracking errors distribution suggested that the 6th, 3rd, and 5th order battery model is the optimal choice in the SOC range of 90% to 100%, 90% to 10%, and 10% to 0%, respectively. Further data analysis showed that increasing the model order after the optimal one barely improves the accuracy, possibly due to the limited sampling precision and the additional white noise in the collected data sets. One limitation of our research is that our data only refers to one kind of lithium ion battery (LiFePO₄). Clearly, the data sets of LiFePO₄ batteries are not enough to make generalizations about all types of lithium ion batteries. However, the innovative approaches proposed in our search are generic, and the same approaches can used for conducting similar research on other types of lithium ion batteries.

In contrast to other research studies [22,23,24,25,26], we used the data-driven PBSID identification algorithm instead of physicochemical parameter measurements to obtain the model parameters. This difference mainly relates to the time-effective and application-oriented advantages involved in the PBSID algorithm: it provides the researchers with the ability to identify battery model with multiple inputs and multiple outputs and does not require any physicochemical related information.

In addition, it should be noted that the identification of the battery model was performed under random high rate charging and discharging experimental conditions with wide SOC ranges. It is expected that the accuracy of the identified models will improve considerably when applied on a more realistic operating condition, such as cycle charge and discharge with 10–90% DOD (depth of discharge) in EV and PHEV applications.

7. Conclusions

To summarize, in this paper we investigated the electrode-level modeling, cell-level model reduction, and the system-level model identification of the lithium-ion battery. This work demonstrated that the lithium-ion battery is a complicated electrochemical system with multi electrode-level physicochemical processes such as the mass and charge conservations as well as the electrochemical kinetics. We showed that it is possible to build a reduced 9th order battery model through cell-level physicochemical and mathematical theories including the volume-average analysis method and small-signal analysis method.

The system-level predictor-based subspace identification algorithm was presented and its effectiveness for the estimation of lithium-ion battery model was shown. This data-driven identification technique does not require any physicochemical processes related information, which makes the proposed modeling and identification method generic and applicable to all type of batteries. The effectiveness and robustness of the proposed methods were shown in an experimental study, where the algorithms were used on the data sets of hybrid pulse test and UDDS driving cycle test. It can be concluded that the PBSID algorithm performs very well with high precision and good robustness.

A statistical study of the tracking error distribution in different SOC ranges was conducted based on the identified battery models with different model orders. The comparison results revealed that the model order should be further reduced to achieve an optimal balance between high-precision and low-complexity. We showed that the 6th, 3rd, and 5th order battery model is the optimal choice in the SOC range of 90% to 100%, 90% to 10%, and 10% to 0%, respectively. Future work in this subject will involve the application of the proposed approaches to data from other types of lithium ion batteries.