Online Broadband Impedance Identification for Lithium-Ion Batteries Based on a Nonlinear Equivalent Circuit Model

Abstract

Models play a crucial role in explaining internal processes, estimating states, and managing lithium-ion batteries. Electrochemical models can effectively illustrate the battery’s mechanism; however, their complexity renders them unsuitable for onboard use in electric vehicles. On the other hand, equivalent circuit models (ECMs) utilize a simple set of circuit elements to simulate voltage–current characteristics. This approach is less complex and easier to implement. However, most ECMs do not currently account for the nonlinear impact of operating conditions on battery impedance, making it difficult to obtain accurate wideband impedance characteristics of the battery when used in online applications. This article delves into the intrinsic mechanism of batteries and discusses the influence of nonstationary conditions on impedance. An ECM designed for non-steady state conditions is presented. Online adaptive adjustment of model parameters is achieved using the forgetting factor recursive least squares (FFRLS) algorithm and varied parameters approach (VPA) algorithm. Experimental results demonstrate the impressive performance of the model and parameter identification method, enabling the accurate acquisition of online impedance.

1. Introduction

In recent years, electric vehicles have been considered a solution that can balance high performance, energy conservation, and environmental protection requirements. Lithium-ion batteries are widely used as vehicle power batteries due to their high energy density and long lifespan. During the operation of electric vehicles, the battery management system (BMS) is required to detect the temperature, voltage, and current of the battery and estimate the battery state to ensure normal operation. The lithium-ion battery is a complex and strongly coupled nonlinear system; moreover, ad under vehicle conditions, the battery works in an unstable state, and its internal process changes with time. Therefore, it is very difficult to estimate the state and response of the battery in real time. Establishing a suitable battery model is considered one of the key methods available to achieve this goal, and there has been considerable research in this field [1,2,3,4,5,6,7,8,9,10].

The battery model includes the electrochemical model and the equivalent circuit model (ECM). The electrochemical model focuses on the internal electrode processes of the battery, including the diffusion process of solid and liquid phases, the conduction process of solid and liquid phases, and the electrochemical process of the electrode interface. Newman et al. [11] summarized the relevant physical and chemical laws of these electrode processes and established a porous electrode model for batteries for the first time, systematically elucidating the electrode processes and mechanisms inside lithium-ion batteries. Electrochemical models are widely used in the study of aging [12], lithium plating [13,14,15], and in the growth of SEI film [16]. Although electrochemical models can comprehensively and accurately describe the internal processes of batteries, and because the parameters of the model also have clear physical meanings, the large number of partial differential equations and parameters that require identification makes it difficult to achieve online control-oriented applications; therefore, they are usually used in simulation experiments. Another commonly used battery model is the ECM. Unlike mechanistic-derived electrochemical models, the ECM uses a combination of equivalent elements such as inductance, resistance, capacitance, the constant phase element (CPE), and the Warburg element to fit the impedance of the battery and simulate its external behavior. These elements typically correspond to each electrode process; therefore, they require a wide frequency range of battery impedance for fitting. The ECM is composed of only a few equivalent elements, and the model structure is simpler, making it more suitable for control-oriented applications in vehicles. For example, Meng et al. [17] developed a battery charging control method based on the ECM, taking the battery charging efficiency and battery health status into account, and achieved good results. In the field of ECM research, Randles [18] proposed an ECM, namely the Randles model, to describe battery impedance. The model consists of resistors, capacitors, and a Warburg element. Eddine et al. [19] considered the charge transfer and diffusion processes of batteries, improved the elements in the Randles model, and established a fractional order ECM. Wang et al. [20] analyzed the electrochemical impedance spectroscopy (EIS) characteristics of lithium-ion batteries and introduced CPEs to improve the Randles model, enabling it to maintain higher accuracy. Most studies on ECMs have simplified the electrode processes to linear, with the parameters of equivalent elements being fixed values, without considering the impact of current on the electrode processes. This may lead to significant errors in the model under high current conditions, making it difficult to accurately describe the broadband impedance characteristics of the battery.

Another major challenge in the application of battery models is the identification of model parameters. In this field, the Levenberg–Marquardt (L–M) algorithm [21,22] and intelligent optimization algorithms such as the particle swarm optimization (PSO) algorithm [23,24] are commonly used. The intelligent optimization algorithm has strong universality and good applicability to the model structure and can overcome the problem of the L–M algorithm easily falling into local optima. However, the computational complexity of the intelligent optimization algorithm is large, making it difficult to apply online. Therefore, algorithms such as the deviation compensation recursive least squares algorithm and decoupling weighted recursive least squares algorithm have been proposed for the online identification of battery model parameters [25,26]. When fitting the impedance of a battery, a fractional order ECM is used to reflect the electrode process characteristics more accurately. Fractional order elements, i.e., CPEs, are used to replace ideal capacitors to describe the impedance characteristics at medium to high frequencies. Unlike integer-order models, fractional-order models require a special method to convert them from the frequency domain to the time domain for identification. Cai et al. [27] established the fractional-order ECM of the battery, achieved the discretization of the fractional-order elements in the model in the time domain according to the G-L definition of fractional-order calculus, and identified the model parameters based on the double unscented Kalman filter algorithm. Wu et al. [28] studied the fractional order model of lithium-ion batteries and identified the model parameters based on the L–M method and time-domain data. However, the fractional order was fixed to 0.5 or 1 during the identification process. Wang et al. [20] established a fractional-order ECM with a variable order; however, the fractional order needs to be determined using an offline algorithm such as PSO. Due to the nonlinearity and time-varying battery impedance, the order of the fractional-order model changes with operating conditions and states, resulting in continuous changes in the time-domain structure of the fractional-order model. Therefore, it is necessary to identify and adjust the order of the fractional order online. In existing research involving online identification methods for fractional order model parameters, the order is mostly set to a fixed value or obtained through offline methods, lacking an adaptive adjustment of the order to the operating conditions.

