A Comparison of Battery Equivalent Circuit Model Parameter Extraction Approaches Based on Electrochemical Impedance Spectroscopy

1. Introduction

Battery Management System (BMS) is a crucial component in electric vehicles (EVs), renewable energy systems, and other applications that use rechargeable batteries to ensure safe and efficient operation. The primary function of a BMS is to monitor the battery’s state of charge (SOC) and state of health (SOH) [1]. Electrochemical impedance spectroscopy (EIS) was proposed by Heaviside 1894 [2], it has become an essential tool in the study and development of batteries, this allows researchers to better understand the fundamental electrochemical mechanisms such that correct strategies can be implemented in BMS for improving battery performance and lifespan.

EIS has been used to study the ion transport properties and electrode/electrolyte interfacial behavior of Li-ion batteries, providing insights into their performance and potential avenues for improvement [3]. In a typical EIS experiment, a small amplitude sinusoidal current/voltage signal is applied to the battery, and the resulting voltage/current response is measured over a range of frequencies [4]. The resulting impedance data can be analyzed using equivalent electrical circuit models to extract information on the underlying electrochemical behavior [5]. EIS can provide detailed information on the electrochemical processes occurring within the battery, including the charge transfer kinetics, ion transport properties, and electrode/electrolyte interfacial behavior [6]. Pastor-Fernández et al. [7] conducted the battery aging identification and quantification research by analyzing the EIS of four parallel Li-ion cells.

The EIS can be characterized using an equivalent circuit model (ECM), which represents the battery as a combination of resistive, inductive, and capacitive components; then, the ECM parameters can be identified by fitting the ECM to the measured EIS data [8]. There are different ECM models relevant to different types of batteries, this requires prior knowledge of the battery chemistry [9]; furthermore, by iteratively adjusting the ECM parameters, the best fit can be obtained.

There are various approaches to fit the ECM model to the measured EIS data. For instance, nonlinear least squares (NLS) approach can be used to estimate the parameters of a nonlinear model; this aims to optimize the nonlinear function such that the difference between experimental data and the estimations based on the ECM model can be minimized [10]. Xu et al. [11] applied a NLS approach to fit the EIS data of a supercapacitor using the ECM model, this approach can identify the ECM parameters by minimizing the error between measurements and predictions.

Besides, Genetic Algorithm (GA) is a population-based optimization approach that can also be used in fitting the ECM to EIS data; however, the computational complexity of this approach will increase significantly with the number of ECM parameters, which means that a large number of iterations are needed; furthermore, the selection of population size, mutation rate, and crossover rate requires continuous tuning to reach the optimal estimation [12].

The complex nonlinear least squares (CNLS) approach is widely used to fit the ECM model to EIS data. Pastor-Fernández et al. [7] applied the CNLS algorithm to fit ECM to the EIS data measured from four Lion-ion batteries. Feng et al. [13] applied CNLS approach to estimate ECM parameters using the EIS data collected from a battery cell at different SOC level and temperature. The drawback of this approach is that the fitting accuracy can easily be affected by the initial guess of the ECM parameters; for instance, the optimization algorithm may converge to a local minimum instead of converging to a global minimum if the initial guess is selected inappropriately; this will lead to inaccurate ECM parameter estimation. Also, CNLS approach requires the specification of ECM models, such as the number of components and the arrangement of RC circuits, which leads to extra work being done before the fitting process; furthermore, CNLS is a computationally expensive approach, especially when the size of EIS measurements are large [14].

Ghadi [15] applied the least squares (LS) approach to fit the EIS data to identify ECM parameters by assuming the solid electrolyte interphase (SEI) arc and charge transfer (CT) arc as semicircle, and the solid electrolyte interphase resistance ${\mathrm{R}}_{\mathrm{SEI}}$ and charge transfer resistance ${\mathrm{R}}_{\mathrm{CT}}$ is the diameter of the SEI arc and the CT arc, respectively; the merit of this approach is that the estimation of the parameters can be expressed in closed form; however, the main drawback is that the accuracy of this approach is not sufficient. One of the improvements is to apply the exhaustive search (ES) approach to identify more accurate estimations of ECM parameters with the assistance of LS approach, in this paper, the ES approach will be explained in detail.

Furthermore, we proposed a novel NLS approach which only needs to define the objective function, then, randomly choose the initial guess in each Monte-Carlo Run and the estimated parameters that can reach the lowest fitting error will be selected. While the ES approach can somewhat reach a better fitting accuracy compared with the LS approach, the computational time is still very slow. Besides, the computational time of the NLS approach is much faster than that of the ES approach; in addition, compared to LS and ES approach, NLS shows higher fitting accuracy.

The contributions of this paper can be summarized as follows:

• This paper compares the performance of the LS, ES, and Monte-Carlo-based NLS approaches to identifying battery ECM parameters.
• Compared to the LS approach presented in [15], the ES and the NLS approach can significantly boost the fitting accuracy of EIS measurements.
• This paper presents a novel approach to implementing NLS through Monte Carlo runs. At each Monte Carlo run, the initial parameters required for the NLS approach are selected randomly. This approach results in better accuracy and much faster computation time than the ES approach.
• All the methods are validated using both simulated EIS data with different noise levels and real EIS data collected from two different types of Li-ion batteries; the fitting performance of the NLS approach outweighs other approaches in all cases.

The remainder of the paper is organized as follows: Section 2 describes the analysis of battery ECM parameters via EIS in the frequency domain. Section 3 describes the algorithms to estimate ECM parameters using the least squares approach. Section 4 describes algorithms of exhaustive search, the Monte-Carlo-based nonlinear least squares approach is explained in Section 5. Implementation procedure is explained in Section 6. Results are discussed in Section 7. Section 8 concludes the paper.

2. Analysis of ECM Parameters in Frequency Domain

EIS is a widely used technique to investigate the impedance response of the battery. To measure the EIS, a small perturbation current with a wide range of frequencies (0.01 Hz to 10 kHz) is injected to the battery; then, by using the discrete Fourier transform (DFT), the measured voltage and current in the time domain can be converted to the frequency domain; thus, the impedance in the frequency domain can be analyzed [16,17]. The battery EIS can then be represented by the real and imaginary part of the impedance on the complex plane to form the Nyquist plot [18,19]. This plot represents the impedance spectrum of the battery at a range of frequencies; the ECM parameters can be estimated by fitting the EIS data with suitable fitting algorithms.

The frequency domain approach uses the Adaptive Randles (AR) ECM shown in Figure 1. The AR-ECM consists of the following components [17]:

