Robust Approach to Battery Equivalent-Circuit-Model Parameter Extraction Using Electrochemical Impedance Spectroscopy

Abstract

Electrochemical impedance spectroscopy (EIS) is a well-established method of battery analysis, where the response of a battery to either a voltage or current excitation signal spanning a wide frequency spectrum is measured and analyzed. State-of-the-art EIS analysis is limited to high-precision measurement systems within laboratory environments. In order to be relevant in practical applications, EIS analysis needs to be carried out with low-cost sensors, which suffer from high levels of measurement noise. This article presents an approach to estimate the equivalent circuit model (ECM) parameters of a Li-Ion battery pack based on EIS measurements in the presence of high levels of noise. The proposed algorithm consists of a fast Fourier transform, feature extraction, curve fitting, and least-squares estimation. The results of the proposed parameter-estimation algorithm are compared to that of recent work for objective performance comparison. The error analysis of the proposed approach, in comparison to the existing approach, demonstrated significant improvement in parameter estimation accuracy in low signal-to-noise ratio (SNR) regions. Results show that the proposed algorithm significantly outperforms the previous method under high-measurement-noise scenarios without requiring a significant increase in computational resources.

1. Introduction

Li-ion batteries (LIBs) have seen widespread adoption due to their higher energy and power densities along with longer lifespan [1]. They are widely used in electrified vehicles, consumer electronics, and stationary energy-storage systems. However, their performance is highly dependent on operating conditions (i.e., temperature, load current, etc.). Meanwhile, one of the key challenges in an energy-storage system [2] is to understand the system’s availability from the end-user perspective; this requires knowledge of different battery parameters to monitor, control, and forecast the system’s behaviour [3,4]. Batteries are managed with battery management systems (BMS), which are designed to enhance their performance, prolong their lifespan, and ensure their safety. The performance of a BMS depends largely on its ability to estimate equivalent-circuit-model (ECM) parameters, which characterize the battery behaviour. Estimated battery characterization parameters are useful in determining the state of health (SOH) and remaining useful life (RUL), which characterize battery performance [5,6]. For battery fuel gauging, parameters such as open-circuit voltage (OCV), state of charge (SOC), depth of discharge (DOD), and capacity need to be estimated [7]. The BMS uses the measured voltage, current and temperature information to keep the battery within safe limits by the online estimation of ECM parameters [8].

The state and parameter estimation of LIB is based on a battery model. In practice, ECMs are the most commonly used battery models because of their simple structure and the availability of computationally efficient approaches to estimate their parameters [9]. In these models, electric circuit elements, such as resistance, capacitance and voltage source, are used to represent the electrical characteristics of LIBs [10,11]. Widely used ECMs include the R-int model, RC model, Thevenin model, enhanced Thevenin model, PNGV model and a plethora of combined models featuring RC circuits of different orders [12,13].

The ECM parameter-estimation approaches fall into two categories; time-domain and frequency-domain approaches. Time-domain approaches [7,14] use the available voltage and current measurements during battery usage to extract the ECM parameters. Consequently, time-domain approaches are suitable for real-time parameter estimation [15]. On the other hand, frequency-domain approaches require special excitation signals spanning a wide range of frequencies. These signals can be superimposed on existing currents during charging as well as discharging. Electrochemical impedance spectroscopy (EIS), introduced by Heaviside in 1894 [16], is a widely studied frequency-domain approach to characterize the electrical and electrochemical behaviour of batteries [17]. The measured voltage and current signals are converted to the frequency domain through the fast Fourier transform (FFT) technique to compute the frequency response, generally known as a Nyquist plot [18,19] of the battery impedance, where the real and imaginary parts of the impedance are plotted on the x-axis and y-axis, respectively. Fitting the Nyquist plot from EIS to a LIB ECM allows the battery kinetic processes to be modelled [20]. In the literature [5,21], the adaptive Randles ECM (AR-ECM) is the most commonly used circuit model for the derivation and interpretation of impedance spectrum from EIS [22].

Significant work has been performed in the literature utilizing the EIS approach for battery state and parameter estimation. Pattipati et al. [23] used a frequency-selective nonlinear least squares estimation technique to estimate the Randles circuit parameters from the EIS data. They also used support vector machines (SVMs) to characterize the SOH and estimate the SOC of the battery. The algorithm proposed in [24] uses the EIS technique to identify and quantify the battery aging process over time. The method is applied to four parallelized Lithium-ion cells for 500 cycles. The parameter estimation strategy proposed in [18] used a frequency sweep method for input excitation. It employed precise impedance measurement techniques, consisting of offset clipping, amplifying circuits, FFT, peak amplitude computation, and phase difference estimation.

In [25], the authors presented an online approach to estimate the uncertainty in parameter estimation using the EIS measurements. A complex nonlinear least squares method was applied to the EIS measurements in [26], to determine the ECM parameters that best fit the measurements. In [27], the authors proposed a distribution-of-relaxation-time-based method to merge the time- and frequency-domain analysis. It was described that the dependency of the estimation algorithm on boundary conditions and initial starting values was alleviated using this approach. Another work proposed the use of a recursive least square algorithm for estimating the equivalent circuit model parameters [28]. In [29], a simplified model from the low-frequency region was considered for simpler online estimation of battery ECM parameters. Optimization-based algorithms were used in [30] to estimate the parameters of a second-order RC model. However, the algorithms did not consider measurement noise and the estimation was performed for the entire SOC region. Thus, the effect of the variation in the parameters with SOC was not taken into account in this work. A polynomial fitting of parameters on the EIS data was used to determine the SOC effect on ECM parameter estimation in [31]. Islam et al. [22,32] presented a circuit parameter-extraction algorithm for an LIB charging system using online EIS [22] and compared their method with a hysteresis-free ECM parameter algorithm in [32]. None of the existing work considered the effect of measurement noise on the performance of the ECM parameter estimation algorithm. The focus of the present work is to develop ECM parameter-estimation algorithms which are robust against measurement noise.

In this paper, an improved online method for ECM parameter estimation based on the Nyquist plot is presented. The key contributions of this paper are as follows. First, various factors affecting the accuracy of frequency-domain ECM parameter estimation are discussed. Specifically, it is demonstrated that the ECM parameter-estimation accuracy deteriorates with increased measurement noise. Then, novel approaches are developed to estimate each parameter representing the AR-ECM of a battery. A systematic performance comparison is performed, in terms of the percent of error versus SNR, for each parameter to be estimated. The proposed approach showed significant performance improvement over the existing method.

The remainder of this paper is organized as follows: Section 2 describes the EIS parameter approach proposed in [22]. Section 3 analyzes how approximation and measurement noise affect the EIS parameter-estimation accuracy. Section 4 presents an improved method for parameter estimation in the presence of measurement noise. Section 5 shows how feature points are extracted. The results from the proposed approach and the existing approach in [22] are compared in Section 6. The paper concludes in Section 7, and future research directions are explained in Section 8.

2. Impedance Response of an LIB

Randles circuit models are one of the widely used ECM in battery analysis [33]. Figure 1a shows a diagram of the AR-ECM. This model consists of the following elements [22]:

• Voltage source, $E_{\mathrm{cell}}$

• Stray inductance, L

• Ohmic resistance, $R_{\mathsf{\Omega }}$

• Solid electrolyte interface (SEI) resistance, $R_{\mathrm{SEI}}$

• SEI capacitance, $C_{\mathrm{SEI}}$

• Charge-transfer(CT) resistance, $R_{\mathrm{CT}}$

• Double-layer (DL) capacitance, $C_{\mathrm{DL}}$

• Warburg impedance, Z_w

where the Warburg impedance [34], defined as

(1)$\begin{array}{c}\hfill Z_{(}j\omega )=(1-j)\frac{\sigma }{\sqrt{\omega }}\end{array}$

describes a phenomenon observed at very low frequencies. It can be noticed that as the frequency decreases, the Warburg impedance increases. When there is a significant frequency in the system, which is the case in practical applications, the Warburg impedance effectively becomes zero. On the other hand, at low frequencies, only the effects of resistors and the Warburg element remain. The impedance spectrum also marks the following frequencies at critical points: ${\omega }_{k0}$, ${\omega }_{k1}$, ${\omega }_{k2}$, ${\omega }_{k3}$, and ${\omega }_{k4}$. The lowest and highest frequencies of the spectrum are ${\omega }_{k0}$ and ${\omega }_{k1}$, respectively. The frequency ${\omega }_{k1}$ denotes the point where the diffusion arc meets the CT arc; ${\omega }_{k2}$ denotes the point where CT Arc meets the SEI arc; and ${\omega }_{k3}$ denotes the point where RL arc meets the SEI arc.

In EIS, an AC perturbation signal (current or voltage) is applied to a battery and its response (voltage or current) is recorded; by jointly analyzing the applied signal and its response, the parameters of the battery ECM can be estimated. Figure 1b shows the typical response of a battery in the frequency domain—this impedance plot is generally known as the Nyquist plot. In order to compute the Nyquist plot, the voltage and current signals, measured in the time domain, are converted to the frequency domain through the discrete Fourier transform (DFT) technique. Let us assume that $z_v\left(t\right)$ and $z_c\left(t\right)$ are the measured voltage and current, respectively, from the battery over a certain time window L, i.e.,

(2)$\begin{array}{c}\hfill \begin{array}{cc}\hfill z_v\left(t\right)& =v\left(t\right)+n_v\left(t\right),t=1,2,\dots ,L\hfill \\ \hfill z_c\left(t\right)& =i\left(t\right)+n_c\left(t\right),t=1,2,\dots ,L\hfill \end{array}\end{array}$

where t indicates time, $v\left(t\right)$ is the true voltage, and $i\left(t\right)$ is the true current. The voltage and current measurement noise $n_v\left(t\right)$ and $n_c\left(t\right)$ are assumed to be zero-mean i.i.d. with standard deviation ${\sigma }_v$ and ${\sigma }_c$, respectively. The Fourier transform of the voltage and current measurements in (2) are defined as