To establish a suitable ECM, the accurate impedance of the battery must be obtained over a wide frequency range. Cai et al. [29] used the M-sequence to measure the impedance of the battery, obtaining the impedance measurement in the frequency range of 1000 Hz~1 Hz in about 1 s. The impedance measurement speed is much faster than the traditional impedance measurement methods, and the real-time impedance information of the battery can be obtained in vehicle applications; Secondly, it is necessary to analyze the obtained impedance to establish an ECM. In existing research, the battery model does not consider the nonlinear effect of current on the electrode processes and cannot accurately describe the battery impedance characteristics under non-steady state conditions. Finally, it is necessary to identify the model parameters. In existing research, most online identification methods for models do not consider the adaptive adjustment of the order of fractional order elements with operating conditions, as well as the identification of the nonlinear impact of current on impedance.

Considering the research gap, the main contribution of this article is as follows: (1) By combining the control equation of the basic electrode processes inside the battery, an ECM is established, considering the influence of current under non-steady state operating conditions; (2) According to the definition of fractional calculus and the discretization method, the battery model is discretization. On this basis, the ECM parameters are identified, and broadband impedance is obtained using the forgetting factor recursive least squares (FFRLS) algorithm.

The remainder of this article is organized as follows: Section 2 studies the relationship between battery impedance characteristics and current bias through experiments. Section 3 establishes a battery model based on Section 2. In Section 4, the discretization of the battery model is realized, and the parameters of the model are adaptively identified based on FFRLS and the varied parameters approach (VPA) algorithms. Section 5 verifies the accuracy of the battery model and model parameter identification method through experiments. Lastly, Section 6 concludes the article.

2. Experimental Section

The experiment was conducted on the Samsung INR18650-29E batteries with a nominal capacity of 2.75 Ah. To clearly understand the impedance characteristics of the battery for modeling, broadband impedance was quickly measured, referring to the broadband impedance fast measurement method based on M-sequence excitation, as proposed in our previous work [29]. Readers can find more details there.

Figure 1 shows the measured impedance in the frequency range of 1000 Hz~1 Hz at different temperatures and SOC. As the charging rate increases, the arc in the mid-frequency region of the impedance spectrum gradually contracts, and the same pattern occurs during discharge. Comparing Figure 1c,e,f, it can be seen that, the lower the temperature is, the greater the impedance will be; moreover, as the charging rate increases, the more obvious the contraction of the arc in the mid-frequency range of the impedance spectrum becomes.

3. Modeling

To explain and describe the above experimental results, an impedance model of the battery is derived from the electrode process mechanism. This model starts with the Butler–Volmer equation. It describes the relationship between the electrode reaction overpotential ${\eta }_i$ and the Faraday current density $j_{\mathrm{fd},i}$, as shown in Equation (1).

$$j_{\mathrm{fd},i}=j_{0,i}\left[\exp \left(\frac{\alpha \mathrm{F}}{\mathrm{R}T}{\eta }_i\right)-\exp \left(\frac{-\left(1-\alpha \right)\mathrm{F}}{\mathrm{R}T}{\eta }_i\right)\right]$$

(1)

where $j_{0,i}$ is the exchange current density; $\mathrm{R}$ is the universal gas constant; $\mathrm{F}$ is the Faraday constant; $\alpha$ is the electrode reaction transfer coefficient; $T$ is the temperature.

For the Faraday current at the electrode interface participating in electrochemical reactions, considering the influence of reaction overpotential and reactant concentration, the Faraday current density can be expressed as Equation (2).

$$Z_{\mathrm{fd},i}=R_{\mathrm{ct},i}-R_{\mathrm{ct},i}\left(\frac{\partial j_{\mathrm{fd},i}}{\partial c_{\mathrm{s},\mathrm{surf},i}}\right)\frac{\Delta c_{\mathrm{s},\mathrm{surf},i}}{\Delta j_{\mathrm{fd},i}}=R_{\mathrm{ct},i}+Z_{\mathrm{d},i}$$

(2)

From the above equation, the Faraday impedance mainly consists of two parts: the charge transfer resistance $R_{\mathrm{ct},i}$ and diffusion impedance $Z_{\mathrm{d},i}$. This article mainly analyzes the charge transfer resistance under non-steady state conditions.

Under quasi-steady state conditions, there is no DC bias, and the amplitude of the excitation current is very small. The exponential term in Equation (1) can be expanded using Taylor expansion near $\eta =0$ to obtain the charge transfer resistance under quasi-steady state conditions, as shown in Equation (3), with the subscript ss indicating quasi-steady state conditions.

$$R_{\mathrm{ct},\mathrm{ss},i}={\left(\frac{\partial {\eta }_i}{\partial j_{\mathrm{fd},i}}\right)}_{\mathrm{ss}}=\frac{\mathrm{R}T}{\mathrm{F}{j}_{0,i}}$$

(3)

However, when the charging and discharging current of the battery is high, that is, under non-steady state conditions, the Faraday current density at the electrode interface is high, and the first-order Taylor approximation of Equation (1) mentioned above is no longer applicable. For Equation (1), assuming that the electrode reaction transfer coefficients of the positive and negative electrodes are equal, i.e., $\alpha =0.5$ [30], according to Equation (1), Equation (4) can be obtained.

$$\frac{\mathrm{F}}{2\mathrm{R}T}{\eta }_i=\ln \left(\frac{j_{\mathrm{fd},i}}{2j_{0,i}}+\sqrt{{\left(\frac{j_{\mathrm{fd},i}}{2j_{0,i}}\right)}^2+1}\right)$$

(4)

According to the definition of charge transfer resistance $R_{\mathrm{ct},\mathrm{ss},i}$ under quasi-steady state conditions in Equation (3), the deformation sorting of Equation (4) can obtain the charge transfer resistance $R_{\mathrm{ct},\mathrm{ts},i}$ under non-steady state conditions, as shown in Equation (5).