• Voltage Source, $\mathrm{EMF}$
• Stray Inductance, L
• Ohmic Resistance, $R_{\mathsf{\Omega }}$
• Solid Electrolytic Interface (SEI) Resistance, $R_{SEI}$
• SEI Capacitance, $C_{SEI}$
• Charge Transfer (CT) resistance, $R_{CT}$
• Double Layer (DL) Capacitance, $C_{DL}$
• Warburg Impedance, $Z_{\mathrm{W}}$

Figure 1 shows the Nyquist plot relevant to the AR-ECM. According to this figure, the AC impedance $Z\left(w\right)$ corresponding to the AR-ECM can be written as [15]

$$\begin{array}{cc}\hfill Z\left(\omega \right)& \triangleq Z\left(j\omega \right)\hfill \\ \hfill & =j\omega L+R_{\mathsf{\Omega }}+{\displaystyle \frac{1}{\frac{1}{R_{\mathrm{SEI}}}+j\omega C_{\mathrm{SEI}}}}\hfill \\ \hfill & +{\displaystyle \frac{1}{\frac{1}{R_{\mathrm{CT}}+Z_{\mathrm{w}}\left(j\omega \right)}+j\omega C_{\mathrm{DL}}}}\hfill \\ \hfill & =\underset{Z_{\mathrm{RL}}}{\underbrace{j\omega L+R_{\mathsf{\Omega }}}}+\underset{Z_{\mathrm{SEI}}}{\underbrace{{\displaystyle \frac{R_{\mathrm{SEI}}}{1+j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}}}}}\hfill \\ \hfill & +\underset{Z_{\mathrm{CT}\&\mathrm{DF}}}{\underbrace{{\displaystyle \frac{R_{\mathrm{CT}}+Z_{\mathrm{w}}\left(j\omega \right)}{1+j\omega \left(R_{\mathrm{CT}}+Z_{\mathrm{w}}\left(j\omega \right)\right)C_{\mathrm{DL}}}}}}\hfill \end{array}$$

(1)

Where $Z_{\mathrm{RL}}$ denotes the impedance in the RL arc, $Z_{\mathrm{SEI}}$ denotes the impedance in the SEI arc, and $Z_{\mathrm{CT}\&\mathrm{DF}}$ denotes the impedance in the CT arc and DF arc.

3. Least Squares Approach

To solve the problem of ECM parameter estimation, Ghadi [15] applied the LS algorithm to fit the EIS measurements; furthermore, this approach can express the estimation of ECM parameters in closed form. In this section, an improved LS approach to AR-ECM parameter estimation is presented.

Figure 2 shows the impedance spectrum/Nyquist plot corresponding to the AR-ECM shown in Figure 1. Each data point in the Nyquist plot is obtained through the procedure as shown in Figure 3, where $z_v\left(t\right)$ and $z_c\left(t\right)$ are the measured voltage and current in the time domain while injecting sinusoidal current to the battery at different frequencies; $Z_v\left(\omega \right)$ and $Z_c\left(\omega \right)$ are the Fourier transform of the corresponding voltage and current measurements; and the real part and imaginary part of the measured impedance can be defined as:

$$\begin{array}{c}\hfill z_r\left(k\right)=\mathrm{Re}\left(Z\left({\omega }_k\right)\right),z_i\left(k\right)=-\mathrm{Im}\left(Z\left({\omega }_k\right)\right)\end{array}$$

(2)

where, ${\omega }_{k_{DF1}}\le {\omega }_k\le {\omega }_{k_{RL2}}$.

It can be observed that the Nyquist plot needs to be divided into four parts to see how it is directly related to the AR-ECM. The feature points of the Nyquist plot are indicated by index $k_{DF1}$ and $k_{DF2}$ in the DF arc; are indicated by $k_{CT1}$ and $k_{CT2}$ in the CT arc; are indicated by $k_{SEI1}$ and $k_{SEI2}$ in the SEI arc; and are indicated by $k_{RL1}$ and $k_{RL2}$ in the RL arc. Different parts of the Nyquist plot represent the battery’s impedance at different frequencies [20]. In this paper, to keep the consistency of nomenclature, we define:

• $k_{DF1}$ is the index of the first data point in the DF arc, in this paper, we define $k_{DF1}=1$.
• $k_{DF2}$ is selected such that the data points from $k_{DF1}$ to $k_{DF2}$ follows the linear line.
• $k_{CT1}$ is selected at the beginning of the CT arc such that the data points start to follow the arc.
• $k_{CT2}$ is selected at the end of the CT arc such that data points between $k_{CT1}and{k}_{CT2}$ follow the CT arc.
• Similarly, $k_{SEI1}$ is selected at the beginning of the SEI arc.
• $k_{SEI2}$ is selected at the end of the SEI arc such that data points between $k_{SEI1}and{k}_{SEI2}$ follow the SEI arc.
• $k_{RL1}$ is selected at the beginning of the RL arc.
• $k_{RL2}$ is selected at the end of the RL arc.

3.1. Estimation of Ohmic Resistance and Stray Inductance

As shown in Figure 2, based on the impedance measurements in the RL arc, the estimation of ohmic resistance $R_{\mathsf{\Omega }}$ can be estimated as

$$\begin{array}{cc}\hfill {\widehat{R}}_{\mathsf{\Omega }}& ={\displaystyle \frac{1}{k_{RL2}-k_{RL1}+1}}\sum _{k=k_{RL1}}^{k_{RL2}}{z}_r\left(k\right)\hfill \end{array}$$

(3)

and stray inductance L can be estimated using the improved method:

$$\begin{array}{cc}\hfill \widehat{L}& ={\displaystyle \frac{z_i\left(k_{RL2}\right)}{{\omega }_{k_{RL2}}}}\hfill \end{array}$$

(4)

3.2. Estimation of Diffusion Arc’s Gradient m

Considering the imaginary part of the measured impedance $z_i\left(k\right)$ and the real part of the measured impedance $z_r\left(k\right)$ in the Diffusion arc, it can be represented with a linear model [21]:

$$\begin{array}{c}\hfill z_i\left(k\right)=m{z}_r\left(k\right)+a\end{array}$$

(5)

Assuming the measurements are from $k_{DF1}$ to $k_{DF2}$ as shown in Figure 2, the following can be written as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill z_i\left(k_{DF1}\right)=& m{z}_r\left(k_{DF1}\right)+a\hfill \\ \hfill z_i\left(k_{DF1+1}\right)=& m{z}_r\left(k_{DF1+1}\right)+a\hfill \\ \hfill & ⋮\hfill \\ \hfill z_i\left(k_{DF2}\right)=& m{z}_r\left(k_{DF2}\right)+a\hfill \end{array}\end{array}$$

(6)

Equation (6) can be written in the matrix form