(3)$\begin{array}{c}\hfill \begin{array}{cc}{Z}_c\left(f_{\omega }\right)\hfill & =\mathrm{FFT}\left(z_c\left(k\right)\right)={\displaystyle \sum _{k=1}^L}{z}_c\left(k\right)e^{\frac{-i2\pi k{f}_{\omega }}{L}}\hfill \\ \hfill & =I\left(f_{\omega }\right)+N_c\left(f_{\omega }\right)\hfill \end{array}\end{array}$

(4)$\begin{array}{c}\hfill \begin{array}{cc}{Z}_v\left(f_{\omega }\right)\hfill & =\mathrm{FFT}\left(z_v\left(k\right)\right)={\displaystyle \sum _{k=1}^L}{z}_v\left(k\right)e^{\frac{-i2\pi k{f}_{\omega }}{L}}\hfill \\ \hfill & =V\left(f_{\omega }\right)+N_v\left(f_{\omega }\right)\hfill \end{array}\end{array}$

where $Z_v\left(f_{\omega }\right)$ and $Z_c\left(f_{\omega }\right)$ are voltage and current measurement in the frequency domain at the frequency $f_{\omega }$. Here, $V\left(f_{\omega }\right)$ and $I\left(f_{\omega }\right)$ indicate the Fourier transforms of the noiseless voltage and current, respectively. The Fourier transforms of the voltage and current measurement noises are given by $N_v\left(f_{\omega }\right)$ and $N_c\left(f_{\omega }\right)$, respectively. The impedance at frequency $f_{\omega }$ is now written as

$\begin{array}{cc}\hfill Z\left(f_{\omega }\right)=\frac{Z_v\left(f_{\omega }\right)}{Z_c\left(f_{\omega }\right)}& =\frac{V\left(f_{\omega }\right)+N_v\left(f_{\omega }\right)}{I\left(f_{\omega }\right)+N_c\left(f_{\omega }\right)}\hfill \\ & =\left(V\left(f_{\omega }\right)+N_v\left(f_{\omega }\right)\right)\left(\frac{1}{I\left(f_{\omega }\right)+N_c\left(f_{\omega }\right)}\right)\hfill \end{array}$

The term $\left(\frac{1}{I\left(f_{\omega }\right)+N_c\left(f_{\omega }\right)}\right)$ will be expanded using Taylor’s series expansion, i.e., for any $f(x+h)$, the Taylor’s series expansion about ‘h’ is given by

$\begin{array}{c}\hfill f(x+h)=f\left(x\right)+h{f}^{\prime }\left(x\right)+\frac{h^2}{2!}{f}^{″}\left(x\right)+\frac{h^3}{3!}{f}^{\prime \prime \prime }\left(x\right)+....\end{array}$

Using first-order Taylor’s expansion for $\left(\frac{1}{I\left(f_{\omega }\right)+N_c\left(f_{\omega }\right)}\right)$, where, $f(x+h)=\left(\frac{1}{I\left(f_{\omega }\right)+N_c\left(f_{\omega }\right)}\right)$, $f\left(x\right)=\frac{1}{I\left(f_{\omega }\right)}$, and $h=N_c\left(f_{\omega }\right)$, it could be written as

(5)$\begin{array}{c}\hfill \frac{1}{I\left(f_{\omega }\right)+N_c\left(f_{\omega }\right)}\approx \frac{1}{I\left(f_{\omega }\right)}-\frac{N_c\left(f_{\omega }\right)}{I{\left(f_{\omega }\right)}^2}\end{array}$

$Z\left(f_{\omega }\right)$ is now

(6)$\begin{array}{cc}\hfill Z\left(f_{\omega }\right)& =\left(V\left(f_{\omega }\right)+N_v\left(f_{\omega }\right)\right)\left(\frac{1}{I\left(f_{\omega }\right)}-\frac{N_c\left(f_{\omega }\right)}{I{\left(f_{\omega }\right)}^2}\right)\hfill \end{array}$

Expanding the terms in impedance in (6)

$\begin{array}{c}\hfill Z\left(f_{\omega }\right)=\frac{V\left(f_{\omega }\right)}{I\left(f_{\omega }\right)}-\frac{V\left(f_{\omega }\right)N_c\left(f_{\omega }\right)}{I{\left(f_{\omega }\right)}^2}+\frac{N_v\left(f_{\omega }\right)}{I\left(f_{\omega }\right)}-\frac{N_v\left(f_{\omega }\right)N_c\left(f_{\omega }\right)}{I{\left(f_{\omega }\right)}^2}\end{array}$

This is re-written in the following format:

(7)$\begin{array}{c}\hfill Z\left(f_{\omega }\right)=\frac{V\left(f_{\omega }\right)}{I\left(f_{\omega }\right)}+N_z\left(f_{\omega }\right)\end{array}$

where

$\begin{array}{c}\hfill N_z\left(f_{\omega }\right)=-\frac{V\left(f_{\omega }\right)N_c\left(f_{\omega }\right)}{I{\left(f_{\omega }\right)}^2}+\frac{N_v\left(f_{\omega }\right)}{I\left(f_{\omega }\right)}-\frac{N_v\left(f_{\omega }\right)N_c\left(f_{\omega }\right)}{I{\left(f_{\omega }\right)}^2}\end{array}$

where it can be shown that the noise $N_z\left(f_{\omega }\right)$ is zero-mean.

Remark 1.

The voltage and current measurements in (2) are assumed to represent the AC components. It is also assumed that the DC component remains the same during the FFT. That is, the DC component is assumed to be the same for$t=1,\dots ,L$. Under this assumption, the effect of the DC component can be simply removed by subtracting the mean before performing the FFT. In realistic application, this assumption may not hold, i.e., the value of the DC component may not remain the same for all values for$t=1,\dots ,L$. The result of this violation becomes an estimation error. This is demonstrated later using experimental studies.

Using the AR-ECM shown in Figure 1, the AC impedance $Z\left(f_{\omega }\right)$ can be written as

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

where the parameters are as indicated in Figure 1a. Likewise, the qualitative impedance plot shown in Figure 1b has four branches associated with four specific electrochemical processes [10]. In the first branch, denoted in this paper as the ‘RL Arc’, the effect of the inductive behaviour at high frequencies ($\omega >{\omega }_{k_3}$) as well as the ohmic resistance ($R_{\mathsf{\Omega }}$) can be seen. The second branch (${\omega }_{k_2}<\omega <{\omega }_{k_3}$) consists of a semi-circle, denoted in this paper as the ‘SEI Arc’, related to the SEI. The effect of the DL capacitance and CT resistance at the electrodes represents the second semi-circle (denoted as ‘CT Arc’) in the third branch (${\omega }_{k_1}<\omega <{\omega }_{k_2}$). Finally, the constant slope (denoted as ‘Diffusion Arc’) in the impedance plot in the last branch represents the diffusion processes in the active material of the electrodes; it has a significant effect at very low frequencies (${\omega }_0<\omega <{\omega }_{k_1}$) only.

For fast interpretation of the impedance spectrum, an algorithm is proposed in [32] that obtains the AR-ECM parameters from an impedance spectrum. The remainder of this section summarizes the EIS approach presented in [32].

3. Nonlinear LS Approach for ECM Parameter Estimation

The ECM parameter estimation in the frequency domain can be formally stated as follows: given the frequency response of the system $Z\left(f_{\omega }\right)$ at the frequencies ${\omega }_1,{\omega }_2,\dots ,{\omega }_L$, the ECM parameters that minimize the following cost function are estimated:

(9)$\begin{array}{c}\hfill \widehat{\mathsf{\Theta }}=arg\underset{\mathsf{\Theta }}{min}{∥Z\left({\omega }_i\right)-\widehat{Z}\left({\omega }_i\right)∥}^2,i=1,\dots ,L\end{array}$

where $∥·∥$ denotes the second norm and

(10)$\begin{array}{c}\hfill \mathsf{\Theta }=\{R_{\mathsf{\Omega }},L,R_{\mathrm{SEC}},C_{\mathrm{SEC}},R_{\mathrm{CT}},C_{\mathrm{DL}},Z_W\}\end{array}$

The above optimization can be solved using a non-linear least squares approach. In this paper, the Matlab command ‘lsqnonlin’ was used to perform non-linear least squares optimization.

4. Prior Approach for ECM Parameter Estimation

This section summarizes a previous approach presented in [32] for ECM parameter estimation using EIS measurements.

Consider the impedance in (8) at very high frequencies $\omega >{\omega }_3.$ The capacitive reactance approaches zero at very high frequencies and only the inductive and ohmic resistance remain dominant. Hence, one can write

(11)$\begin{array}{c}\hfill Z\left(\omega \right)\approx j\omega L+R_{\mathsf{\Omega }}\omega >{\omega }_{k_3}\end{array}$

From this, the estimates of L and $R_{\mathsf{\Omega }}$ can be obtained as

(12)$\begin{array}{cc}\hfill {\widehat{R}}_{\mathsf{\Omega }}& =\min \left(z_r\left(k\right)\right)k>k_3\hfill \end{array}$

(13)$\begin{array}{cc}\hfill \widehat{L}& =\frac{\left|\min \left(-z_i\left(k\right)\right)\right|}{\omega }k>k_3\hfill \end{array}$

At very low frequencies, the Warburg impedance (1) becomes dominant. The resistive part of the Warburg impedance is written as

(14)$\begin{array}{c}\hfill R_{=}\frac{\sigma }{\sqrt{\omega }}\end{array}$

Consequently, the real part of $Z\left(\omega \right)$, when $\omega <{\omega }_{k_1}$, can be written as

(15)$\begin{array}{c}\hfill z_r\left(k\right)\approx R_{\mathsf{\Omega }}+R_{\mathrm{SEI}}+R_{\mathrm{CT}}+R_{\mathrm{w}}k<k_1\end{array}$

From (15), $\sigma$ can be calculated as the slope of $z_r\left(k\right)$ versus $\frac{1}{{\sqrt{\omega }}_k}$ such that $k<k_1$. To compute $\sigma$, two low frequencies are chosen and the corresponding resistance value is taken from the impedance plot. Let us select these two frequencies as follows