$$R_{\mathrm{ct},\mathrm{ts},i}=\frac{{\eta }_i}{j_{\mathrm{fd},i}}=\frac{\mathrm{R}T}{\mathrm{F}{j}_{0,i}}\cdot \frac{\ln \left(\frac{j_{\mathrm{fd},i}}{2j_{0,i}}+\sqrt{{\left(\frac{j_{\mathrm{fd},i}}{2j_{0,i}}\right)}^2+1}\right)}{\frac{j_{\mathrm{fd},i}}{2j_{0,i}}}=R_{\mathrm{ct},\mathrm{ss},i}\cdot \frac{\ln \left(\frac{j_{\mathrm{fd},i}}{2j_{0,i}}+\sqrt{{\left(\frac{j_{\mathrm{fd},i}}{2j_{0,i}}\right)}^2+1}\right)}{\frac{j_{\mathrm{fd},i}}{2j_{0,i}}}$$

(5)

where the subscript ts represents non-steady state conditions, indicating that its value will change with the magnitude of the loading current.

It can be seen from Equation (5) that, unlike the quasi-steady state condition, the charge transfer resistance $R_{\mathrm{ct},\mathrm{ts},i}$ under the unsteady state condition is related to the Faraday current density $j_{\mathrm{fd},i}$ and the exchange current density $j_{0,i}$. When the Faraday current density $j_{\mathrm{fd},i}=0$, $R_{\mathrm{ct},\mathrm{ts},i}$ is equal to the charge transfer resistance $R_{\mathrm{ct},\mathrm{ts},i}$ under quasi-steady state conditions.

The battery impedance model established in this article mainly focuses on the impedance characteristics in the frequency range of 1000 Hz~1 Hz, as shown in Figure 1. This frequency range can basically cover the charge transfer process at different temperatures. In addition to Faraday impedance, the non-Faraday process at the electrode interface includes the charging and discharging processes of SEI film equivalent capacitance and electric double-layer capacitance, which are described using corresponding equivalent resistance and equivalent capacitance.

Due to the dispersion effect of porous electrodes, the characteristics of the electric double-layer capacitance in the impedance spectrum deviate from the impedance characteristics of the ideal capacitor [31]. Therefore, the impedance model established in this article replaces the ideal capacitor with a CPE, thereby more accurately describing the impedance characteristics of the electric double-layer capacitance. The expression for the impedance of a CPE, $Z_{\mathrm{dl}}$, is shown in Equation (6).

$$Z_{\mathrm{dl}}=\frac{1}{T_{\mathrm{dl}}{\left(\mathrm{j}\omega \right)}^{\alpha }}$$

(6)

where $T_{\mathrm{dl}}$ and $\alpha$ are parameters of CPEs, $\omega$ is the angular frequency. $\alpha$ is the fractional order that satisfies $0\le \alpha \le 1$.

The impedance of SEI film can be equivalent to the impedance of SEI film resistance $R_{\mathrm{film}}$ and film capacitance $C_{\mathrm{film}}$ in parallel. The impedance of the SEI film is shown as a circular arc in the high-frequency region in the impedance spectrum. However, it can be seen from the impedance spectra obtained from the battery under different states of charge and discharge in Figure 1 that the circular arc in the high-frequency region is not obvious. Moreover, unlike the charge transfer resistance, the SEI film impedance does not change with the magnitude of the loading current [32]. To reduce the complexity of the model and facilitate online application, the impedance of the SEI film is ignored in the impedance model. The impedance model established in this article under non-steady state conditions is shown in Figure 2 [20].

Since the expression of the charge transfer resistance shown in Equation (5) is related to the Faraday current density $j_{\mathrm{fd}}$, while the actual measured data are current, it is necessary to establish a relationship between current density and current magnitude. Introducing the concept of electrode active surface area $A$, the relationship between $j_{\mathrm{fd}}$ and $A$ satisfies $I_{\mathrm{R}}=j_{\mathrm{fd}}A$, where $I_{\mathrm{R}}$ is the current flowing through the branch of the charge transfer resistance. The active surface area of the electrode can be calculated using Equation (7).

$$A_i=S{L}_i{a}_{s,i}=\frac{3{\epsilon }_{\mathrm{s},i}S{L}_i}{R_i}$$

(7)

where $R_i$ is the radius of the solid phase particle; ${\epsilon }_{\mathrm{s},i}$ is the volume fraction of the solid phase component; $S$ is the cross-sectional area of the electrode; $L_i$ is the thickness of the positive or negative electrode area.

By introducing the relationship between the active surface area $A$ of the electrode and the current density $j_{\mathrm{fd}}$ into Equation (5), the expression for the influence of current on the charge transfer resistance under non-steady state conditions can be obtained [33], as shown in Equation (8).

$$R_{\mathrm{ct}}\left(I_{\mathrm{R}}\right)=\frac{\eta }{I_{\mathrm{R}}}=\frac{\mathrm{R}T}{A\mathrm{F}{j}_0}\cdot \frac{\ln \left(\frac{I_{\mathrm{R}}}{2A{j}_0}+\sqrt{{\left(\frac{I_{\mathrm{R}}}{2A{j}_0}\right)}^2+1}\right)}{\frac{I_{\mathrm{R}}}{2A{j}_0}}=R_{\mathrm{ct},0}\cdot \frac{\ln \left(k_{\mathrm{I}}{I}_{\mathrm{R}}+\sqrt{{\left(k_{\mathrm{I}}{I}_{\mathrm{R}}\right)}^2+1}\right)}{k_{\mathrm{I}}{I}_{\mathrm{R}}}$$

(8)

where $R_{\mathrm{ct},0}=\frac{\mathrm{R}T}{A\mathrm{F}{j}_0}$ is the charge transfer resistance under quasi-steady state conditions, that is, the charge transfer resistance at $I_{\mathrm{R}}=0$; $k_{\mathrm{I}}$ is the current dependence coefficient of the charge transfer resistance, which satisfies $k_I=\frac{1}{2A{j}_0}$.