$$\begin{array}{c}\hfill \underset{\mathbf{y}}{\underbrace{\left[\begin{array}{c}{z}_i\left(k_{DF1}\right)\\ z_i(k_{DF1}+1)\\ ⋮\\ z_i\left(k_{DF2}\right)\end{array}\right]}}=\underset{\mathbf{A}}{\underbrace{\left[\begin{array}{cc}{z}_r\left(k_{DF1}\right)& 1\\ z_r(k_{DF1}+1)& 1\\ ⋮\\ z_r\left(k_{DF2}\right)& 1\end{array}\right]}}\underset{\mathbf{k}}{\underbrace{\left[\begin{array}{c}m\\ a\end{array}\right]}}\end{array}$$

(7)

m and a can be estimated using the LS approach:

$$\begin{array}{c}\hfill \widehat{\mathbf{k}}={\left({\mathbf{A}}^T\mathbf{A}\right)}^{-1}\left({\mathbf{A}}^T\mathbf{y}\right)\end{array}$$

(8)

$$\widehat{m}=\widehat{\mathbf{k}}\left(1\right),\widehat{a}=\widehat{\mathbf{k}}\left(2\right)$$

(9)

Algorithm 1 estimates quantities (list) based on the following impedance values:

$$\begin{array}{cc}\hfill {\mathbf{z}}_r& =[z_r\left(1\right),z_r\left(2\right),\dots ,z_r\left(n\right)]\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\mathbf{z}}_i& =[z_i\left(1\right),z_i\left(2\right),\dots ,z_i\left(n\right)]\hfill \end{array}$$

(11)

In this paper, the uppermost bound of the DF arc is denoted as $k_{DF}^H$, the lowest bound of the CT arc is denoted as $k_{CT}^L$, and the lowest bound of the SEI arc is denoted as $k_{SEI}^L$, these boundaries can be identified by applying moving average filter (MAF) to process the measured impedance data via Algorithm 1, the filtered EIS data is shown in Figure 4a. The algorithms presented in this paper are written by utilizing MATLAB syntax. Algorithm 1 uses the following MATLAB commands: ’smooth’, ’length’, ’find’,’min’, ’break’, and ’continue’.

Figure 3. The procedure to obtain battery EIS [22].

Figure 3. The procedure to obtain battery EIS [22].

Algorithm 1 Boundary identification.

The gradient m of the Diffusion arc can be estimated by fitting the Diffusion arc with the linear model mentioned in (5) and searching out the best fit using Algorithm 2. The fitting process is also demonstrated in Figure 4b,c. Algorithm 2 uses the following MATLAB commands: ’mean’, ’find’, and ’max’.

Algorithm 2 Diffusion arc fitting.

3.3. First Estimation of Warburg Coefficient $\sigma$

From the observation of EIS measurements in [20], it was found that the gradient of the Diffusion arc varies with the SOC level; in addition, gradients may be different even at the same SOC level of batteries from two different manufacturers; therefore, an improved method to represent Warburg impedance is defined mathematically as

$$\begin{array}{c}\hfill Z_w\left(j\omega \right)=(1-jm){\displaystyle \frac{\sigma }{\sqrt{\omega }}}\end{array}$$

(12)

where $\sigma$ is the Warburg coefficient, m is the gradient of fitted DF arc, and j is $\sqrt{-1}$.

It must be emphasized here that in [15] the gradient was assumed to be $m=1$. In this paper, it is proposed to estimate the gradient m to achieve better EIS fitting.

It can be shown based on (1) that the Warburg impedance is significant only at lower frequencies. In Figure 2, impedance measurements from $k_{DF1}$ to $k_{DF2}$ are selected to estimate the Warburg coefficient (we define $k_{DF1}=1$, and $k_{DF2}$ is obtained via Algorithm 2). Considering the real part of the impedance $z_r$ in the diffusion arc:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill z_r\left(k_{DF1}\right)-z_r\left(k_{DF2}\right)=\sigma & \left({\displaystyle \frac{1}{\sqrt{{\omega }_{k_{DF1}}}}}-{\displaystyle \frac{1}{\sqrt{{\omega }_{k_{DF2}}}}}\right)\hfill \\ \hfill z_r\left(2\right)-z_r\left(k_{DF2}-1\right)=\sigma & \left({\displaystyle \frac{1}{\sqrt{{\omega }_2}}}-{\displaystyle \frac{1}{\sqrt{{\omega }_{(k_{DF2}-1)}}}}\right)\hfill \\ & ⋮\hfill \\ \hfill z_r\left(n\right)-z_r\left(k_{DF2}-n+1\right)& =\sigma \left({\displaystyle \frac{1}{\sqrt{{\omega }_n}}}\right.\hfill \\ & \left.-{\displaystyle \frac{1}{\sqrt{{\omega }_{(k_{DF2}-n+1)}}}}\right)\hfill \end{array}\end{array}$$

(13)

where $n=\lfloor {\displaystyle \frac{k_{DF2}-k_{DF1}+1}{2}}\rfloor$.

The observation model corresponding to (13) is

$$\begin{array}{c}\hfill \tilde{\mathbf{z}}=\mathbf{b}\sigma \end{array}$$

(14)

where

$$\tilde{\mathbf{z}}=\left[\begin{array}{c}{z}_r\left(1\right)-z_r\left(k_{DF2}\right)\\ z_r\left(2\right)-z_r\left(k_{DF2}-1\right)\\ ⋮\\ z_r\left(n\right)-z_r\left(k_{DF2}-n+1\right)\end{array}\right],$$

$$\mathbf{b}=\left[\begin{array}{c}\left({\displaystyle \frac{1}{\sqrt{{\omega }_1}}}-{\displaystyle \frac{1}{\sqrt{{\omega }_{k_{DF2}}}}}\right)\\ \left({\displaystyle \frac{1}{\sqrt{{\omega }_2}}}-{\displaystyle \frac{1}{\sqrt{{\omega }_{(k_{DF2}-1)}}}}\right)\\ ⋮\\ \left({\displaystyle \frac{1}{\sqrt{{\omega }_n}}}-{\displaystyle \frac{1}{\sqrt{{\omega }_{(k_{DF2}-n+1)}}}}\right)\end{array}\right]$$

and the LS estimate of $\sigma$ is

$$\begin{array}{c}\hfill {\widehat{\sigma }}_{LS}={\left({\mathbf{b}}^T\mathbf{b}\right)}^{-1}\left({\mathbf{b}}^T\tilde{\mathbf{z}}\right)\end{array}$$

(15)

3.4. Estimation of $R_{\mathrm{SEI}}$ and $C_{\mathrm{SEI}}$