(16)$\begin{array}{cc}\hfill {\omega }_a& ={\omega }_0\hfill \end{array}$

(17)$\begin{array}{cc}\hfill {\omega }_b& =\omega \mathrm{s}.\mathrm{t}.{\omega }_0<\omega <{\omega }_{k_1}\hfill \end{array}$

Then, the Warburg coefficient can be written as follows:

(18)$\begin{array}{c}\hfill \widehat{\sigma }=\frac{\left(\sqrt{{\omega }_a{\omega }_b}\right)(z_r\left(a\right)-z_r\left(b\right))}{\sqrt{{\omega }_a}-\sqrt{{\omega }_b}}\end{array}$

Now, let us consider the two arcs (SEI Arc and CT Arc) in the Nyquist plot to determine the value of $R_{\mathrm{SEI}}$, $C_{\mathrm{SEI}}$, $R_{}$ and $C_{}$. First, consider the CT Arc which occurs in lower frequencies, i.e., ${\omega }_{k_1}<\omega <{\omega }_{k_2}$. The (Faradaic) impedance due to $R_{\mathrm{CT}}$ and $C_{\mathrm{DL}}$, in this region, is

(19)$\begin{array}{c}\hfill \begin{array}{cc}{Z}_{\mathrm{F}}\left(\omega \right)\hfill & =\frac{1}{\frac{1}{R_{\mathrm{CT}}+Z_{(}j\omega )}+j\omega C_{\mathrm{DL}}}\hfill \\ \hfill & =\frac{\frac{1}{R_{\mathrm{CT}}+Z_{(}j\omega )}-j\omega C_{\mathrm{DL}}}{{\left(\frac{1}{R_{\mathrm{CT}}+Z_{(}j\omega )}\right)}^2+{\omega }^2{C^2}_{\mathrm{DL}}}{\omega }_{k_1}<\omega <{\omega }_{k_2}\hfill \end{array}\end{array}$

By ignoring the effect of Warburg impedance, the impedance corresponding to the CT Arc can be written as

(20)$\begin{array}{c}\hfill Z_{\mathrm{CT}}\left(\omega \right)\approx \frac{\frac{1}{R_{\mathrm{CT}}}-j\omega C_{\mathrm{DL}}}{{\left(\frac{1}{R_{\mathrm{CT}}}\right)}^2+{\omega }^2{C^2}_{\mathrm{DL}}}{\omega }_{k_1}<\omega <{\omega }_{k_2}\end{array}$

At the peak of CT Arc, one can observe (see Figure 1b)

(21)$\begin{array}{cc}\hfill \frac{1}{\omega C_{\mathrm{DL}}}& =R_{\mathrm{CT}}\mathrm{at}\omega ={\omega }_{\mathrm{CT},\mathrm{peak}}\hfill \end{array}$

(22)$\begin{array}{cc}\hfill \left|\mathrm{Im}\left(Z_{\mathrm{CT}}\left(\omega \right)\right)\right|& =\frac{1}{2}{R}_{\mathrm{CT}}\mathrm{at}\omega ={\omega }_{\mathrm{CT},\mathrm{peak}}\hfill \end{array}$

and the following two estimates can be obtained:

(23)$\begin{array}{cc}\hfill {\widehat{R}}_{\mathrm{CT}}& =2\left|\mathrm{Im}\left(Z_{\mathrm{CT}}\left(\omega \right)\right)\right|\mathrm{at}\omega =2\pi f_{\mathrm{CT},\mathrm{peak}}\hfill \end{array}$

(24)$\begin{array}{cc}\hfill {\widehat{C}}_{\mathrm{DL}}& =\frac{1}{\omega {\widehat{R}}_{\mathrm{CT}}}\mathrm{at}\omega =2\pi f_{\mathrm{CT},\mathrm{peak}}\hfill \end{array}$

where $f_{\mathrm{CT},\mathrm{peak}}$ denotes the frequency corresponding to the peak of the CT Arc. Let us now consider the SEI Arc in the range of ${\omega }_{k_2}<\omega <{\omega }_{k_3}$. In this region, the impedance is given by

(25)$\begin{array}{c}\hfill \begin{array}{cc}{Z}_{\mathrm{SEI}}\left(\omega \right)\hfill & =\frac{1}{\frac{1}{R_{\mathrm{SEI}}}+j\omega C_{\mathrm{SEI}}}\hfill \\ \hfill & =\frac{\frac{1}{R_{\mathrm{SEI}}}-j\omega C_{\mathrm{SEI}}}{\frac{1}{{R^2}_{\mathrm{SEI}}}+{\omega }^2{C^2}_{\mathrm{SEI}}}{\omega }_{k_2}<\omega <{\omega }_{k_3}\hfill \end{array}\end{array}$

At the peak of SEI Arc, we have

(26)$\begin{array}{c}\hfill \frac{1}{\omega C_{\mathrm{SEI}}}=R_{\mathrm{SEI}}\mathrm{at}\omega ={\omega }_{\mathrm{SEI},\mathrm{peak}}\end{array}$

Based on the above observation, we have

(27)$\begin{array}{c}\hfill \left|\mathrm{Im}\left(Z_{\mathrm{SEI}}\left(\omega \right)\right)\right|=\frac{1}{2}{R}_{\mathrm{SEI}}\mathrm{at}\omega ={\omega }_{\mathrm{SEI},\mathrm{peak}}\end{array}$

and the following two estimates can be obtained:

(28)$\begin{array}{cc}\hfill {\widehat{R}}_{\mathrm{SEI}}& =2\left|\mathrm{Im}\left(Z_{\mathrm{SEI}}\left(\omega \right)\right)\right|\mathrm{at}\omega =2\pi f_{\mathrm{SEI},\mathrm{peak}}\hfill \end{array}$

(29)$\begin{array}{cc}\hfill {\widehat{C}}_{\mathrm{SEI}}& =\frac{1}{\omega {\widehat{R}}_{\mathrm{SEI}}}\mathrm{at}\omega =2\pi f_{\mathrm{SEI},\mathrm{peak}}\hfill \end{array}$

where $f_{\mathrm{SEI},\mathrm{peak}}$ denotes the frequency corresponding to the peak of the SEI Arc. Despite its simplicity and low processing time, the algorithm reviewed in this section has some drawbacks due to the effect of approximation and due to effect of measurement noise. These effects are described next. Despite its simplicity and low processing time, the algorithm reviewed in this section has some drawbacks, which will be analyzed in the next section.

In this section, a previous approach presented in [32] to estimate the ECM parameters of a battery based on the EIS Nyquist spectrum was reviewed. The reviewed approach exploits the geometric properties of the Nyquit spectrum for parameter estimation. However, the existing approach relies on identifying feature points, such as the peak of the SEI Arc, for parameter estimation. This approach is likely to suffer when there is significant measurement noise. The next section discusses how the performance of the existing approach may degrade in the presence of measurement noise.

5. Performance Analysis

5.1. Effect of Approximation

In Section 2, peak values of SEI Arc and CT Arc are found and the value of RC elements are calculated using (22) to (29). However, it should be noted that, for these calculations, the effect of estimated Warburg impedance on CT Arc is neglected. Figure 2 illustrates how this approximation will affect the parameter estimation. For the case shown in Figure 2, the frequency at which the CT Arc reaches its maximum value is changed by more than 50 percent (from 0.59 to 0.91); that will affect the values of $R_{\mathrm{CT}}$ and $C_{\mathrm{DL}}$ according to (23).

5.2. Effect of Measurement Noise

With increasing levels of measurement noise, the approach summarized in Section 2 becomes less accurate. In order to analyze the performance, let us first define the SNR as

(30)$\begin{array}{cc}\hfill \mathrm{SNR}& =10log\left(\frac{P_{\mathrm{Signal}}}{P_{\mathrm{Noise}}}\right)\hfill \end{array}$

where the unit of SNR is decibels (dB). Figure 3 shows the Nyquist plots at four different SNR values (0, 5, 15, 30 dB). The algorithm summarized in Section 2 is applied to estimate the AR-ECM parameters in each case. The estimated parameters are then used to generate the Nyquist plot. Ideally, both plots should coincide. It can be noticed that with increasing noise, the discrepancies become prominent.

This section discussed the difficulties of applying geometric approaches, reviewed in Section 4, for ECM parameter estimation when there is significant measurement noise. It is shown that noise affects the accuracy of feature-point detection needed for parameter estimation. This section also introduced SNR as a tool to objectively quantify the level of measurement noise in practical applications.

6. Improved Approach for Parameter Estimation

Let us denote the real and imaginary parts of the frequency response at ${\omega }_k$ as

(31)$\begin{array}{c}\hfill \begin{array}{cc}{z}_r\left(k\right)\hfill & =z_r\left({\omega }_k\right)=\mathrm{Re}\left(Z\left({\omega }_k\right)\right)\hfill \\ z_i\left(k\right)\hfill & =z_i\left({\omega }_k\right)=\mathrm{Im}\left(Z\left({\omega }_k\right)\right)\hfill \end{array}\end{array}$

The frequencies at different branches of the Nyquist plot are defined as follows:

• Warburg: ${\omega }_0,{\omega }_1,\dots ,{\omega }_{k_1}$

• CT: ${\omega }_{k_1+1},{\omega }_{k_1+2},\dots ,{\omega }_{k_2}$

• SEI: ${\omega }_{k_2+1},{\omega }_{k_2+2},\dots ,{\omega }_{k_3}$

• RL: ${\omega }_{k_3+1},{\omega }_{k_3+2},\dots ,{\omega }_{k_4}$

That is, in total, there are $k_4$ frequency pairs at which the impedance measurements were computed. The real and imaginary parts of the impedance measurements at each of the above frequencies are denoted as follows:

• Warburg:

  (32)$\begin{array}{c}\hfill [z_r\left(0\right),z_i\left(0\right)],[z_r\left(1\right),z_i\left(1\right)],\dots ,[z_r\left(k_1\right),z_i\left(k_1\right)]\end{array}$