In summary, the impedance expression of the impedance model shown in Figure 2 is shown in Equation (9). This model considers the nonlinear effect of the current on the charge transfer resistance and still has a high accuracy under high-rate charging and discharging conditions.

$$Z=R_0+\frac{R_{\mathrm{ct}}(I_{\mathrm{R}})}{1+R_{\mathrm{ct}}(I_{\mathrm{R}})T_{\mathrm{dl}}{s}^{\alpha }}$$

(9)

4. Methodology

4.1. Discretization of Impedance Model

The impedance model established contains CPE with fractional order characteristics. Fractional order elements cannot be discretized directly through bilinear transformation. Fractional order elements should be converted from frequency domain to time domain discrete form according to the relevant theory of fractional calculus. According to the circuit structure of the battery impedance model shown in Figure 2, the transfer function expression of impedance in the frequency domain can be obtained, as shown in Equation (10).

$$\frac{U_{\mathrm{d}}(s)-U_{\mathrm{OCV}}(s)}{I(s)}=R_0+\frac{R_{\mathrm{ct}}(I_{\mathrm{R}})}{1+R_{\mathrm{ct}}(I_{\mathrm{R}})T_{\mathrm{dl}}{s}^{\alpha }}$$

(10)

where $U_{\mathrm{OCV}}$ is the open circuit voltage of the battery; $U_{\mathrm{d}}$ is the terminal voltage of the battery; $I$ is the current of the battery; $R_0$ is the ohmic internal resistance; $R_{\mathrm{ct}}\left(I_{\mathrm{R}}\right)$ is the charge transfer resistance.

According to Equation (8), because the relationship between $R_{\mathrm{ct}}$ and the current is complex, in order to facilitate the discretization of the model, it is assumed that there is a linear relationship between the overpotential and current, that is, the current is unbiased and the amplitude is small. At this time, $R_{\mathrm{ct}}\left(I_{\mathrm{R}}\right)$ degenerates to $R_{\mathrm{ct},0}$ under quasi-steady state conditions. The following algorithm first identifies $R_{\mathrm{ct},0}$ and, based on this, identifies the relationship between the current and the charge transfer resistance under non-steady state conditions; then, it corrects the quasi-steady state charge transfer resistance $R_{\mathrm{ct},0}$. The next discretization process uses $R_{\mathrm{ct},0}$ instead of $R_{\mathrm{ct}}\left(I_{\mathrm{R}}\right)$.

For the system in Figure 2, the system input is $u\left(t\right)=I\left(t\right)$, which is the battery current, and the system output is $y\left(t\right)=U_{\mathrm{d}}\left(t\right)-U_{\mathrm{OCV}}\left(t\right)$. Since the observed value of the system is the terminal voltage, the open circuit voltage is also taken as the quantity $e$ to be identified; moreover, the above relationship is incorporated into Equation (10), and its discretization processing is performed to obtain the difference equation of the recursive relationship between the battery terminal voltage and the battery current, as shown in Equation (11).

$$\begin{array}{ll}{U}_{\mathrm{d}}(k)& =a{U}_{\mathrm{d}}(k-1)+bI(k)+cI(k-1)+\left(1-a\right)e\\ & =a{U}_{\mathrm{d}}(k-1)+bI(k)+cI(k-1)+d\end{array}$$

(11)

where $d=\left(1-a\right)e$. Assuming order $\alpha$ is known, the model parameters can be obtained from the identified values $a,b,c,d$, as shown in Equation (12).