As shown in Figure 2, to fit the SEI arc precisely, we select feature points that lie between $k_{SEI1}$ and $k_{SEI2}$. Let us denote the impedance measurements in the SEI arc as

$$\begin{array}{cc}\hfill & s_r\left(k\right)\triangleq z_r\left(k\right)\mathrm{s}.\mathrm{t}.k_{SEI1}\le k\le k_{SEI2}\hfill \\ \hfill & s_i\left(k\right)\triangleq z_i\left(k\right)\mathrm{s}.\mathrm{t}.k_{SEI2}\le k\le k_{SEI2}\hfill \end{array}$$

(16)

The estimation of $R_{\mathrm{SEI}}$ is to fit the SEI arc using a semicircle with its centre lying on the real axis; the coordinate of this semicircle’s center can be denoted as ($x_s$, 0); the radius of the semicircle can be denoted as $r_s$; thus, the measurements in (16) should satisfy the equation of the semicircle [20]:

$$\begin{array}{c}\hfill {\left(s_r\left(k\right)-x_s\right)}^2+{\left(s_i\left(k\right)-0\right)}^2={r_s}^2\end{array}$$

(17)

$$\begin{array}{c}\hfill {s_r\left(k\right)}^2-2x_s{s}_r\left(k\right)+x_s^2+{s_i\left(k\right)}^2={r_s}^2\end{array}$$

(18)

Let $c=-2x_s$ and $d=x_s^2-{r_s}^2$, thus

$$\begin{array}{c}\hfill {r_s}^2={\displaystyle \frac{c^2}{4}}-d\end{array}$$

(19)

$$\begin{array}{c}\hfill r_s=\sqrt{{\displaystyle \frac{c^2}{4}}-d}\end{array}$$

(20)

And (18) can be rewritten as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {s_r\left(k\right)}^2+{s_i\left(k\right)}^2+c{s}_r\left(k\right)+d=0\end{array}\end{array}$$

(21)

In the matrix form, (21) can be written as

$$\begin{array}{cc}\hfill & \underset{\mathbf{z}}{\underbrace{\left[\begin{array}{c}-\left(s_r{\left(k_{SEI1}\right)}^2+s_i{\left(k_{SEI1}\right)}^2\right)\\ -\left(s_r{(k_{SEI1}+1)}^2+s_i{(k_{SEI1}+1)}^2\right)\\ -\left(s_r{(k_{SEI1}+2)}^2+s_i{(k_{SEI1}+2)}^2\right)\\ ⋮\\ -\left(s_r{\left(k_{SEI2}\right)}^2+s_i{\left(k_{SEI2}\right)}^2\right)\end{array}\right]}}\hfill \\ \hfill & =\underset{\mathbf{B}}{\underbrace{\left[\begin{array}{cc}{s}_r\left(k_{SEI1}\right)& 1\\ s_r(k_{SEI1}+1)& 1\\ s_r(k_{SEI1}+2)& 1\\ ⋮\\ s_r\left(k_{SEI2}\right)& 1\end{array}\right]}}\underset{{\mathbf{x}}_{\mathrm{SEI}}}{\underbrace{\left[\begin{array}{c}c\\ d\end{array}\right]}}+\underset{\mathbf{n}}{\underbrace{\left[\begin{array}{c}{n}_v\left(1\right)\\ n_v\left(2\right)\\ ⋮\\ n_v\left(n\right)\end{array}\right]}}\hfill \end{array}$$

(22)

Using LS algorithm, the estimate of ${\widehat{\mathbf{x}}}_{\mathrm{SEI}}$ will be given by

$$\begin{array}{c}\hfill {\widehat{\mathbf{x}}}_{\mathrm{SEI}}={\left({\mathbf{B}}^T\mathbf{B}\right)}^{-1}\left({\mathbf{B}}^T\mathbf{z}\right)\end{array}$$

(23)

The estimates of c and d are:

$$\widehat{c}={\widehat{\mathbf{x}}}_{\mathrm{SEI}}\left(1\right),\widehat{d}={\widehat{\mathbf{x}}}_{\mathrm{SEI}}\left(2\right)$$

(24)

From Figure 2, ${\mathrm{R}}_{\mathrm{SEI}}$ is the diameter of the SEI arc; thus, by substituting the values of c and d in (20), the estimate of $R_{\mathrm{SEI}}$ is

$$\begin{array}{c}\hfill {\widehat{R}}_{\mathrm{SEI}}=2{\widehat{r}}_s=2\sqrt{{\displaystyle \frac{{\widehat{c}}^2}{4}}-\widehat{d}}\end{array}$$

(25)

The estimated centre of the semicircle can then be expressed as

$$({\widehat{x}}_s,0)=(-{\displaystyle \frac{\widehat{c}}{2}},0)$$

(26)

The fitting accuracy of the SEI arc can be evaluated as [23]

$$error=\sqrt{{\displaystyle \frac{{\sum }_{k=k_{SEI1}}^{k_{SEI2}}{d_k}^2}{k_{SEI2}-k_{SEI1}+1}}}$$

(27)

Where $d_k$ is the geometrical distance between the actual EIS data point and predicted EIS data point, which is defined as

$$d_k=\sqrt{{\left[s_r\left(k\right)-{\widehat{x}}_s\right]}^2+{\left[s_i\left(k\right)-0\right]}^2}-{\widehat{r}}_s$$

(28)

It can be shown in (1) that when the frequency is very high, the impedance in CT arc and Diffusion arc will be minimal so that it is negligible; thus, we assume the $Z_{CT\&DF}$ term will be zero, that is

$$\begin{array}{c}\hfill Z=Z_{\mathrm{RL}}+Z_{\mathrm{SEI}}+0\end{array}$$

(29)

Therefore, the impedance in SEI arc can be expressed as:

$$\begin{array}{cc}\hfill Z_{\mathrm{SEI}}& =Z-Z_{\mathrm{RL}}\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\displaystyle \frac{R_{\mathrm{SEI}}}{1+j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}}}& =Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill 1+j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}& ={\displaystyle \frac{R_{\mathrm{SEI}}}{Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}}}\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}& ={\displaystyle \frac{R_{\mathrm{SEI}}}{Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}}}-1\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill j\omega C_{\mathrm{SEI}}& ={\displaystyle \frac{1}{Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}}}-{\displaystyle \frac{1}{R_{\mathrm{SEI}}}}\hfill \end{array}$$

(34)

Take the imaginary part on both sides of the above equation,

$$\begin{array}{cc}\hfill C_{\mathrm{SEI}}& =\left({\displaystyle \frac{1}{\omega }}\right)\mathrm{Im}\left({\displaystyle \frac{1}{Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}}}-{\displaystyle \frac{1}{R_{\mathrm{SEI}}}}\right)\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill C_{\mathrm{SEI}}& =\left({\displaystyle \frac{1}{\omega }}\right)\mathrm{Im}\left({\displaystyle \frac{1}{Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}}}\right)\hfill \end{array}$$