• CT:

  (33)$\begin{array}{c}\hfill [z_r(k_1+1),z_i(k_1+1)],[z_r(k_1+2),z_i(k_1+2)],\dots ,[z_r\left(k_2\right),z_i\left(k_2\right)]\end{array}$

• SEI:

  (34)$\begin{array}{c}\hfill [z_r(k_2+1),z_i(k_2+1)],[z_r(k_2+2),z_i(k_2+2)],\dots ,[z_r\left(k_3\right),z_i\left(k_3\right)]\end{array}$

• RL:

  (35)$\begin{array}{c}\hfill [z_r(k_3+1),z_i(k_3+1)],[z_r(k_3+2),z_i(k_3+2)],\dots ,[z_r\left(k_4\right),z_i\left(k_4\right)]\end{array}$

6.1. Estimation of Warburg Coefficient

Several frequencies in the Diffusion Arc can be written as

(36)$\begin{array}{c}\hfill \begin{array}{cc}\hfill z_r\left(0\right)-z_r\left(k_1\right)& =\sigma \left(\frac{1}{\sqrt{{\omega }_0}}-\frac{1}{\sqrt{{\omega }_{k_1}}}\right)\hfill \\ \hfill z_r\left(1\right)-z_r(k_1-1)& =\sigma \left(\frac{1}{\sqrt{{\omega }_{k_1}}}-\frac{1}{\sqrt{{\omega }_{k_1-1}}}\right)\hfill \\ \hfill & ⋮\hfill \\ \hfill z_r\left(n\right)-z_r(k_1-n)& =\sigma \left(\frac{1}{\sqrt{{\omega }_n}}-\frac{1}{\sqrt{{\omega }_{k_1-n}}}\right)\hfill \end{array}\end{array}$

where $n<k_1/2$. The observations in (36) were selected in such a way that the quantity $z_r\left(i\right)-z_r\left(j\right)$ could be as high as possible—this strategy is designed to reduce the effect of noise in the observations.

The observations (36) can be written in matrix form as

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

where

(38)$\begin{array}{c}\hfill \tilde{\mathbf{z}}=\left[\begin{array}{c}{z}_r\left(0\right)-z_r\left(k_1\right)\\ z_r\left(1\right)-z_r(k_1-1)\\ ⋮\\ z_r\left(n\right)-z_r(k_1-n)\end{array}\right],\mathbf{b}=\left[\begin{array}{c}\left(\frac{1}{\sqrt{{\omega }_0}}-\frac{1}{\sqrt{{\omega }_{k_1}}}\right)\\ \left(\frac{1}{\sqrt{{\omega }_{k_1}}}-\frac{1}{\sqrt{{\omega }_{k_1-1}}}\right)\\ ⋮\\ \left(\frac{1}{\sqrt{{\omega }_n}}-\frac{1}{\sqrt{{\omega }_{k_1-n}}}\right)\end{array}\right]\end{array}$

The least-square estimate of $\sigma$ is

(39)$\begin{array}{c}\hfill \widehat{\widehat{\sigma }}=\frac{\left({\mathbf{b}}^T\tilde{\mathbf{z}}\right)}{\left({\mathbf{b}}^T\mathbf{b}\right)}\end{array}$

6.2. Estimation of CT Components

Let us denote an impedance measurement in the CT Arc as

(40)$\begin{array}{c}\hfill \begin{array}{c}\hfill z_r\triangleq z_r\left(k\right)\mathrm{s}.\mathrm{t}.k_1<k\le k_2\\ \hfill z_i\triangleq z_i\left(k\right)\mathrm{s}.\mathrm{t}.k_1<k\le k_2\end{array}\end{array}$

The measurements in (40) will satisfy the following circular equation

(41)$\begin{array}{c}\hfill \begin{array}{c}{z}_r^2+z_i^2+a{z}_r+b=0\hfill \\ {\left(z_r+\frac{a}{2}\right)}^2+z_i^2=\frac{a^2}{4}-b\hfill \end{array}\end{array}$

where it was assumed that the center of the circle lies on the real axis (see Figure 1). The center of the circle (41) is denoted as $(x_{\mathrm{CT}},0)$ where

(42)$\begin{array}{c}\hfill x_{\mathrm{CT}}=-\frac{a}{2}\end{array}$

and the radius of the circle (41) is

(43)$\begin{array}{c}\hfill r_{\mathrm{CT}}=\sqrt{\frac{a^2}{4}-b}\end{array}$

Based on (41), one can notice that the argument of the square root is always positive. It is now easy to see that the estimate of the resistance $R_{\mathrm{CT}}$ is

(44)$\begin{array}{cc}\hfill {\widehat{\widehat{R}}}_{\mathrm{CT}}& =2\sqrt{\frac{{\widehat{a}}^2}{4}-\widehat{b}}\hfill \end{array}$

where $\widehat{a}$ and $\widehat{b}$ are estimates of a and b, respectively.

In order to estimate a and b, the pairs of impedance measurements shown in (33) can be substituted in (41) to obtain the following sets of equations

(45)$\begin{array}{c}\hfill \begin{array}{cc}\hfill -(z_r{(k_1+1)}^2+z_i{(k_1+1)}^2)& =a{z}_r(k_1+1)+b\hfill \\ \hfill -(z_r{(k_1+2)}^2+z_i{(k_1+2)}^2)& =a{z}_r(k_1+2)+b\hfill \\ \hfill & ⋮\hfill \\ \hfill -(z_r{\left(k_2\right)}^2+z_i{\left(k_2\right)}^2)& =a{z}_r\left(k_2\right)+b\hfill \end{array}\end{array}$

and the above can be written in matrix form as

(46)$\begin{array}{c}\hfill \mathbf{z}=\mathbf{B}{\mathbf{x}}_{\mathrm{CT}}\end{array}$

where

(47)$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathbf{z}=\left[\begin{array}{c}-(z_r{(k_1+1)}^2+z_i{(k_1+1)}^2)\\ -(z_r{(k_1+2)}^2+z_i{(k_1+2)}^2)\\ ⋮\\ -(z_r{\left(k_2\right)}^2+z_i{\left(k_2\right)}^2)\end{array}\right],\mathbf{B}=\left[\begin{array}{cc}{z}_r(k_1+1)& 1\\ z_r(k_1+2)& 1\\ ⋮\\ z_r\left(k_2\right)& 1\end{array}\right],{\mathbf{x}}_{\mathrm{CT}}=\left[\begin{array}{c}a\\ b\end{array}\right]\end{array}\end{array}$

The least-square estimate of ${\mathbf{x}}_{\mathrm{CT}}$ is

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

and the estimates of a and b are

(49)$\begin{array}{cc}\hfill \widehat{a}& ={\widehat{\mathbf{x}}}_{\mathrm{CT}}\left(1\right)\hfill \end{array}$

(50)$\begin{array}{cc}\hfill \widehat{b}& ={\widehat{\mathbf{x}}}_{\mathrm{CT}}\left(2\right)\hfill \end{array}$

which will be substituted in (44) to estimate $R_{\mathrm{CT}}$. Then, for estimating $C_{\mathrm{DL}}$, the inverse of (19) is calculated as:

(51)$\begin{array}{c}\hfill \begin{array}{cc}\frac{1}{Z_{\mathrm{F}}\left(\omega \right)}\hfill & =\frac{1}{R_{\mathrm{CT}}+Z_{(}j\omega )}+j\omega C_{\mathrm{DL}}\hfill \\ \hfill & =\frac{1}{R_{\mathrm{CT}}+\frac{\sigma }{\sqrt{\omega }}-j\frac{\sigma }{\sqrt{\omega }}}+j\omega C_{\mathrm{DL}}\hfill \\ \hfill & =\frac{R_{\mathrm{CT}}+\frac{\sigma }{\sqrt{\omega }}+j\frac{\sigma }{\sqrt{\omega }}}{{\left(R_{\mathrm{CT}}+\frac{\sigma }{\sqrt{\omega }}\right)}^2+\frac{{\sigma }^2}{\omega }}+j\omega C_{\mathrm{DL}}\hfill \end{array}\end{array}$

Hence,

(52)$\begin{array}{cc}\hfill \mathrm{Im}\left(\frac{1}{Z_{\mathrm{F}}\left(\omega \right)}\right)& =\frac{\frac{\sigma }{\sqrt{\omega }}}{{\left(R_{\mathrm{CT}}+\frac{\sigma }{\sqrt{\omega }}\right)}^2+\frac{{\sigma }^2}{\omega }}+\omega C_{\mathrm{DL}}\hfill \end{array}$

and by substituting the estimate of $R_{\mathrm{CT}}$ we obtain (for $\omega ={\omega }_k$)

(53)$\begin{array}{c}\hfill \begin{array}{c}\hfill {\overline{C}}_{\mathrm{DL}}\left(k\right)=\frac{1}{{\omega }_k}\left(\mathrm{Im}\left(\frac{1}{Z_{\mathrm{F}}\left({\omega }_k\right)}\right)-\frac{\frac{\sigma }{\sqrt{{\omega }_k}}}{{\left({\widehat{\widehat{R}}}_{\mathrm{CT}}+\frac{\sigma }{\sqrt{{\omega }_k}}\right)}^2+\frac{{\sigma }^2}{{\omega }_k}}\right)\\ \hfill k=k_1+1,k_1+2,\dots ,k_2\end{array}\end{array}$

Finally, all estimates of ${\overline{C}}_{\mathrm{DL}}\left(k\right)$ are averaged to obtain

(54)$\begin{array}{c}\hfill {\widehat{\widehat{C}}}_{\mathrm{DL}}=\frac{1}{k_2-k_1}\sum _{k=k_1+1}^{k_2}{\overline{C}}_{\mathrm{DL}}\left(k\right)\end{array}$

6.3. Estimation of SEI Components

Let us denote an impedance measurement in the SEI Arc as