$$\left\{\begin{array}{l}{R}_0=-\frac{c}{a}\\ R_{\mathrm{ct},0}=-\frac{ab+c}{a-\alpha }+\frac{c}{a}+b\\ T_{\mathrm{dl}}=\frac{a^2{T}^{\alpha }}{\alpha c+ab\alpha }\\ e=\frac{d}{1-a}\end{array}\right.$$

(12)

where $T$ is the sampling period. As the sampling frequency of the voltage and current used in the parameter identification algorithm in this article is 10 kHz, the sampling period $T$ is 0.1 ms.

4.2. Online Identification Algorithm for Impedance Models

The FFRLS algorithm can effectively solve the problem of data saturation during the recursive process of the recursive least squares (RLS) algorithm, prevent estimation divergence, and maintain the fast convergence speed of the algorithm. The recursive calculation equations of the FFRLS algorithm are shown in Equations (13)–(15).

$$K_k=\frac{P_{k-1}{\phi }_k}{\lambda +{\phi }_k^T{P}_{k-1}{\phi }_k}$$

(13)

$${\widehat{\theta }}_k={\widehat{\theta }}_{k-1}+K_k(y_k-{\phi }_k^T{\widehat{\theta }}_{k-1})$$

(14)

$$P_k=\frac{1}{\lambda }\left(P_{k-1}-K_k{\phi }_k^T{P}_{k-1}\right)$$

(15)

where $k$ represents the time; ${\widehat{\theta }}_k$ is the estimated parameter value at time $k$; $\phi \left(k\right)$ is the data vector; $y_k-{\phi }_k^T{\widehat{\theta }}_{k-1}$ is the deviation caused by the parameter estimation value at time $k-1$ used for $k$-time prediction; $K_k$ is the gain matrix; $P_k$ is the covariance matrix; forgetting factor $\lambda$ indicates the speed of forgetting old data, with a general value range of $0.95\le \lambda \le 1$.

A smaller $\lambda$ causes a faster forgetting speed, which accelerates the convergence speed of the algorithm; however, if $\lambda$ is too small, this can lead to the poor stability of the algorithm. The closer $\lambda$ is to 1, the higher the stability of the algorithm becomes. When $\lambda =1$, it indicates that the old data has not been forgotten, and the FFRLS algorithm is equivalent to the recursive least squares algorithm.

Combining the discretization results of the impedance model shown in Equation (11), the data vector can be obtained as $\phi \left(k\right)={\left[U_{\mathrm{d}}\left(k-1\right),I\left(k\right),I\left(k-1\right),1\right]}^T$, the model parameter vector to be identified is $\theta ={\left[R_0,R_{\mathrm{ct},0},T_{\mathrm{dl}},e\right]}^T$, the parameter vector used for the FFRLS algorithm recursive process is ${\theta }_{\mathrm{d}}={\left[a,b,c,d\right]}^T$, and the relation between ${\theta }_{\mathrm{d}}$ and $\theta$ is shown in Equation (12). The parameter identification process based on the FFRLS algorithm is shown in Figure 3.

In the parameter identification process based on the FFRLS algorithm mentioned above, it is assumed that the order $\alpha$ of the CPE is known. However, the order of the CPE will change with operating conditions and battery status; therefore, an adaptive identification of the order is necessary. When the order of the CPE changes, the battery terminal voltage can be expressed in the form of Equation (16).

$$\begin{array}{ll}{\widehat{U}}_{\mathrm{d}}(\alpha )=& \frac{\alpha R_{\mathrm{ct},0}{T}_{\mathrm{dl}}}{R_{\mathrm{ct},0}{T}_{\mathrm{dl}}+T^{\alpha }}{U}_{\mathrm{d}}(k-1)+\frac{\left(R_0+R_{\mathrm{ct},0}\right)T^{\alpha }+R_0{R}_{\mathrm{ct},0}{T}_{\mathrm{dl}}}{R_{\mathrm{ct},0}{T}_{\mathrm{dl}}+T^{\alpha }}I(k)\\ & -\frac{\alpha R_0{R}_{\mathrm{ct},0}{T}_{\mathrm{dl}}}{R_{\mathrm{ct},0}{T}_{\mathrm{dl}}+T^{\alpha }}I(k-1)+\frac{(1-\alpha )R_{\mathrm{ct},0}{T}_{\mathrm{dl}}+T^{\alpha }}{R_{\mathrm{ct},0}{T}_{\mathrm{dl}}+T^{\alpha }}e\end{array}$$

(16)

Based on this, the objective function $F\left(\alpha \right)=\min \left|U_{\mathrm{d}}\left(k\right)-{\widehat{U}}_{\mathrm{d}}\left(\alpha \right)\right|$ can be constructed. When the $F\left(\alpha \right)$ obtained at a certain value of order $\alpha$ is the minimum value of all repeated processes at time $k$, the order value is the optimal identification result at time $k$ and is used in the recursive iteration at the next time. Set the optimal order obtained from the previous moment as the optimal order obtained from the current moment as the initial value ${\alpha }_k$ for the current order optimization. The number of trial and error attempts for the order is $N$, and the step size for each trial and error is $\Delta \alpha$. Set ${\alpha }_k$ as the center value of order $\alpha$ during the repeated verification process at the current moment, and the range of change for $\alpha$ at the current moment is ${\alpha }_k-\frac{N\Delta \alpha }{2}\le \alpha \le {\alpha }_k+\frac{N\Delta \alpha }{2}$. The larger the number of trial and error attempts $N$, the smaller the step size $\Delta \alpha$ for each trial and error is, leading to more accurate results of the order in the identification; however, the amount of computation will also increase accordingly. This article sets the initial values of $\Delta \alpha =0.01$, $N=20$, and $\alpha =0.8$.

Due to the small sampling period of voltage and current in the parameter identification algorithm in this article, executing the order trial and error process at each sampling time will greatly increase the computational workload. Therefore, the minimum time interval for performing order adjustment is set to 1 s, and a terminal voltage error threshold is set. The order identification process is executed only when the terminal voltage error exceeds the threshold at a certain time. The process of the order adaptive identification algorithm is shown in Figure 4.

Under the quasi-steady state condition of an unbiased small current, the overpotential and current exhibit an approximately linear relationship. Under this premise, the identification of the transfer resistance using the FFRLS algorithm is accurate. When the charging and discharging current of the battery is large, the overpotential and current exhibit a complex nonlinear relationship, and the charge transfer resistance is affected by the current. Therefore, it is necessary to correct $R_{\mathrm{ct},0}$ according to Equation (8). At this time, it is necessary to identify the current dependence coefficient $k_{\mathrm{I}}$ of the charge transfer resistance. Due to $k_{\mathrm{I}}=\frac{1}{2A{j}_0}$, determined by the exchange current density and the active surface area of the electrode, both of which change with the aging state and temperature of the battery, the parameter $k_{\mathrm{I}}$ must adapt to changes in battery states. Therefore, a first in first out (FIFO) data queue is needed to collect eligible current and voltage data online. After the queue has collected enough data, parameter optimization algorithms are executed on the voltage and current data in the queue. The optimization goal is to minimize the root mean square error between the measured terminal voltage and the model terminal voltage, thereby obtaining the value of parameter $k_{\mathrm{I}}$.

In Equation (8), the current flowing through the charge transfer resistance needs to be obtained. However, according to the connection relationship between the charge transfer resistance and the electric double-layer capacitance described in Figure 5, the obtained battery current flows through two branches. During battery charging and discharging, the electric double-layer capacitance also undergoes corresponding charging and discharging processes, and the current flowing through the equivalent capacitor $I_{\mathrm{dl}}\left(t\right)$ also constantly changes. When the capacitor is fully charged or discharged, $I_{\mathrm{dl}}\left(t\right)$ is close to 0. At this point, all battery current flows through the charge transfer resistance, i.e., $I_{\mathrm{R}}\left(t\right)=I\left(t\right)$, and the battery current can be used for parameter identification.

Using Simulink to simulate the circuit model shown in Figure 5, by changing the values of $R_{\mathrm{ct}}$ and $C_{\mathrm{dl}}$ for multiple sets of simulations, it was found that the time required for $I_{\mathrm{dl}}$ to decrease to $1/10,000$ of $I$ is related to the time constant of the RC circuit, which is about seven times the time constant. Therefore, it is believed that, after 1 s of charging current loading, the capacitor is fully charged and all battery current flows through the charge transfer resistance.

The algorithm designed in this article judges the current data at each moment. When the charging/discharging rate is greater than or equal to 0.1C (C-rate) and the duration exceeds 1 s, the data meets the conditions for entering the queue; Continue to make judgments on the following data. When the charge/discharge ratio is greater than or equal to 0.2C, store the measured current and voltage data at the current time, as well as the model parameters identified by the FFRLS algorithm, in the queue. When the current direction at a certain moment is opposite to the current direction of the data in the queue, the data at that moment is discarded, and the queue condition count is set to 0. Starting from the next moment, the queue condition judgment is re-executed. Therefore, the data in the queue is obtained by continuing to charge or discharge while ensuring that the electric double-layer capacitor is fully charged or discharged; moreover, the charging or discharging rate is not less than 0.2C, which can be used to identify the nonlinear relationship between the charge transfer resistance and current. A parameter optimization algorithm is executed when the amount of data in the queue reaches the set value. While executing the algorithm, data that meet the requirements still enter the queue from the end of the queue. If the amount of data in the queue exceeds the capacity of the queue, the earliest data will be discarded from the beginning of the queue.

After sufficient data are stored in the queue, the VPA algorithm is used to identify the current dependence coefficient $k_{\mathrm{I}}$ of the charge transfer resistance. VPA is a pattern search algorithm. Waag et al. [34] first proposed VPA as a model parameter optimization method and applied the parameter identification results to battery diagnosis. Moreover, the FFRLS algorithm has already estimated the parameter vector $\theta ={\left[R_0,R_{\mathrm{ct},0},T_{\mathrm{dl}},e\right]}^T$ of the model at each moment, based on which the VPA only needs to identify $k_I$ as a parameter. For each sampling time in the queue, the transfer resistance $R_{ct,K_I}$ under non-steady state is calculated according to Equation (17).

$$R_{\mathrm{ct},k_{\mathrm{I}}}=R_{\mathrm{ct},0}\cdot \frac{\ln \left(k_{\mathrm{I}}I+\sqrt{{\left(k_{\mathrm{I}}I\right)}^2+1}\right)}{k_{\mathrm{I}}I}$$

(17)

Since the data in the queue are obtained after the equivalent electric double-layer capacitor is fully charged or emptied, it meets the requirements of $I_{\mathrm{dl}}\left(t\right)=0$. Therefore, the predicted terminal voltage $U_{\mathrm{pred}}$ of the model at each time can be calculated using Equation (18).

$$U_{\mathrm{pred}}=\left(R_0+R_{\mathrm{ct},k_{\mathrm{I}}}\right)I(t)+e$$

(18)

The root mean square error (RMSE) of the measurement terminal voltage $U_{\mathrm{meas}}$ and the model prediction terminal voltage $U_{\mathrm{pred}}$ of the data in the queue are selected as the cost functions for parameter optimization, as shown in Equation (19), where $N$ is the number of data in the queue. The optimization goal of the algorithm is to minimize the cost function.

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum _{k=1}^N{\left(U_{\mathrm{pred}}-U_{\mathrm{meas}}\right)}^2}}$$