(36)

Substituting the expression for $R_{\mathsf{\Omega }}$ and L from (3) and (4), respectively, in (36) at $\omega ={\omega }_k$$(k_{SEI1}\le k\le k_{SEI2})$:

$$\begin{array}{cc}\hfill {\tilde{C}}_{\mathrm{SEI}}\left(k\right)& =\left({\displaystyle \frac{1}{{\omega }_k}}\right)\mathrm{Im}\left({\displaystyle \frac{1}{Z\left({\omega }_k\right)-j{\omega }_k\widehat{L}-{\widehat{R}}_{\mathsf{\Omega }}}}\right)\hfill \end{array}$$

(37)

Finally, average all the estimates ${\tilde{C}}_{\mathrm{SEI}}\left(k\right)$ to obtain the final estimate

$$\begin{array}{c}\hfill {\widehat{C}}_{\mathrm{SEI}}={\displaystyle \frac{1}{k_{SEI2}-k_{SEI1}+1}}\sum _{k=k_{SEI1}}^{k_{SEI2}}{\tilde{C}}_{\mathrm{SEI}}\left(k\right)\end{array}$$

(38)

Using the LS approach to identify $R_{SEI}$ and $C_{SEI}$ via the automatic selection of feature points are fully described in Algorithm 3. In addition, Figure 5a shows the RMSE of the fitted SEI arc in each iteration and Figure 5b shows the SEI arc, which is selected since it can reach the best fit. Algorithm 3 uses the following MATLAB commands: ’floor’, ’length’.

Algorithm 3 Estimate $R_{SEI}$ and $C_{SEI}$ via automatic feature detection.

3.5. Estimation of $R_{\mathrm{CT}}$ and $C_{\mathrm{DL}}$

It can be observed in Figure 2 that to fit the CT arc using a semicircle precisely, we need to select feature points that lie between $k_{CT1}$ to $k_{CT2}$; therefore, the impedance measurements in the CT arc can be denoted as:

$$\begin{array}{cc}\hfill & c_r\left(k\right)\triangleq z_r\left(k\right)\mathrm{s}.\mathrm{t}.k_{CT1}\le k\le k_{CT2}\hfill \\ \hfill & c_i\left(k\right)\triangleq z_i\left(k\right)\mathrm{s}.\mathrm{t}.k_{CT1}\le k\le k_{CT2}\hfill \end{array}$$

(39)

Assuming that the centre of the semicircle lies on the real axis, which is noted as ($x_c$, 0), the radius of the semicircle can be noted as $r_c$; thus, the measurements in (39) should satisfy the equation of the semicircle [20]:

$$\begin{array}{c}\hfill {\left(c_r\left(k\right)-x_c\right)}^2+{\left(c_i\left(k\right)-0\right)}^2={r_c}^2\end{array}$$

(40)

$$\begin{array}{c}\hfill {c_r\left(k\right)}^2-2x_c{c}_r\left(k\right)+x_c^2+{c_i\left(k\right)}^2={r_c}^2\end{array}$$

(41)

Let $g=-2x_c$ and $h=x_c^2-{r_c}^2$, thus

$$\begin{array}{c}\hfill {r_c}^2={\displaystyle \frac{g^2}{4}}-h\end{array}$$

(42)

$$\begin{array}{c}\hfill r_c=\sqrt{{\displaystyle \frac{g^2}{4}}-h}\end{array}$$

(43)

And (41) can be rewritten as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {c_r\left(k\right)}^2+{c_i\left(k\right)}^2+g{c}_r\left(k\right)+h=0\end{array}\end{array}$$

(44)

In the matrix form, (44) can be written as

$$\begin{array}{cc}\hfill & \underset{\mathbf{p}}{\underbrace{\left[\begin{array}{c}-\left(c_r{\left(k_{CT1}\right)}^2+c_i{\left(k_{CT1}\right)}^2\right)\\ -\left(c_r{(k_{CT1}+1)}^2+c_i{(k_{CT1}+1)}^2\right)\\ -\left(c_r{(k_{CT1}+2)}^2+c_i{(k_{CT1}+2)}^2\right)\\ ⋮\\ -\left(c_r{\left(k_{CT2}\right)}^2+c_i{\left(k_{CT2}\right)}^2\right)\end{array}\right]}}\hfill \\ \hfill & =\underset{\mathbf{C}}{\underbrace{\left[\begin{array}{cc}{c}_r\left(k_{CT1}\right)& 1\\ c_r(k_{CT1}+1)& 1\\ c_r(k_{CT1}+2)& 1\\ ⋮\\ c_r\left(k_{CT2}\right)& 1\end{array}\right]}}\underset{{\mathbf{x}}_{\mathrm{CT}}}{\underbrace{\left[\begin{array}{c}g\\ h\end{array}\right]}}+\underset{\mathbf{n}}{\underbrace{\left[\begin{array}{c}{n}_v\left(1\right)\\ n_v\left(2\right)\\ ⋮\\ n_v\left(n\right)\end{array}\right]}}\hfill \end{array}$$

(45)

From (45), ${\mathbf{x}}_{\mathrm{CT}}$ can be estimated using LS algorithm

$${\widehat{\mathbf{x}}}_{\mathrm{CT}}={\left({\mathbf{C}}^T\mathbf{C}\right)}^{-1}\left({\mathbf{C}}^T\mathbf{p}\right)$$

(46)

Thus, the estimates of a and b are:

$$\widehat{g}={\widehat{\mathbf{x}}}_{\mathrm{CT}}\left(1\right),\widehat{h}={\widehat{\mathbf{x}}}_{\mathrm{CT}}\left(2\right)$$

(47)

As shown in Figure 2, $R_{\mathrm{CT}}$ is the diameter of the CT arc; thus, by substituting the values of a and b in (43), the estimate of $R_{\mathrm{CT}}$ is

$$\begin{array}{c}\hfill {\widehat{R}}_{\mathrm{CT}}=2{\widehat{r}}_c=2\sqrt{{\displaystyle \frac{{\widehat{g}}^2}{4}}-\widehat{h}}\end{array}$$

(48)

The estimated centre of the semicircle can then be expressed as

$$({\widehat{x}}_c,0)=(-{\displaystyle \frac{\widehat{g}}{2}},0)$$

(49)

The fitting accuracy of the CT arc can be evaluated as [23]

$$error=\sqrt{{\displaystyle \frac{{\sum }_{k=k_{CT1}}^{k_{CT2}}{d_k}^2}{k_{CT2}-k_{CT1}+1}}}$$