(55)$\begin{array}{c}\hfill \begin{array}{c}\hfill y_r\triangleq z_r\left(k\right)\mathrm{s}.\mathrm{t}.k_2<k\le k_3\\ \hfill y_i\triangleq z_i\left(k\right)\mathrm{s}.\mathrm{t}.k_2<k\le k_3\end{array}\end{array}$

The measurements in (55) will satisfy the following circular equation

(56)$\begin{array}{c}\hfill y_r^2+y_i^2+c{y}_r+d=0\end{array}$

where it is assumed that the center of the circle lies on the real axis (see Figure 1). The center of the circle (41) can be denoted as $(x_{\mathrm{SEI}},0)$, where

(57)$\begin{array}{c}\hfill x_{\mathrm{SEI}}=-\frac{c}{2}\end{array}$

and the radius of the circle (41) is

(58)$\begin{array}{c}\hfill r_{\mathrm{SEI}}=\sqrt{\frac{c^2}{4}-d}\end{array}$

By the same reasoning in Section 6.2, the argument of the square root can be shown to be always positive. It now is easy to see that the estimate of the resistance $R_{\mathrm{CT}}$ is

(59)$\begin{array}{cc}\hfill {\widehat{\widehat{R}}}_{\mathrm{SEI}}& =2\sqrt{\frac{{\widehat{c}}^2}{4}-\widehat{d}}\hfill \end{array}$

where $\widehat{c}$ and $\widehat{d}$ are estimates of c and d, respectively.

In order to estimate c and d, the pairs of impedance measurements shown in (34) can be substituted in (56) to obtain the following sets of equations

(60)$\begin{array}{c}\hfill \begin{array}{cc}\hfill -(y_r{(k_2+1)}^2+y_i{(k_2+1)}^2)& =c{y}_r(k_2+1)+d\hfill \\ \hfill -(y_r{(k_2+2)}^2+y_i{(k_2+2)}^2)& =c{y}_r(k_2+2)+d\hfill \\ \hfill & ⋮\hfill \\ \hfill -(y_r{\left(k_3\right)}^2+y_i{\left(k_3\right)}^2)& =c{y}_r\left(k_3\right)+d\hfill \end{array}\end{array}$

The above can be written in matrix form as

(61)$\begin{array}{c}\hfill \mathbf{y}=\mathbf{A}{\mathbf{x}}_{\mathrm{SEI}}\end{array}$

where

(62)$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathbf{y}=\left[\begin{array}{c}-(y_r{(k_2+1)}^2+y_i{(k_2+1)}^2)\\ -(y_r{(k_2+2)}^2+y_i{(k_2+2)}^2)\\ ⋮\\ -(y_r{\left(k_3\right)}^2+y_i{\left(k_3\right)}^2)\end{array}\right],\\ \hfill \mathbf{A}=\left[\begin{array}{cc}{y}_r(k_2+1)& 1\\ y_r(k_2+2)& 1\\ ⋮\\ y_r\left(k_3\right)& 1\end{array}\right],{\mathbf{x}}_{\mathrm{SEI}}=\left[\begin{array}{c}c\\ d\end{array}\right]\end{array}\end{array}$

The least-square estimate of ${\mathbf{x}}_{\mathrm{SEI}}$ is

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

and the estimates of c and d are

(64)$\begin{array}{cc}\hfill \widehat{c}& ={\widehat{\mathbf{x}}}_{\mathrm{SEI}}\left(1\right)\hfill \end{array}$

(65)$\begin{array}{cc}\hfill \widehat{d}& ={\widehat{\mathbf{x}}}_{\mathrm{SEI}}\left(2\right)\hfill \end{array}$

which will be substituted in (59) to estimate $R_{\mathrm{SEI}}$. Then, $C_{\mathrm{SEI}}$ is simply driven from the inverse of (25):

(66)$\begin{array}{cc}\hfill \frac{1}{Z_{\mathrm{SEI}}\left(\omega \right)}& =\frac{1}{R_{\mathrm{SEI}}}+j\omega C_{\mathrm{SEI}}\hfill \end{array}$

From the above, we obtain (for $\omega ={\omega }_k$)

(67)$\begin{array}{c}\hfill {\overline{C}}_{\mathrm{SEI}}\left({\omega }_k\right)=\left(\frac{1}{{\omega }_k}\right)\mathrm{Im}\left(\frac{1}{Z_{\mathrm{SEI}}\left({\omega }_k\right)}\right)\end{array}$

Finally, all possible estimates of ${\overline{C}}_{\mathrm{SEI}}$ are averaged to obtain

(68)$\begin{array}{c}\hfill {\widehat{\widehat{C}}}_{\mathrm{SEI}}=\frac{1}{k_3-k_2}\sum _{k=k_2+1}^{k_3}{\overline{C}}_{\mathrm{SEI}}\left(k\right)\end{array}$

6.4. Estimation of Resistance and Inductance

Finally, for a better estimation of the ohmic impedance and the inductance, the values are averaged over a range of high frequencies (${\omega }_{k_3}<\omega \le {\omega }_{k_4}$).

(69)$\begin{array}{c}\hfill {\widehat{\widehat{R}}}_{\mathsf{\Omega }}=\frac{1}{k_4-k_3}\sum _{k=k_3+1}^{k_4}{z}_r\left(k\right)\end{array}$

(70)$\begin{array}{c}\hfill \widehat{\widehat{L}}=\frac{1}{(k_4-k_3)}\sum _{k=k_3+1}^{k_4}\frac{z_i\left(k\right)}{{\omega }_k}\end{array}$

in which $z_r\left(k\right)=\mathrm{Re}\left(Z\left({\omega }_k\right)\right)$ and $z_i\left(k\right)=\mathrm{Im}\left(Z\left({\omega }_k\right)\right)$.

7. Feature-Point Extraction

This section presents robust approaches to estimate the feature points $k_1,k_2,$ and $k_3$ (see Figure 1).

7.1. Estimation of $k_1$

Starting from $k_1$, the Nyquist plot changes from Diffusion Arc to CT Arc. In order to estimate $k_1$, the fact that the Diffusion Arc is (approximately) linear is considered. Thereby, the impedance measurements satisfy a linear model

(71)$\begin{array}{c}\hfill z_i\left(k\right)=a_0{z}_r\left(k\right)+a_{1}0\le k\le k_1\end{array}$

By considering $k+1$ measurements, (71) can be written in vector form as follows

(72)$\begin{array}{c}\hfill {\mathbf{z}}_i^k=\left[\begin{array}{c}{\mathbf{z}}_r^k\underset{¯}{1}\end{array}\right]{\mathbf{z}}_d\\ \hfill {\mathbf{z}}_i^k={\mathbf{A}}_d^k{\mathbf{b}}_d\end{array}$

where subscript d denotes the Diffusion Arc and

(73)$\begin{array}{c}\hfill {\mathbf{z}}_i^k=\left[\begin{array}{c}{z}_i\left(0\right)\\ z_i\left(1\right)\\ ⋮\\ z_i\left(k\right)\end{array}\right],{\mathbf{z}}_r^k=\left[\begin{array}{c}{z}_r\left(0\right)\\ z_r\left(1\right)\\ ⋮\\ z_r\left(k\right)\end{array}\right],\underset{¯}{1}=\left[\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}\right],{\mathbf{b}}_d=\left[\begin{array}{c}{a}_0\\ a_1\end{array}\right]\end{array}$

The least-square estimate of ${\mathbf{b}}_d$ is obtained as

(74)$\begin{array}{c}\hfill {\widehat{\mathbf{b}}}_d=\left[\begin{array}{c}\widehat{a_0}\\ \widehat{a_1}\end{array}\right]={\left({\left({\mathbf{A}}_d^k\right)}^T{\mathbf{A}}_d^k\right)}^{-1}{\left({\mathbf{A}}_d^k\right)}^T{\mathbf{z}}_i^k\end{array}$

The sum of fitting error squares is

(75)$\begin{array}{c}\hfill s_r^k=\sum _{n=0}^k{\left(z_i\left(n\right)-\widehat{a_1}-\widehat{a_0}{z}_r\left(n\right)\right)}^2\end{array}$

The sum of error squares around the mean is

(76)$\begin{array}{c}\hfill s_t^k=\sum _{n=0}^k{\left(z_i\left(n\right)-{\overline{z}}_i\right)}^2\end{array}$

where $\overline{z_i}=\frac{1}{k+1}{\sum }_{n=0}^k{z}_i\left(n\right)$.

The correlation of fitting is given by

(77)$\begin{array}{c}\hfill r^k=\sqrt{\left|\frac{s_t^k-s_r^k}{s_t^k}\right|}\end{array}$

Figure 4 shows $r^k$ against k for two different cases of SNR.

It must be expected that $r^k$ starts to drop for $k>k_1$ because the linearity is lost after $k>k_1$. Considering the fact that $r^k$ demonstrates how the data fit the selected model, a threshold of 98% on $r^k$ is used to estimate the value of $k_1$.

7.2. Estimation of $k_2$

The next change in the Nyquist plot occurs between CT Arc and SEI Arc at ${\omega }_{k_2}$. When estimating $k_2$, the semicircle nature of the CT Arc is taken into account. Accordingly, the impedance measurements are in agreement with the following quadratic formula

(78)$\begin{array}{c}\hfill z_i^2\left(k\right)=a_0{z}_r^2\left(k\right)+a_1{z}_r\left(k\right)+a_2{k}_1<k\le k_2\end{array}$

By considering $k+1$ measurements, (78) can be written in vector form as follows

(79)$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathbf{z}}_i^k=\left[\begin{array}{c}{\left({\mathbf{z}}_r^k\right)}^2{\mathbf{z}}_r^k\underset{¯}{1}\end{array}\right]{\mathbf{b}}_c\\ \hfill {\mathbf{z}}_i^k={\mathbf{A}}_c^k{\mathbf{b}}_c\end{array}\end{array}$

where subscript c denotes the CT Arc and