(19)

In each iteration of VPA, there are three variations in the value of $k_{\mathrm{I}}$ based on the optimal parameter $k_{\mathrm{I},\mathrm{ref}}$ obtained from the previous iteration: $k_{\mathrm{I},\mathrm{ref}}$, $\left(1+k_{\mathrm{vary}}\right)k_{\mathrm{I},\mathrm{ref}}$, $\left(1-k_{\mathrm{vary}}\right)k_{\mathrm{I},\mathrm{ref}}$. Where $0<k_{\mathrm{vary}}<1$, the variation range of parameter $k_{\mathrm{I}}$ is determined, and setting $k_{\mathrm{vary}}=0.8$ in the initial stage can obtain a relatively wide search range. For the three variants of $k_{\mathrm{I}}$, the RMSE of the terminal voltage in the queue is calculated based on Equations (17)–(19). The optimal parameter is the $k_{\mathrm{I}}$ variation form that minimizes the RMSE, which serves as the reference value for the next iteration. If the optimal parameter obtained in a certain iteration is still $k_{\mathrm{I},\mathrm{ref}}$, it indicates that the search range of $k_{\mathrm{I}}$ is too large. Therefore, in the next iteration, halve $k_{\mathrm{vary}}$ to narrow the range of $k_{\mathrm{I}}$ changes. Repeat the above iterative process until the terminal voltage RMSE meets the preset accuracy or reaches the maximum number of iterations. The maximum number of iterations depends on the processor’s computing power and allowed computing time. The larger the number of iterations, the more accurate the estimation result of $k_{\mathrm{I}}$, ensuring that the algorithm can converge to the accurate value of $k_{\mathrm{I}}$. Figure 6 shows the convergence of $k_{\mathrm{I}}$ during the execution of a VPA algorithm, where the initial value is $k_{\mathrm{I},0}=0.01$. From Figure 6, it can be seen that the optimal value of $k_{\mathrm{I}}$ is about 0.0205; moreover, the points near the optimal value are relatively dense, indicating that the algorithm gradually decreases $k_{\mathrm{vary}}$ during the iteration process, the search range gradually approaches the optimal value, and the search range gradually decreases.

In summary, the identification process for the current dependence coefficient $k_{\mathrm{I}}$ of the charge transfer resistance is as follows.

Step 1: Initialize the capacity of the queue $M=2000$, the initial value of $k_{\mathrm{I}}$, $k_{\mathrm{I},0}=0.01$, the parameter determining the range of changes of $k_{\mathrm{I}}$: $k_{\mathrm{vary}}=0.8$.

Step 2: After executing the FFRLS recursive process at each sampling time, judge the current data: if the charging/discharging ratio is greater than or equal to 0.1C, add one to the count. Due to the sampling frequency $f_{\mathrm{s}}=10\mathrm{kHz}$, when the count exceeds 10,000, it indicates that it has been continuously charged or discharged for 1 s and that the following data have met the admission conditions.