(50)

Where $d_k$ is the geometrical distance between the actual EIS data point and predicted EIS data point, which is defined as [23]

$$d_k=\sqrt{{\left[c_r\left(k\right)-{\widehat{x}}_c\right]}^2+{\left[c_i\left(k\right)-0\right]}^2}-{\widehat{r}}_c$$

(51)

Based on (1),

$$\begin{array}{c}\hfill Z=Z_{\mathrm{RL}}+Z_{\mathrm{SEI}}+Z_{\mathrm{CT}\&\mathrm{DF}}\end{array}$$

(52)

Therefore, the impedance in CT arc and DF arc can be expressed as:

$$\begin{array}{cc}\hfill Z_{\mathrm{CT}\&\mathrm{DF}}& =Z-Z_{\mathrm{RL}}-Z_{\mathrm{SEI}}\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{R_{\mathrm{CT}}+Z_{\mathrm{w}}\left(j\omega \right)}{1+j\omega \left(R_{\mathrm{CT}}+Z_{\mathrm{w}}\left(j\omega \right)\right)C_{\mathrm{DL}}}}\hfill \\ & =Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}-{\displaystyle \frac{R_{\mathrm{SEI}}}{1+j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}}}\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill & j\omega \left(R_{\mathrm{CT}}+Z_{\mathrm{w}}\left(j\omega \right)\right)C_{\mathrm{DL}}\hfill \\ \hfill & ={\displaystyle \frac{R_{\mathrm{CT}}+Z_{\mathrm{w}}\left(j\omega \right)}{Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}-\frac{R_{\mathrm{SEI}}}{1+j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}}}}-1\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill j\omega C_{\mathrm{DL}}& ={\displaystyle \frac{1}{Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}-\frac{R_{\mathrm{SEI}}}{1+j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}}}}\hfill \\ \hfill & -{\displaystyle \frac{1}{R_{\mathrm{CT}}+Z_{\mathrm{w}}\left(j\omega \right)}}\hfill \end{array}$$

(56)

Take the imaginary part on both sides of (56), and substitute $Z_w\left(jw\right)$ with the expression given in (13), we obtain

$$\begin{array}{cc}\hfill & C_{\mathrm{DL}}\hfill \\ \hfill & =\left({\displaystyle \frac{1}{\omega }}\right)\mathrm{Im}\left({\displaystyle \frac{1}{Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}-\frac{R_{\mathrm{SEI}}}{1+j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}}}}\right.\hfill \\ \hfill & \left.-{\displaystyle \frac{1}{R_{\mathrm{CT}}+(1-jm)\frac{\sigma }{\sqrt{\omega }}}}\right)\hfill \end{array}$$

(57)

Substituting L, $R_{\mathsf{\Omega }}$, $R_{\mathrm{SEI}}$, $C_{\mathrm{SEI}}$, $R_{\mathrm{CT}}$ and $\sigma$ with the estimations given in (4), (3), (25), (38), (48) and (15) respectively, in the above equation at $\omega ={\omega }_k$$(k_{CT1}\le k\le k_{CT2})$:

$$\begin{array}{cc}\hfill & {\tilde{C}}_{\mathrm{DL}}\left(k\right)\hfill \\ \hfill & =\left({\displaystyle \frac{1}{{\omega }_k}}\right)\mathrm{Im}\left({\displaystyle \frac{1}{Z\left({\omega }_k\right)-j{\omega }_k\widehat{L}-{\widehat{R}}_{\mathsf{\Omega }}-\frac{{\widehat{R}}_{\mathrm{SEI}}}{1+j{\omega }_k{\widehat{R}}_{\mathrm{SEI}}{\widehat{C}}_{\mathrm{SEI}}}}}\right.\hfill \\ \hfill & \left.-{\displaystyle \frac{1}{{\widehat{R}}_{\mathrm{CT}}+(1-j\widehat{m})\frac{\widehat{\sigma }}{\sqrt{{\omega }_k}}}}\right)\hfill \end{array}$$

(58)

Finally, average all the estimates ${\tilde{C}}_{\mathrm{DL}}\left(k\right)$ to obtain the final estimate

$$\begin{array}{c}\hfill {\widehat{C}}_{\mathrm{DL}}={\displaystyle \frac{1}{k_{CT2}-k_{CT1}+1}}\sum _{k=k_{CT1}}^{k_{CT2}}{\tilde{C}}_{\mathrm{DL}}\left(k\right)\end{array}$$

(59)

Using the LS approach to identify $R_{CT}$ and $C_{DL}$ via the automatic selection of feature points is shown in Algorithm 4. Furthermore, Figure 6a shows the RMSE of the fitted CT arc in each iteration, and Figure 6b shows the CT arc selected since it can reach the best fit. Algorithm 4 uses the following MATLAB commands: ’floor’ and ’round’.

Algorithm 4 Estimate $R_{CT}$ and $C_{DL}$ via automatic feature detection.

3.6. Evaluation of the General Fitting Accuracy

In the complex plane, the absolute value of the error between the measured EIS and estimated EIS is actually the distance between measured EIS data points and estimated EIS data points as shown in Figure 7, $d_1$, where $d_2$, ... $d_n$ are the distances between measured EIS data point $(z_r\left(k\right),z_i\left(k\right))$ and estimated EIS data point $({\widehat{z}}_r\left(k\right),{\widehat{z}}_i\left(k\right))$, which can be used to evaluate the fitting accuracy, the distance $d_k$ is represented as.

$$d_k=\sqrt{{[z_r\left(k\right)-{\widehat{z}}_r\left(k\right)]}^2+{[z_i\left(k\right)-{\widehat{z}}_i\left(k\right)]}^2}$$

(60)

where n is the number of measurements, $k\in 1,2,...,n$.

Therefore, the evaluation of EIS fitting can be expressed as

$$MAE={\displaystyle \frac{1}{N}}\sum _{k=1}^N{d}_k={\displaystyle \frac{1}{N}}\sum _{k=1}^N\left|Z\left({\omega }_k\right)-\widehat{Z}\left({\omega }_k\right)\right|$$

(61)

where $Z\left({\omega }_k\right)$ is the measured impedance at ${\omega }_k$, $\widehat{Z}\left({\omega }_k\right)$ is the impedance estimation at ${\omega }_k$, which is computed based on (1) with the estimated ECM parameters, N is the total number of measurements, $|·|$ denotes the absolute value of the complex number.

The goal is to fit the EIS measurements such that the fitted EIS can achieve the lowest MAE.

The percentage error of the estimated parameters can be expressed as:

$$\begin{array}{cc}\hfill & \mathrm{Percentage}\mathrm{Error}(\%)\hfill \\ \hfill & =\left|{\displaystyle \frac{\mathrm{True}\mathrm{Value}-\mathrm{Predicted}\mathrm{Value}}{\mathrm{True}\mathrm{Value}}}\right|\times 100\%\hfill \end{array}$$

(62)

4. Exhaustive Search Approach

The Exhaustive Search (ES) approach aims to improve the goodness of fitting by searching the optimal value of each parameter based on the initial estimations of two RC components and gradient m identified via LS approach, this process is shown in Figure 8.

4.1. Second Estimation of Warburg Coefficient $\sigma$

Assume $R_{\mathsf{\Omega }}$, L, $R_{SEI}$, $C_{SEI}$, $R_{CT}$, $C_{DL}$ and m are given, based on (54),

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{\frac{1}{R_{\mathrm{CT}}+Z_{\mathrm{w}}\left(j\omega \right)}+j\omega C_{\mathrm{DL}}}}=Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}-{\displaystyle \frac{R_{\mathrm{SEI}}}{1+j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}}}\end{array}$$

(63)

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{R_{\mathrm{CT}}+Z_{\mathrm{w}}\left(j\omega \right)}}+j\omega C_{\mathrm{DL}}={\displaystyle \frac{1}{Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}-\frac{R_{\mathrm{SEI}}}{1+j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}}}}\end{array}$$

(64)

thus,

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{R_{\mathrm{CT}}+Z_{\mathrm{w}}\left(j\omega \right)}}={\displaystyle \frac{1}{Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}-\frac{R_{\mathrm{SEI}}}{1+j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}}}}-j\omega C_{\mathrm{DL}}\end{array}$$

(65)

then,

$$\begin{array}{c}\hfill R_{\mathrm{CT}}+Z_{\mathrm{w}}\left(j\omega \right)={\displaystyle \frac{1}{\frac{1}{Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}-\frac{R_{\mathrm{SEI}}}{1+j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}}}-j\omega C_{\mathrm{DL}}}}\end{array}$$

(66)

then, the Warburg impedance can also be expressed as:

$$\begin{array}{c}\hfill Z_{\mathrm{w}}\left(j\omega \right)={\displaystyle \frac{1}{\frac{1}{Z\left(\omega \right)-j\omega L-R_{\mathsf{\Omega }}-\frac{R_{\mathrm{SEI}}}{1+j\omega R_{\mathrm{SEI}}{C}_{\mathrm{SEI}}}}-j\omega C_{\mathrm{DL}}}}-R_{\mathrm{CT}}\end{array}$$

(67)

Take the real part on both sides of the above equation, at $\omega ={\omega }_k$$(k_{DF1}\le k\le k_{DF2})$,

$$\begin{array}{c}\hfill \mathrm{Re}\left(Z_{\mathrm{w}}\left(j{\omega }_k\right)\right)=\mathrm{Re}\left({\displaystyle \frac{1}{\frac{1}{Z\left({\omega }_k\right)-j{\omega }_{k}L-R_{\mathsf{\Omega }}-\frac{R_{\mathrm{SEI}}}{1+j{\omega }_k{R}_{\mathrm{SEI}}{C}_{\mathrm{SEI}}}}-j{\omega }_k{C}_{\mathrm{DL}}}}-R_{\mathrm{CT}}\right)\end{array}$$

(68)

The real part of the Warburg impedance can be noted as:

$$\begin{array}{cc}\hfill & W_r\left(k\right)\triangleq \mathrm{Re}\left(Z_{\mathrm{w}}\left(j{\omega }_k\right)\right)\mathrm{s}.\mathrm{t}.k_{DF1}\le k\le k_{DF2}\hfill \end{array}$$

(69)

Take the real part on both sides of (12), we get

$$\begin{array}{c}\hfill W_r\left(k\right)={\displaystyle \frac{\sigma }{\sqrt{{\omega }_k}}}\end{array}$$

(70)

Thus,

$$\begin{array}{c}\hfill {\tilde{\sigma }}_k^{\mathrm{ES}}=W_r\left(k\right)\sqrt{{\omega }_k}\end{array}$$

(71)

Finally, average all the estimates ${\tilde{\sigma }}_k^{\mathrm{ES}}$ to obtain the final estimate

$$\begin{array}{c}\hfill {\widehat{\sigma }}_{\mathrm{ES}}={\displaystyle \frac{1}{k_{DF2}-k_{DF1}+1}}\sum _{k=k_{DF1}}^{k_{DF2}}{\tilde{\sigma }}_k^{\mathrm{ES}}\end{array}$$

(72)

4.2. Specify the Range of Parameters for Exhaustive Search

As presented in Section 3, the rough estimations of ECM parameters are given; therefore, we can assign the lower and upper bound for each parameter such that the exhaustive search approach can identify the most suitable parameter between the boundary. As shown in Algorithm 5, the range of each ECM parameter is assigned based on the empirical coefficient. In this paper, the size of possible values in each ECM parameter is restricted to 20 such that the computational time is within the acceptable range. Algorithm 5 uses the MATLAB command ’linspace’.

Algorithm 5 Set the range of ECM parameters estimated by LS approach.

4.3. Implement Exhaustive Search

Algorithm 6 describes the ES approach that can be applied to precisely identify the ECM parameters; where, the input ${\widehat{R}}_{\mathsf{\Omega }}$ is estimated via (3), and $\widehat{L}$ is estimated via (4). Figure 9 shows the MAE evaluated from the initial iteration throughout the end of the exhaustive search process; by finding the lowest MAE in this figure, one can identify the best estimation of ECM parameters via the ES approach.

Algorithm 6 Exhaustive Search Approach.

5. Nonlinear Least Squares Approach

The concept of implementing the NLS approach based on the Monte-Carlo run is shown in Figure 10. This approach randomly selects initial guesses of the ECM parameters in each Monte Carlo run to fit the EIS data in different cases.

5.1. Objective Function

The goal is to find the optimized ECM parameters to minimize the following function.

$$\begin{array}{c}\hfill \widehat{\mathbf{x}}=arg\underset{\mathbf{x}}{min}\sum _{k=1}^N|\widehat{Z}\left({\omega }_k\right)-Z\left({\omega }_k\right)|\end{array}$$

(73)

Where,