(80)$\begin{array}{c}\hfill {\mathbf{z}}_i^k=\left[\begin{array}{c}{z}_i\left(0\right)\\ z_i\left(1\right)\\ ⋮\\ z_i\left(k\right)\end{array}\right],{\mathbf{z}}_r^k=\left[\begin{array}{c}{z}_r\left(0\right)\\ z_r\left(1\right)\\ ⋮\\ z_r\left(k\right)\end{array}\right],\underset{¯}{1}=\left[\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}\right],{\mathbf{b}}_c=\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right]\end{array}$

The least-square estimate of ${\mathbf{b}}_c$ is obtained as

(81)$\begin{array}{c}\hfill {\hat{\mathbf{b}}}_c=\left[\begin{array}{c}\widehat{a_0}\\ \widehat{a_1}\\ \widehat{a_2}\end{array}\right]={\left({\left({\mathbf{A}}_c^k\right)}^T{\mathbf{A}}_c^k\right)}^{-1}{\left({\mathbf{A}}_c^k\right)}^T{\mathbf{z}}_i^k\end{array}$

The sum of fitting error squares is

(82)$\begin{array}{c}\hfill s_r^k=\sum _{n=0}^k{\left(z_i^2\left(n\right)-\widehat{a_2}-\widehat{a_1}{z}_r\left(n\right)-\widehat{a_0}{z}_r^2\left(n\right)\right)}^2\end{array}$

The sum of error squares around the mean is

(83)$\begin{array}{c}\hfill s_t^k=\sum _{i=0}^k{\left(z_i^2\left(n\right)-\overline{{\left(z_i\right)}^2}\right)}^2\end{array}$

where $\overline{{\left(z_i\right)}^2}=\frac{1}{k+1}{\sum }_{n=0}^k{z}_i^2\left(n\right)$.

The correlation of fitting is given by

(84)$\begin{array}{c}\hfill r^k=\sqrt{\left|\frac{s_t^k-s_r^k}{s_t^k}\right|}\end{array}$

Figure 5 shows $r^k$ against k for two different cases of SNR.

It is obvious that at the first data points on the left side of the $k_1$ point in the Nyquist plot, the measured data does not fit a circular model. The point at which $r^k$ starts to increase is approximately where the CT Arc fits a circular model and is labelled as $k_{2,low}$. In addition, it must be expected that $r^k$ starts to drop for $k>k_2$ because the CT Arc semi-circle ends at $k_2$. With this consideration, $k_2$ is estimated. The points between $k_{2,low}$ and $k_2$ are used for CT-Arc parameter estimation.

7.3. Estimation of $k_3$, $k_4$

Here, $k_3$ is easily found as the x-intercept of the Nyquist plot, i.e, if $z_i\left(k\right)\ge 0$ and $z_i(k+1)<0$ then $k=k_3$.

Finally, considering the fact that ${\omega }_{k_4}$ is the maximum frequency used in the Nyquist plot, $k_4$ will be the number of frequencies through which the Nyquist plot is made. Figure 6 shows how the extracted $k_3$ and $k_4$ are used to detect the SEI Arc and the RL Arc.

In summary, this section presented the details of the proposed approach in this paper for ECM parameter estimation based on EIS. Here, curve-fitting approaches are used to detect specific geometric sections of the Nyquit spectrum and to estimate ECM parameters. The next section presents detailed performance analysis of the proposed approach through objective comparisons.

8. Results

8.1. Demonstration Using Simulated Data

In order to obtain the impedance spectrum, simulation was performed for a 1000mAh lithium-ion battery, the parameters of which are taken from [22]. For a better comparison with the previous approach [22], a 200 mA (C/5 Rate) DC charging current and an AC perturbation current with a peak of 70 mA were used to be the same as that used in [22]. Perturbation signals can be in different forms, e.g., square wave, multi-sine wave, linear or logarithmic chirps. In this paper, a series of 901 single-sine waves was used with different frequencies in the range of 0.01 Hz to 10 kHz ($k=$ 901 frequencies) for better resolution and a more precise demonstration of an industrial LIBs. The battery charging current, $I_{\mathrm{B}}\left(t\right)$, is represented as

(85)$\begin{array}{c}\hfill I_{\mathrm{B}}\left(t\right)=\left\{\begin{array}{cc}{I}_{\mathrm{dc}}+I_{\mathrm{m}}sin\left({\omega }_{1}t\right),\hfill & t<t_1\hfill \\ I_{\mathrm{dc}}+I_{\mathrm{m}}sin\left({\omega }_{2}t\right),\hfill & t_1\le t<t_2\hfill \\ ⋮\hfill \\ I_{\mathrm{dc}}+I_{\mathrm{m}}sin\left({\omega }_{k}t\right),\hfill & t_{k-1}\le t<t_k\hfill \end{array}\right.\end{array}$

where $I_{\mathrm{dc}}$ is the DC current, $I_{\mathrm{m}}$ is the peak of the perturbation current and each $t_k$ is selected such that $t_k-t_{k-1}=\frac{1}{f_k}$, i.e., each frequency is made to have one full cycle of data.

To measure the noise effects, it is assumed that the measured current and voltage are corrupted with Gaussian noise of the same noise variances, i.e.,

(86)$\begin{array}{c}\hfill {\sigma }_v={\sigma }_i=I_m{10}^{\left(-\frac{\mathrm{SNR}}{20}\right)}\end{array}$

where SNR varies from 0 to 50 dB (0, 10, 20, 30, 40, 50 dB). The Nyquist plot was derived for each case (each level of noise) and both the previous approach [32] and the proposed algorithm in this article were applied to estimate the ECM parameters. The performance of each algorithm is quantified in terms of the normalized percentage mean square error, simply referred to hereafter as Error (%). For example, the Error (%) of estimating the Ohmic resistance by the proposed approach is defined as

The Nyquist plot is derived for each case (each level of noise) and algorithms explained in Section 3, Section 4 and Section 6 are applied to estimate ECM parameters. The performance of each algorithm is quantified in terms of the normalized percentage mean square error, simply referred to hereafter as Error%. For example, the Error (%) of estimating the Ohmic resistance by the proposed approach is defined as

(87)$\begin{array}{c}\hfill \mathrm{Error}(\%)=\frac{\left|R_{\mathsf{\Omega }}-{\widehat{\widehat{R}}}_{\mathsf{\Omega }}\right|}{R_{\mathsf{\Omega }}}100\end{array}$

Each reported error measure was averaged over 100 Monte-Carlo runs.

The results of the approach proposed in this paper is compared with the previous approach, presented in [32]. This approach is briefly described in Section 4. In the comparison figures presented in this section, the previous approach [32] is labelled as “Previous Approach”.

Figure 7, Figure 8, Figure 9 and Figure 10 present the percentage error for each of the approaches discussed in this paper. Figure 7 compares the performance in estimating the Warburg coefficient ($\sigma$). It shows that the proposed algorithm performs better than the previous approach and non-linear LS estimation at all SNR values. An explanation of the performance loss by the previous approach could be that it used only two points to find the slope of a line to estimate $\sigma$. On the other hand, the proposed algorithm used many pairs of points and resulted in improved estimation errors. Moreover, the non-linear LS approach shows a significant increase in error with respect to SNR. The reason for the failure of the non-linear LS approach could be attributable to the severe non-linearity in the model when it comes to estimating the Warburg coefficient. In addition, it is worth mentioning that the previous approach’s performance in estimating $\sigma$ is highly dependent on the arbitrary frequency selected, but the proposed algorithm is more stable.

Figure 8a summarizes the performance comparison for estimating $R_{\mathrm{CT}}$. The proposed algorithm not only outperforms the previous approach at low SNRs but also shows about 10% lower error at high-SNR regions. Once again, a possible explanation for the loss of performance by the previous approach is that it relied on peak detection and just one single data point for the estimation of two parameters. Figure 8a shows the non-linear LS approach outperforming both the previous approach and the proposed approach in estimating $R_{\mathrm{CT}}$ and $C_{\mathrm{DL}}$. However, when seen as a whole, the advantage of the non-linear LS algorithm is severely lost in the estimation of some other parameters.

Figure 9a shows that, despite the high accuracy of the previous approach in estimating $R_{\mathrm{SEI}}$ at high SNRs, it severely lags behind the proposed algorithm in the presence of a high level of noise. The proposed algorithm estimates $C_{\mathrm{SEI}}$ with lower errors when the SNR is less than 10 dB compared to the previous approach. When the SNRs are higher, however, the methods are very close. This is illustrated in Figure 9b. Moreover, the non-linear LS approach performs poorly in estimating $R_{\mathrm{SEI}}$. However, the non-linear LS algorithms perform exceptionally well in the estimation of $C_{\mathrm{SEI}}$ compared to both the previous and the proposed approaches.

All the algorithms have excellent performance (error rate lower than 1%) in estimating $R_{\mathsf{\Omega }}$. However, the proposed approach has even better performance at lower SNR (error rate lower than 0.2%) compared to the previous approach. The non-linear LS approach performs the same at all SNRs, making it best at lower SNR and worst at higher SNR. Although the proposed algorithm uses the same concept from the previous approach in estimating RL Arc parameters, the reason it works better than the previous approach is the simple averaging of the data between $k_3$ and $k_4$. This is also obvious from Figure 10b, which shows that in the presence of a high level of noise, the proposed algorithm performs much better than the previous approach; however, both approaches are affected by estimation bias at high SNR.

For a better comparison, the extracted parameters from each method were used to generate the Nyquist plot, in order to compare how well the two methods perform in following the measurements. Figure 11a compares the performance of the previous approach and the proposed algorithm in estimating the Nyquist plot at high SNR (SNR = 30 dB). It can be seen clearly that although the previous method (solid blue line) estimated the battery parameters in an acceptable range, the Nyquist plot from the proposed approach (red line) is much closer to the measurement (purple line). However, what makes the proposed method surpass the previous approach is its accuracy in estimating the battery parameters in the presence of a high level of measurement noise. As shown in Figure 11b, the proposed algorithm (red line) performed impressively in estimating the battery parameters from severely corrupted signals (purple line), whereas the Nyquist plot generated by the previous method (blue line) failed to follow the measurement (purple line). In the presence of a high level of noise, the previous approach not only estimates the Warburg and CT Arc with high error rates but also fails to detect the SEI Arc, a lack of which can be observed in Figure 11b (dashed red line).

8.2. Demonstration Using Real Data

Figure 12 shows the experimental setup consisting of an Arbin battery cycler with Gamry interface 5000P EIS device. A relatively new cylindrical LIB (LG INR18650 MJ1) is used in the experiment. The specifications of the battery are shown in Table 1. The SOC of the battery was nearly empty when the experiment was performed. The Gamry EIS devise used zero DC current with a 50mA sinusoid superimposed on it. In a second experiment, right after the first experiment, the experiment was repeated with a 200 mA charging current and a 50 mA sinusoid superimposed on it. The output of the Gamry EIS device is the real and imaginary values of the measured impedance, as shown in Figure 13. The parameter-estimation algorithm proposed in Section 3 and Section 6 was applied to the impedance measurement to estimate the battery AR-ECM parameters. Table 2 shows the estimated values of the AR-ECM parameters based on the nonlinear LS approach. It can be seen in Table 2 that the estimated parameters are slightly different when there is a DC current compared to the zero-mean current experiment. It can also be concluded that the estimates with DC current could be closer to the true value. However, simulation studies showed the ineffectiveness of the nonlinear LS algorithm in the estimation of some of the ECM parameters. Table 3 summarizes the results of ECM parameter estimation using the approach proposed in this paper. As expected, the estimated values are closer to the ones obtained through the nonlinear LS approach.

Figure 14 compares the Nyquist plots of the estimated parameters from the proposed and non-linear approaches to the real data collected at zero-DC current and nonzero-DC current.

8.3. Discussion of Computational Complexity

In this section, the computational performance of the proposed algorithm is discussed. The algorithms were implemented on a Ryzen 7 5800H 8-core CPU which had 4.4 GHz max speed. The operating system was Windows 11. Table 4 shows the average computation time taken by the non-linear LS, previous approach and proposed approach for one run. Given the Nyquist plot, Table 4 shows the time taken by the respective algorithms to compute the parameters.

9. Conclusions

This paper presented an approach to estimate the equivalent-circuit-model (ECM) parameters of a battery based on electrochemical impedance spectroscopy (EIS). The proposed approach exploits the expected geometric features of the Nyquist plot of battery impedance for ECM parameter estimation. The Nyquist plot of the battery is divided into four sections, namely, the ‘Diffusion Arc’, ‘CT Arc’, ‘SEI Arc’ and the ‘RL Arc’; each arc is shown to correspond to different ECM parameters. First, an approach based on the correlation of curve fitting was used to detect the boundaries of each of these four arcs in the Nyquist plot. Once these arcs are detected, specific parameter-estimation approaches, based on least squares estimation, was used to estimate the ECM parameters of the battery. The proposed approach was shown to be superior to an existing method. Further, a non-linear least squares estimation approach was implemented to evaluate the proposed approach. Compared to the non-linear estimation approach, the proposed approach was found to have comparable performance, robustness, and computational simplicity.

There are several limitations in the present work which need to be addressed by future researchers. First, the performance of the present work was not compared against any theoretical error bounds. Indeed, a theoretical performance bound for the ECM parameters of a battery based on EIS is still to be derived. One of the observations from the present work is that the non-linear estimation approaches offer better estimation accuracy for some parameters, whereas the proposed approach offers robustness. Hence, a hybrid approach could be developed to combine the benefits of both a non-linear approach and the proposed linear approach in this paper. Further, EIS techniques, in general, are time-consuming, especially at the low-frequency spectrum, i.e., in estimating the parameters of the Diffusion Arc; hence, further techniques must be explored to rapidly estimate the parameters of the diffusion arc. One of the advantages of the proposed approach is that the parameter estimation for each arc can be performed separately; for this, an approach must be developed to detect each arc of the Nyquist plot separately so that the parameter-estimation approach can be decoupled.

Author Contributions

Conceptualization, B.B. and K.R.P.; methodology, B.B. and K.R.P.; software, M.A.; validation M.A.; formal analysis, M.A. and B.B.; data curation, M.A.; writing—original draft preparation, B.B., and M.A.; writing—review and editing, B.B., and M.A.; visualization, B.B and M.A.; supervision, B.B. and B.S.; project administration, B.B.; funding acquisition, B.B. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. AR-ECM and related impedance spectrum. (a) Adaptive Randles equivalent circuit model (AR-ECM). (b) Impedance spectrum (Nyquist Plot).

Figure 2. Effect of approximation on performance. The red line shows the Nyquist plot in the absence of Warburg impedance and the dashed line shows it when the Warburg impedance is present. The presence of Warburg impedance affected the frequency at the peak; this will affect the estimated value of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] (see (24)).