Step 3: Continue to assess the current data. If the charge/discharge ratio is greater than or equal to 0.2C, store the measured current, terminal voltage, and model parameter $\theta ={\left[R_0,R_{\mathrm{ct},0},T_{\mathrm{dl}},e\right]}^T$ estimated using the FFRLS algorithm at the current time in the queue.

Step 4: Start executing the $k_{\mathrm{I}}$ identification algorithm when 1000 sets of data have been stored in the queue. Firstly, based on the initial values of $R_{\mathrm{ct},0}$ and $k_{\mathrm{I}}$ estimated using the FFRLS algorithm at each moment in the queue, the corresponding $R_{\mathrm{ct},k_{\mathrm{I}}}$ is calculated according to Equation (17). Calculate the RMSE of the terminal voltage in the current queue according to Equations (18) and (19).

Step 5: Set the initial value of $k_{\mathrm{I}}$ as the reference value $k_{\mathrm{I},\mathrm{ref}}$ for the first iteration, so that $k_{\mathrm{I}}$ is $k_{\mathrm{I},\mathrm{ref}}$, $\left(1+k_{\mathrm{vary}}\right)k_{\mathrm{I},\mathrm{ref}}$, $\left(1-k_{\mathrm{vary}}\right)k_{\mathrm{I},\mathrm{ref}}$, and calculate the RMSE of the terminal voltage of the data in the current queue according to Equations (17)–(19).

Step 6: Among the three variants of $k_{\mathrm{I}}$ relative to the reference value, the $k_{\mathrm{I}}$ value that minimizes the terminal voltage RMSE is the optimal parameter $k_{\mathrm{I},\mathrm{best}}$ for this iteration and is set as the reference value for the next iteration of $k_{\mathrm{I}}$. If $k_{\mathrm{I},\mathrm{best}}$ and the reference value $k_{\mathrm{I},\mathrm{ref}}$ of this iteration are the same, then $k_{\mathrm{vary}}$ will be halved in the next iteration.

Step 7: Repeat the iteration process of steps 5 and 6 until the terminal voltage RMSE meets the preset accuracy or reaches the maximum number of iterations. The execution of this VPA algorithm is completed, and the optimal value $k_{\mathrm{I},\mathrm{best}}$ of $k_{\mathrm{I}}$ is output, and $k_{\mathrm{I},\mathrm{best}}$ is used as the initial value of $k_{\mathrm{I}}$ for the next algorithm execution.

The identification process of the current dependence coefficient $k_{\mathrm{I}}$ of the charge transfer resistance is shown in Figure 7. At each sampling time, after executing the FFRLS recursive process, only a simple assessment of the current data needs to be made. After storing the eligible data in the queue, the next recursive process can be executed. The identification process of $k_{\mathrm{I}}$ based on data in the queue is executed in parallel with the FFRLS algorithm; therefore, the identification algorithm of $k_{\mathrm{I}}$ will not reduce the overall parameter identification speed. The identification result $k_{\mathrm{I},\mathrm{best}}$ of $k_{\mathrm{I}}$ output after each algorithm execution will be used for the correction of $R_{\mathrm{ct},0}$ from the current time until the $k_{\mathrm{I}}$ identification algorithm is executed again to obtain a new identification result $k_{\mathrm{I},\mathrm{best}}$ and to update the old value.

5. Results and Discussion

5.1. Verification of Impedance Model

Based on the battery model shown in Figure 2, commercial software Zview^® (version 3.1) was used to fit the broadband impedance measured based on M-sequence excitation when the battery in Figure 1 was charged and discharged at different rates and at different temperatures and SOC. The frequency range was 1000 Hz~1 Hz. The impedance fitting results are shown in Figure 8.

From the fitting results shown in Figure 8, it can be seen that the impedance model established in this article can accurately fit the circular arcs in the frequency band of the impedance spectrum. To further demonstrate the accuracy of the model, the RMSE between the output impedance of the model and the measured impedance is calculated as the fitting error of the model. The calculation method is shown in Equation (20)

$$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle \sum _{i=1}^n\left({\left({Z^{\prime }}_i-{z^{\prime }}_i\right)}^2+{\left({Z^{″}}_i-{z^{″}}_i\right)}^2\right)}}{n}}$$

(20)

where $n$ is the number of frequency points, ${Z^{\prime }}_i$ and ${Z^{″}}_i$ are the real and imaginary parts of the impedance calculated theoretically at a certain frequency; $z_i{}^{\prime }$ and $z_i{}^{″}$ are the real and imaginary parts of the impedance obtained experimentally at that frequency [29]. The fitting errors of the model under different battery states are listed in

Table 1. Fitting errors of the impedance model under different battery states.

Temperature	SOC	Charge/Discharge Rate	RMSE/mΩ	Temperature	SOC	Charge/Discharge Rate	RMSE/mΩ
25 °C	50%	0.3C	0.425	25 °C	70%	0.3C	0.458
		0.6C	0.478			0.6C	0.427
		−0.3C	0.450	5 °C	50%	0.3C	0.406
		−0.6C	0.493			0.6C	0.379
	30%	0.3C	0.462	35 °C		0.3C	0.495
		0.6C	0.397			0.6C	0.484

.

From Table 1, it can be seen that the fitting error between the measured impedance of the battery under different SOC and temperature states, the charging and discharging conditions at different rates, and the impedance output of the model within the same frequency range, is within 0.5 mΩ. This indicates that the impedance model has high accuracy in the frequency domain, and that the form of the model can accurately describe the impedance characteristics of the battery.

5.2. Verification of Parameter Identification Methods

The battery is charged at a rate of 0.5C (1.375 A) at 25 °C and 50% SOC. The M-sequence excitation is superimposed on the charging current, and the current excitation is continuous. Figure 9 shows the changes in the excitation current and battery terminal voltage during a certain cycle of the M-sequence during the charging process.