$$\begin{array}{cc}\hfill \mathbf{x}& =[R_{\mathsf{\Omega }},L,R_{SEI},C_{SEI},R_{CT},C_{DL},\sigma ,m]\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill \widehat{Z}\left({\omega }_k\right)& =\widehat{\mathbf{x}}\left(1\right)+j{\omega }_k\widehat{\mathbf{x}}\left(2\right)+{\displaystyle \frac{1}{\frac{1}{\widehat{\mathbf{x}}\left(3\right)}+j{\omega }_k\widehat{\mathbf{x}}\left(4\right)}}\hfill \\ \hfill & +{\displaystyle \frac{1}{\frac{1}{\widehat{\mathbf{x}}\left(5\right)+(1-j\widehat{\mathbf{x}}\left(8\right))\frac{\widehat{\mathbf{x}}\left(7\right)}{\sqrt{{\omega }_k}}}+j\omega \widehat{\mathbf{x}}\left(6\right)}}\hfill \end{array}$$

(75)

and $Z\left({\omega }_k\right)$ is the measured impedance at ${\omega }_k$, N is total number of measurements, $|·|$ denotes the absolute value of the complex number.

The nonlinear least square approach is implemented in MATLAB using the built-in function $lsqnonlin$ to fit the EIS data.

5.2. Initial Guess

Instead of setting the lower bound and upper bound for the NLS approach, in this paper, we randomly select initial guess in each Monte-Carlo runs to try different NLS fitting based on different initial guess, this way, to find the best fit among all cases. In MATLAB the initial guess is defined as:

$$\begin{array}{c}\hfill {\mathbf{x}}_{\mathbf{0}}=[\mathbf{abs}\left(\mathbf{randn}(\mathbf{1},\mathbf{7})\right),\mathbf{1}]\end{array}$$

(76)

5.3. Algorithm Switch

To reach the best fit, the NLS approach will apply ’trust-region-reflective’ algorithm [24,25] and ’levenberg-margquardt’ algorithm [26,27,28] to fit the EIS data with different initial guess, after that, the algorithm that can achieve the lowest mean absolute error (MAE) will be select as the approach that can reach the best fit. The detailed approach is shown in Algorithm 7. This algorithm uses the following MATLAB commands: ’abs’ and ’randn’.

Algorithm 7 NLS Approach.

6. Implementation

This paper implements three ECM parameter estimation approaches in MATLAB R2020a with a 3 GHz Processor and 16 GB RAM.

6.1. Simulate EIS Data

The simulated EIS data was generated using Algorithm 8, where the frequency ranges from 0.01 Hz to 10 kHz, and the number of EIS measurements is 121; this is to keep the conformity with the sampling size we set in the EIS experiment [20]. In addition, Table 1 shows the true ECM parameters for simulating EIS data. To generate various noise levels for the simulated EIS measurements, we implement Gaussian noise with zero-mean, standard deviation ${\sigma }_{noise}$, and independent outcome.

Algorithm 8 uses the MATLAB command ’linspace’, ’real’, ’imag’, and ’randn’.

Algorithm 8 EIS simulation.

6.2. Collect Real EIS Data

The impedance data is measured from two Li-ion batteries: LG 18650 and Molicel 21700. In addition, the specifications of LG and Molicel batteries are shown in [20] Table 1. The data are collected using the Arbin battery cycler (Model: LBT21084UC), which has 16 channels that can operate in parallel. In this experiment, eight channels were used to collect data simultaneously at room temperature (23 °C).

The EIS data are measured by the EIS device (Gamry interface 5000P). We operated the Gamry EIS device and Arbin battery cycler using the software $MitsPro$ provided by Arbin company. The voltage measurement error of the Gamry EIS device is 0.2 mV as specified from [29]. In this paper, we validated LS, ES, and NLS approaches on EIS data collected from one LG and one Molicel battery when the SOC is at 90%, 50% and 10%, while discharging from a fully charged state.

7. Results

In this section, fitting results obtained from the simulated and real EIS data are shown and discussed.

7.1. Estimation Results of ECM Parameters Using Simulated EIS Data

The comparisons of simulated EIS data fitted by LS, ES, and NLS approaches are shown in Figure 11a–d. It can be observed that the LS approach shows insufficient goodness of fitting, whereas ES and NLS approaches generally reach considerable fitting accuracy. Table 2 shows that at any noise level, the computational time of the LS approach is the fastest; however, the fitting accuracy MAE is the lowest. On the contrary, the ES approach has the slowest computational time, but the fitting accuracy is significantly improved compared with the LS approach. The NLS approach reaches the lowest MAE and is considerably faster than the ES approach. Furthermore, with the noise level decreasing, the MAE decreases.

As shown in Figure 12a–d, it is clear that at any noise level, the percentage error of the estimated RC components reaches the highest when using the LS approach and reaches the lowest when using NLS approach, except for that the percentage error of $C_{SEI}$ estimated via NLS is higher than that estimated via ES when ${\sigma }_{noise}$ = 0.1075 $\mathbf{m}\mathsf{\Omega }$; besides, at any noise level, the percentage errors of estimated ECM parameters using NLS approach are well below $5\%$, this shows significantly higher estimation accuracy compared with ES and LS approaches.

7.2. Estimation Results of ECM Parameters Using Real EIS Data

LG 18650 and Molicel 21700 Li-ion batteries were selected to validate whether the LS, ES, and NLS approaches show consistency in fitting real EIS data that are collected at $90\%$, $50\%$, and $10\%$ SOC.

Figure 13a–c show the fitted EIS of LG 18650 battery using LS, ES, and NLS approaches; Figure 13d–f show the fitted EIS of Molicel battery using same approaches. Compare to the LS approach, both ES and NLS approaches show higher fitting accuracy.

In Table 3, it can be observed that when fitting the LG battery’s EIS data at any SOC level, the EIS approach outperforms the ES and LS approaches; furthermore, in terms of the computational time, the ES approach is considerably slower than the NLS approach. Though the LS approach is the fastest, the MAE is the highest among all SOC levels. Additionally, in Table 4, the validation on the Molicel battery shows consistent results.

8. Conclusions and Discussions

This paper presented the LS, ES, and NLS approaches to extract ECM parameters through battery impedance measurements. Compared to the LS approach, the ES and NLS approach can extract ECM parameters more accurately. Though the LS approach shows insufficient goodness of fitting at various noise levels, it can boost the fitting accuracy of the ES approach by offering initial estimations; however, it is worth mentioning that while the ECM contains more than two RC components, the computation time will increase significantly such that this approach will be infeasible.

When fitting the simulated EIS data, both ES and NLS approaches show considerably high accuracy at each noise level, and the fitting accuracy increases as the noise decreases. When fitting the battery EIS measurements, the NLS approach still shows faster and more accurate fitting performance than the ES approach; this result is validated in simulated EIS data.

In future works, we will investigate deploying the NLS approach to the BMS board combined with the rapid EIS measurement hardware to improve the accuracy and computational time for ECM parameters estimation; the BMS can then adopt these precisely estimated ECM parameters for more accurate online SOC/SOH estimation.