Figure 3. Effect of measurement noise on performance. Blue ‘.’ marks show the impedance measurement at various SNR levels. The dashed line shows the re-created Nyquist plot using the estimated values based on the previous approach (Section 4).

Figure 4. Diffusion-arc detection. The top plots show the correlation of fitting (77) along with the threshold. The plots at the bottom show the detected Diffusion Arc in red. (a) Low SNR. (b) High SNR.

Figure 5. CT Arc detection. The top plots show the correlation fitting of the (84) along with threshold. The plots at the bottom show the detected CT arc in red. (a) Low SNR. (b) High SNR.

Figure 6. SEI Arc and RL Arc detection. The top plots show the detected SEI arc. The plots at the bottom depict the detected RL arc. (a) Low SNR. (b) High SNR.

Figure 7. Parameter estimation of the Diffusion Arc. Blue line shows the % error in parameters estimation using previous approach. % error in parameters estimation using the proposed approach and non-linear approach is shown in red and yellow.

Figure 8. Parameter estimation of the CT Arc. Blue line shows the % error in parameters estimation using previous approach. % error in parameters estimation using the proposed approach and non-linear approach is shown in red and yellow. (a) Charge transfer (CT) resistance. (b) Double-layer (DT) capacitance.

Figure 9. Parameter estimation of the SEI Arc. Blue line shows the % error in parameters estimation using previous approach. % error in parameters estimation using the proposed approach and non-linear approach is shown in red and yellow. (a) Solid electrolyte interface (SEI) resistance. (b) SEI capacitance.

Figure 10. Parameter estimation of the RL Arc. Blue line shows the % error in parameters estimation using previous approach. % error in parameters estimation using the proposed approach and non-linear approach is shown in red and yellow. (a) Ohmic resistance. (b) Stray inductance.

Figure 11. Performance comparison. (Estimated parameters are used to reconstruct Nyquist plot). (a) High SNR. (b) Low SNR.

Figure 12. Experimental setup.

Figure 13. Nyquist plot for a real battery. The Nyquist plot in red depicts the impedance response of the battery with 0.2 A DC current. The impedance response of the battery with 0 DC current is shown in blue.

Figure 14. (a) Zero DC current: This figure compares the Nyquist plots of estimated parameters from the proposed and non-linear approaches to the real data (b) Nonzero DC current: this figure compares the Nyquist plots of estimated parameters from the proposed and non-linear approaches to the real data.

Table 1

Battery specifications.

Specification	Value (unit)
$\mathrm{Nominal}\mathrm{capacity}$	$3500\mathrm{mAh}$
$\mathrm{Max}.\mathrm{current}$	$10$
$\mathrm{Nominal}\mathrm{voltage}\left({\mathrm{V}}_{\mathrm{nom}}\right)$	$3.7$
$\mathrm{Height}$	$65\mathrm{mm}$
$\mathrm{Diameter}$	$18\mathrm{mm}$
$\mathrm{Weight}$	$46.5$

Table 2

Estimated parameters (nonlinear LS).

Paramter	Value (0 DC Current)	Value (−0.2 DC Current)
${\mathrm{R}}_{\mathsf{\Omega }}$	$412.5\mathrm{m}\mathsf{\Omega }$	$411.3\mathrm{m}\mathsf{\Omega }$
${\mathrm{R}}_{\mathrm{CT}}$	$2.6848\mathrm{m}\mathsf{\Omega }$	$3.22\mathrm{m}\mathsf{\Omega }$
${\mathrm{C}}_{\mathrm{DL}}$	$1.000\mathrm{F}$	$1.000\mathrm{F}$
$\mathrm{L}$	$1.1588\times {10}^{-6}\mathrm{H}$	$1.1534\times {10}^{-6}\mathrm{H}$
$\sigma$	$0.0584$	$0.0656$

Table 3

Estimated parameters (proposed approach).

Paramter	Value (0 DC Current)	Value (−0.2 DC Current)
${\mathrm{R}}_{\mathsf{\Omega }}$	$410.7\mathrm{m}\mathsf{\Omega }$	$409.7\mathrm{m}\mathsf{\Omega }$
${\mathrm{R}}_{\mathrm{CT}}$	$3.1\mathrm{m}\mathsf{\Omega }$	$3.3\mathrm{m}\mathsf{\Omega }$
${\mathrm{C}}_{\mathrm{DL}}$	$1.0278\mathrm{F}$	$0.8740\mathrm{F}$
$\mathrm{L}$	$1.1142\times {10}^{-6}\mathrm{H}$	$1.1107\times {10}^{-6}\mathrm{H}$
$\sigma$	$0.0020$	$0.0041$

Table 4

Comparison of computational time.

Non-Linear LS (Section 3)	Previous (Section 4)	Proposed (Section 6)
2.188	0.051	0.081

References

1. Vezzini, A. Lithium-ion battery management. Lithium-Ion Batteries; Elsevier: Amsterdam, The Netherlands, 2014; pp. 345-360.