The measured current and voltage data in Figure 9 are used as inputs for the parameter identification algorithm. Figure 10 shows the comparison between the predicted terminal voltage and the measured terminal voltage obtained from the model parameter identification algorithm proposed in this paper. In Figure 10a, the model predicts the terminal voltage based on the FFRLS algorithm and order adaptive identification algorithm, without identifying $k_{\mathrm{I}}$ to correct the charge transfer resistance under high current conditions; in Figure 10b, the VPA algorithm is used to identify the current dependence coefficient $k_{\mathrm{I}}$ of the charge transfer resistance and to correct the charge transfer resistance under high-current conditions based on the result in Figure 10a; Figure 10c shows the comparison results of terminal voltage errors under two scenarios: (a) and (b). It can be clearly seen that the model in Figure 10b can accurately track and measure the terminal voltage, and that the terminal voltage error is smaller. The root mean square error of the terminal voltage in Figure 10a is 2.85 mV, while the root mean square error of the terminal voltage in Figure 10b is 0.228 mV. It can be seen that, after the $k_{\mathrm{I}}$ correction of the charge transfer resistance, the prediction accuracy of the terminal voltage in the parameter identification algorithm is greatly improved.

Figure 11 shows the identification results of the ohmic resistance $R_0$, transfer resistance $R_{\mathrm{ct}}$, and double-layer capacitance $T_{\mathrm{dl}}$ obtained using the above parameter identification algorithm, with the current and voltage signals shown in Figure 9 as inputs.

From Figure 11b, it can be seen that, after $k_{\mathrm{I}}$ correction, the charge transfer resistance value decreases compared to $R_{\mathrm{ct},0}$. Based on the impedance model and model parameter identification results, the impedance spectrum within a frequency range of 1000 Hz~1 Hz was reconstructed and compared with the experimental impedance. As shown in Figure 12, it can be seen that, after $R_{\mathrm{ct},0}$ is corrected by the current dependence coefficient $k_{\mathrm{I}}$ of the charge transfer resistance, the reconstructed impedance spectrum basically coincides with the measured impedance spectrum. The RMSE between the reconstructed impedance and the measured impedance in the frequency range of 1000 Hz~1 Hz is calculated. The RMSE without $k_{\mathrm{I}}$ correction is 0.729 mΩ, and the RMSE after $k_{\mathrm{I}}$ correction is 0.358 mΩ.

In order to further verify the adaptability of the parameter identification algorithm under different operating conditions, the current and voltage data of the battery during charging and discharging at different rates at different temperatures and SOC were used as input data for the identification algorithm. The output results of the model terminal voltage in the time domain and the reconstruction results of the impedance spectrum in the frequency domain were obtained. Figure 13 shows the time-domain and frequency-domain identification results of the battery when discharged at a rate of 1C (2.75 A) at 25 °C and 30% SOC. Figure 14 shows the time-domain and frequency-domain identification results of the battery when charged at a rate of 1C (2.75A) at 5 °C and 50% SOC.

From Figure 13 and Figure 14, it can be seen that, after $k_{\mathrm{I}}$ correction, the terminal voltage output by the model in the time domain can more accurately track the measured terminal voltage. In the frequency domain, the reconstruction results of the impedance spectrum can better fit and calculate the impedance spectrum. To further demonstrate the accuracy of the battery model and model parameter identification algorithm,

Table 2. RMSE of terminal voltage and impedance under different operating conditions.

Temperature	SOC	Charge/ Discharge Rate	Terminal Voltage RMSE/mV	Impedance RMSE/mΩ
25 °C	50%	0.5C	0.228	0.358
		1C	0.154	0.324
		−0.5C	0.286	0.426
		−1C	0.215	0.379
	30%	0.5C	0.239	0.361
		1C	0.173	0.412
	70%	0.5C	0.112	0.344
		1C	0.237	0.466
5 °C	50%	0.5C	0.258	0.438
		1C	0.279	0.435
35 °C	50%	0.5C	0.281	0.472
		1C	0.292	0.493

summarizes the RMSE between the measured terminal voltage and the model output terminal voltage obtained by applying the impedance model and parameter identification algorithm when the battery is charged or discharged at different rates and at different temperatures and SOC, as well as the RMSE between the reconstructed impedance and the experimental impedance within a frequency range of 1000 Hz~1 Hz based on the parameter identification results, representing the identification accuracy in the time domain and frequency domain, respectively.

From Table 2, it can be seen that, under different operating conditions, the terminal voltage RMSE obtained based on the parameter identification algorithm proposed in this article is less than 0.3 mV, and the impedance RMSE is less than 0.5 mΩ, indicating that the algorithm has high identification accuracy in both time and frequency domains under different operating conditions. In terms of model parameter identification ability for different battery cells, Lai et al. [35] achieved an accurate parameter identification of batteries in different aging states using the FFRLS algorithm. Therefore, we have reason to believe that the parameter identification method we used has good properties for different battery cells.

6. Conclusions

Starting from the mechanism of lithium-ion batteries, this article proposes an impedance model suitable for non-steady state conditions. Unlike the impedance model of batteries under quasi-steady state conditions, the impedance model established in this article considers the nonlinear characteristics of batteries. An online parameter identification algorithm for impedance models under non-steady state conditions has been proposed. This includes a preliminary identification algorithm for model parameters based on the FFRLS algorithm, an adaptive identification algorithm for the order of CPEs, as well as for the identification method used for the current dependence coefficient $k_{\mathrm{I}}$ of the charge transfer resistance based on the VPA algorithm. The established model and proposed algorithm were validated. In the time domain, and the RMSE of the output terminal voltage of the model under different operating conditions was less than 0.3 mV. In the frequency domain, the impedance spectrum reconstructed based on the identified model parameters coincides with the impedance spectrum calculated based on the M sequence, and the RMSE of the impedance under different operating conditions was less than 0.5 mΩ.