2. Hu, X.; Zou, C.; Zhang, C.; Li, Y. Technological developments in batteries: A survey of principal roles, types, and management needs. IEEE Power Energy Mag.; 2017; 15, pp. 20-31. [DOI: https://dx.doi.org/10.1109/MPE.2017.2708812]

3. Pattipati, B.; Balasingam, B.; Avvari, G.; Pattipati, K.; Bar-Shalom, Y. Open circuit voltage characterization of lithium-ion batteries. J. Power Sources; 2014; 269, pp. 317-333. [DOI: https://dx.doi.org/10.1016/j.jpowsour.2014.06.152]

4. Lai, X.; Gao, W.; Zheng, Y.; Ouyang, M.; Li, J.; Han, X.; Zhou, L. A comparative study of global optimization methods for parameter identification of different equivalent circuit models for Li-ion batteries. Electrochim. Acta; 2019; 295, pp. 1057-1066. [DOI: https://dx.doi.org/10.1016/j.electacta.2018.11.134]

5. Pastor-Fernández, C.; Widanage, W.; Chouchelamane, G.; Marco, J. A SoH diagnosis and prognosis method to identify and quantify degradation modes in Li-ion batteries using the IC/DV technique. Proceedings of the 6th Hybrid and Electric Vehicles Conference (HEVC 2016); London, UK, 2–3 November 2016.

6. Zhang, D.; Dey, S.; Perez, H.E.; Moura, S.J. Remaining useful life estimation of lithium-ion batteries based on thermal dynamics. Proceedings of the IEEE 2017 American Control Conference (ACC); Seattle, WA, USA, 24–26 May 2017; pp. 4042-4047.

7. Balasingam, B.; Avvari, G.; Pattipati, B.; Pattipati, K.; Bar-Shalom, Y. A robust approach to battery fuel gauging, part I: Real time model identification. J. Power Sources; 2014; 272, pp. 1142-1153. [DOI: https://dx.doi.org/10.1016/j.jpowsour.2014.07.034]

8. Kumar, M.N.S.; Balakrishnan, K. Functional Safety Development of Battery Management System for Electric Vehicles. Proceedings of the IEEE 2019 IEEE Transportation Electrification Conference (ITEC-India); Bengaluru, India, 17–19 December 2019; pp. 1-6.

9. Liu, X.; Hui, S.R. Equivalent circuit modeling of a multilayer planar winding array structure for use in a universal contactless battery charging platform. IEEE Trans. Power Electron.; 2007; 22, pp. 21-29. [DOI: https://dx.doi.org/10.1109/TPEL.2006.886655]

10. Andre, D.; Meiler, M.; Steiner, K.; Wimmer, C.; Soczka-Guth, T.; Sauer, D. Characterization of high-power lithium-ion batteries by electrochemical impedance spectroscopy. I. Experimental investigation. J. Power Sources; 2011; 196, pp. 5334-5341. [DOI: https://dx.doi.org/10.1016/j.jpowsour.2010.12.102]

11. Lai, X.; Zheng, Y.; Sun, T. A comparative study of different equivalent circuit models for estimating state-of-charge of lithium-ion batteries. Electrochim. Acta; 2018; 259, pp. 566-577. [DOI: https://dx.doi.org/10.1016/j.electacta.2017.10.153]

12. He, H.; Xiong, R.; Fan, J. Evaluation of lithium-ion battery equivalent circuit models for state of charge estimation by an experimental approach. Energies; 2011; 4, pp. 582-598. [DOI: https://dx.doi.org/10.3390/en4040582]

13. Mousavi, S.; Nikdel, G. Various battery models for various simulation studies and applications. Renew. Sustain. Energy Rev.; 2014; 32, pp. 477-485. [DOI: https://dx.doi.org/10.1016/j.rser.2014.01.048]

14. He, H.; Xiong, R.; Guo, H. Online estimation of model parameters and state-of-charge of LiFePO4 batteries in electric vehicles. Appl. Energy; 2012; 89, pp. 413-420. [DOI: https://dx.doi.org/10.1016/j.apenergy.2011.08.005]

15. Balasingam, B.; Pattipati, K.R. On the identification of electrical equivalent circuit models based on noisy measurements. IEEE Trans. Instrum. Meas.; 2021; 70, pp. 1-16. [DOI: https://dx.doi.org/10.1109/TIM.2021.3068171]

16. Heaviside, O. Electrical Papers by Oliver Heaviside (In Two Volumes); Macmillon Co: New York, NY, USA, 1894; Volume 268, pp. 68-73.

17. Orazem, M.E.; Tribollet, B. Electrochemical Impedance Spectroscopy; John Wiley & Sons: Hoboken, NJ, USA, 2008; pp. 383-389.

18. Islam, S.R.; Park, S.Y. Precise online electrochemical impedance spectroscopy strategies for Li-ion batteries. IEEE Trans. Ind. Appl.; 2019; 56, pp. 1661-1669. [DOI: https://dx.doi.org/10.1109/TIA.2019.2958555]

19. Waag, W.; Käbitz, S.; Sauer, D.U. Experimental investigation of the lithium-ion battery impedance characteristic at various conditions and aging states and its influence on the application. Appl. Energy; 2013; 102, pp. 885-897. [DOI: https://dx.doi.org/10.1016/j.apenergy.2012.09.030]

20. Choi, W.; Shin, H.C.; Kim, J.M.; Choi, J.Y.; Yoon, W.S. Modeling and applications of electrochemical impedance spectroscopy (EIS) for lithium-ion batteries. J. Electrochem. Sci. Technol.; 2020; 11, pp. 1-13. [DOI: https://dx.doi.org/10.33961/jecst.2019.00528]

21. Zhang, Y.; Wang, C.Y. Cycle-life characterization of automotive lithium-ion batteries with LiNiO2 cathode. J. Electrochem. Soc.; 2009; 156, A527. [DOI: https://dx.doi.org/10.1149/1.3126385]

22. Islam, S.R.; Park, S.Y.; Balasingam, B. Circuit parameters extraction algorithm for a lithium-ion battery charging system incorporated with electrochemical impedance spectroscopy. Proceedings of the 2018 IEEE Applied Power Electronics Conference and Exposition (APEC) IEEE; San Antonio, TX, USA, 4–8 March 2018; pp. 3353-3358.

23. Pattipati, B.; Sankavaram, C.; Pattipati, K. System identification and estimation framework for pivotal automotive battery management system characteristics. IEEE Trans. Syst. Man Cybern. Part C (Appl. Rev.); 2011; 41, pp. 869-884. [DOI: https://dx.doi.org/10.1109/TSMCC.2010.2089979]

24. Pastor-Fernández, C.; Widanage, W.D.; Marco, J.; Gama-Valdez, M.Á.; Chouchelamane, G.H. Identification and quantification of ageing mechanisms in Lithium-ion batteries using the EIS technique. Proceedings of the 2016 IEEE Transportation Electrification Conference and Expo (ITEC); Dearborn, MI, USA, 27–29 June 2016; pp. 1-6.

25. Santoni, F.; De Angelis, A.; Moschitta, A.; Carbone, P. Analysis of the Uncertainty of EIS Battery Data Fitting to an Equivalent Circuit Model. Proceedings of the 2021 IEEE 6th International Forum on Research and Technology for Society and Industry (RTSI); Naples, Italy, 6–9 September 2021; pp. 497-501.

26. Li, H.; Lyu, Z.; Han, M. Robust and fast estimation of equivalent circuit model from noisy electrochemical impedance spectra. Electrochim. Acta; 2022; 422, 140474. [DOI: https://dx.doi.org/10.1016/j.electacta.2022.140474]

27. Wildfeuer, L.; Gieler, P.; Karger, A. Combining the distribution of relaxation times from EIS and time-domain data for parameterizing equivalent circuit models of lithium-ion batteries. Batteries; 2021; 7, 52. [DOI: https://dx.doi.org/10.3390/batteries7030052]

28. Allafi, W.; Uddin, K.; Zhang, C.; Sha, R.M.R.A.; Marco, J. On-line scheme for parameter estimation of nonlinear lithium ion battery equivalent circuit models using the simplified refined instrumental variable method for a modified Wiener continuous-time model. Appl. Energy; 2017; 204, pp. 497-508. [DOI: https://dx.doi.org/10.1016/j.apenergy.2017.07.030]

29. Zheng, Y.; Shi, Z.; Guo, D.; Dai, H.; Han, X. A simplification of the time-domain equivalent circuit model for lithium-ion batteries based on low-frequency electrochemical impedance spectra. J. Power Sources; 2021; 489, 229505. [DOI: https://dx.doi.org/10.1016/j.jpowsour.2021.229505]

30. Sangwan, V.; Sharma, A.; Kumar, R.; Rathore, A.K. Equivalent circuit model parameters estimation of li-ion battery: C-rate, soc and temperature effects. Proceedings of the 2016 IEEE International Conference on Power Electronics, Drives and Energy Systems (PEDES); Trivandrum, India, 14–17 December 2016; pp. 1-6.

31. Wang, Q.K.; He, Y.J.; Shen, J.N.; Hu, X.S.; Ma, Z.F. State of charge-dependent polynomial equivalent circuit modeling for electrochemical impedance spectroscopy of lithium-ion batteries. IEEE Trans. Power Electron.; 2017; 33, pp. 8449-8460. [DOI: https://dx.doi.org/10.1109/TPEL.2017.2780184]

32. Islam, S.; Park, S.Y.; Balasingam, B. Unification of Internal Resistance Estimation Methods for Li-Ion Batteries Using Hysteresis-Free Equivalent Circuit Models. Batteries; 2020; 6, 32. [DOI: https://dx.doi.org/10.3390/batteries6020032]

33. Pastor-Fernández, C.; Uddin, K.; Chouchelamane, G.H.; Widanage, W.D.; Marco, J. A comparison between electrochemical impedance spectroscopy and incremental capacity-differential voltage as Li-ion diagnostic techniques to identify and quantify the effects of degradation modes within battery management systems. J. Power Sources; 2017; 360, pp. 301-318. [DOI: https://dx.doi.org/10.1016/j.jpowsour.2017.03.042]

34. Taylor, S.; Gileadi, E. Physical interpretation of the Warburg impedance. Corrosion; 1995; 51, pp. 664-671. [DOI: https://dx.doi.org/10.5006/1.3293628]