Abstract
The increasing demand for efficient and reliable energy storage solutions has underscored the importance of accurate lithium-ion battery modelling. This project focuses on the estimation of parameters for lithium-ion battery equivalent circuit models from time-domain measurements. The generalized Randles circuit model and its extensions are explored and methods like vector fitting, transfer function estimation and grey-box estimation are used thereafter to perform the parameter estimation. A technique for directly estimating the parameters of the low frequency model including a constant phase element through grey-box estimation is devised. This method is then validated using experimental data from a battery module, achieving an RMSE of approximately \(\qty{4}{mV}\).
Introduction
The advent of portable electronic devices, electric vehicles (EVs), and battery energy storage systems (BESSs) in renewable energy systems has accelerated the demand for efficient and reliable energy storage solutions. Lithium-ion or Li-Ion batteries (LIBs), owing to their high energy density, long cycle life, and favorable performance characteristics, have emerged as the preferred choice for these applications. Accurate modelling of LIBs is crucial for enhancing their performance, extending their lifespan, and ensuring safety in their applications. [1] This process is pivotal for developing accurate battery management systems (BMS) that can predict the state of charge (SOC), state of health (SOH), and other critical parameters of the battery in real-time [2].
One of the fundamental approaches to modelling LIBs is through equivalent circuit models (ECMs). These models simplify the complex electrochemical processes within the battery into electrical components such as resistors, capacitors, and constant phase elements (CPEs). A lot of ECMs have been proposed in the past decades, such as simple internal resistance model, first-order and second-order RC model…, and these models are compared in [1]. The generalised Randles circuit, a widely recognized ECM, captures the essential dynamics of the battery, including charge transfer resistance, double-layer capacitance, and diffusion processes [3]. In [4], the identifiability of generalized Randles circuit model is studied and the findings indicate that it is structurally locally identifiable.
System identification and parameter estimation can be performed in both the time and frequency domains. For methods in the frequency domain, it is possible to excite the cell using sinusoidal signals at different frequencies and obtain the response at each frequency one by one for fitting (i.e., the principle of electrochemical impedance spectroscopy, or EIS), however, this process is relatively slow and is suitable for offline applications in the laboratory; it is also possible to excite the cell at one time using pseudo-random binary sequence or multi-sine signals that incorporate more than one frequency, however, in order to obtain the frequency response it is necessary to perform a discrete Fourier transform (DFT), which introduces estimation bias [5]. Direct parameter estimation in the time domain avoids this, while this places limitations on the components used in the ECM. For example, CPE is an ideal component characterized in the frequency domain, where its impedance is proportional to \(\mathrm{s}^\phi\). Here, \(\mathrm{s}\) represents the Laplace operator and \(\phi\) is a non-integer, making the system a fractional order system and thus challenging to deal in time domain. S.M.M. Alavi et al. employed a method known as simplified refined instrumental variable for continuous-time fractional systems, along with gradient-based optimization, to estimate the order of the fractional term [6].
The objective of this project is to explore the possibility of estimating the parameters of the Li-Ion cell equivalent circuit model from time-domain measurements by using different excitations and estimation methods. The focus is on the generalized Randles circuit model and its extensions, incorporating CPEs to account for low-frequency electrochemical behavior. In this work, we employ advanced techniques such as vector fitting and grey-box modelling to derive the parameters of the ECM from simulated and experimental data. The methodology involves system modelling and simulation, excitation signals generation, parameters estimation, results assessment and comparison, and model validation by experimental data.
The structure of this report is as follows: Section 2 introduces the topic and the tasks of this project. Section 3 details the methodology, including the ECMs used, the methods explored for parameter estimation, and the simulation setup. Section 4 presents the results of the parameter estimation from simulated data and Section 5 validated the model by experimental data. Finally, Section 6 offers a discussion on the findings and their implications, followed by the conclusion in Section 7.
Methodology
Equivalent Circuit Model
Randles Circuit Model
Figure 1 shows generalized Randles circuit model from [6], which has one resistor \(R_0\), one capacitor \({C_\text{w}}\) and several parallel RC branches in series. These components represent different electrochemical processes, \(R_0\) represents contact resistance, RC branches represent charge transfer and electric double layer resistance, and \({C_\text{w}}\) represents diffusion [7].
For the continuous-time linear system in Figure 1, the system input is defined as current \(i(t)\) , and output as port voltage \(v(t)\) and state \(\boldsymbol{x}(t)\). According to Kirchhoff’s voltage and current laws, the state-space model can be established as follows.
\[\boldsymbol{x}(t) = \begin{bmatrix} v_1(t)\\ v_2(t)\\ \vdots \\ v_n(t)\\ v_\text{w}(t) \end{bmatrix} \qquad \boldsymbol{u}(t) = i(t) \qquad \boldsymbol{y}(t) = v(t) \tag{1}\]
\[\begin{cases} \dot{\boldsymbol{x}} = \boldsymbol{Ax}+\boldsymbol{Bu}\\ \boldsymbol{y} = \boldsymbol{Cx} + \boldsymbol{Du} \end{cases} \tag{2}\]
\[\boldsymbol{A} = \begin{bmatrix} -\frac{1}{R_1C_1} & 0 & \ldots & 0 & 0\\ 0 & -\frac{1}{R_2C_2} & \ldots & 0 & 0\\ \vdots& \vdots & \ddots & \vdots &\vdots\\ 0 & 0 & \ldots & -\frac{1}{R_nC_n}& 0\\ 0 & 0 & \ldots & 0 & 0\\ \end{bmatrix} \quad \boldsymbol{B} = \begin{bmatrix} 1/C_1\\ {1}/{C_2} \\ \vdots\\ {1}/{C_n}\\ {1}/{C_\text{w}}\\ \end{bmatrix} \quad \boldsymbol{C} = \begin{bmatrix} 1 & 1 &\ldots&1 &1 \end{bmatrix} \quad \boldsymbol{D} = \begin{bmatrix} R_0 \end{bmatrix} \tag{3}\]
Based on Equation 3, the corresponding transfer function of the system is therefore built by \[ T(s) = \boldsymbol{C}(s\boldsymbol{I}-\boldsymbol{A})^{-1}\boldsymbol{B}+\boldsymbol{D} = \sum_{i=1}^n \frac{\displaystyle\frac{1}{C_i}}{s+\displaystyle\frac{1}{R_iC_i}} + \frac{\displaystyle\frac{1}{C_\text{w}}}{s} + R_0 \tag{4}\]
where \(s\) denotes Laplace operator. For the system, its residues are \(1/C_i\), poles are \(-1/\tau_i\), and direct term \(R_\infty\).
Constant Phase Element
Constant phase element, or CPE, is widely used to model non-ideal electrochemical behavior in batteries, which is defined as: \[ \bar{Z}_\text{CPE} = \frac{1}{(\mathrm{j}\omega)^\phi Q}\quad , \quad \phi\in(0,1) \] if \(\phi = 0\), it is a pure resistor and if \(\phi = 1\), it is a pure capacitor. In the Nyquist plot, it shows a straight line through the origin, which has constant phase characteristic.
However, the existence of fractional order \(\phi\) makes it difficult to simulate the properties of CPE. In [8], [9] , a novel method is discussed to decompose CPE into generalised Randles circuit as Figure 2, which makes it possible to study CPE in time-domain.
The method to decompose CPE is as follows:
Define CPE parameters \(Q\), \(\phi\), max phase ripple \(\Delta\phi\), frequency band of interest from \(\omega_0\) to \(\omega_1\) ;
calculate number of RC branches needed and ratio \(a, b\) ;
\[q = \frac{0.24}{1+\Delta\phi\cdot180/\pi} \qquad\text{and}\qquad n = \left\lceil \frac{\ln{\omega_0} - \ln{\omega_1}}{\ln{q} } \right\rceil\]
\[\begin{cases} a = 10^{-\phi\lg{q}}\\ b = q/a \end{cases}\]
Set initial values of \(R_1\), \(C_1\) and calculate all other values; \[\begin{cases} C_1 = 1\\ R_1 = 1/\omega_0 \end{cases} \Rightarrow\qquad \begin{cases} R_\text{inf} = R_1\displaystyle\frac{a^n}{1-a} \\ C_\text{inf} = C_1\displaystyle\frac{1-b}{b} \end{cases} \qquad \begin{cases} R_k = R_1a^{k-1}\\ C_k = C_1b^{k-1} \end{cases} \text{for } k = 1,2,...,n\]
Calculate gain correction factor \(g\) ;
\[\begin{cases} \omega_\text{avg} = \displaystyle\left(\frac{a}{b}\right)^{0.25}\cdot \frac{1}{C_1R_1q^{\lceil (n/2)-1\rceil}} \\ Z_\text{cpe}(\omega_\text{avg}) = R_\text{inf} + \displaystyle\frac{1}{\mathrm{j}\omega_\text{avg}C_\text{inf}} + \sum_{k=1}^{n} \frac{R_k}{1+\mathrm{j}\omega_\text{avg}R_kC_k} \end{cases} \]
\[ g = \frac{1}{Q\left|Z_\text{cpe}(\omega_\text{avg})\right| \omega_\text{avg}^\phi} \]
Correct all parameters by \(g\) . \[ \begin{cases} R_\text{inf} = R_\text{inf} \times g\\ C_\text{inf} = C_\text{inf} / g \end{cases}\quad \text{and}\qquad \begin{cases} R_k = R_k \times g\\ C_k = C_k / g \end{cases}\quad \text{for } k = 1,2,...,n \]
In this work, \(\Delta\phi\) is always specified as 0, and a MATLAB function is built to finish the steps above:
[Ri, Ci, Rinf, Cinf] = cpe_decomposer(Q,phi,f_0,f_1)where Q and phi are parameters of the CPE and f_0, f_1 indicate the frequency band of interest to decompose CPE (in \(\unit{\Hz}\)).
The following is an example of applying the method. Selecting the CPE parameters \(\phi = 0.50\), Q = \(\num{5.0E4} \mathrm{S\cdot (j\omega)^{0.50}}\), specifying the frequency band of interest as \(\qty{0.01}{\Hz}\) to \(\qty{100}{\Hz}\), and \(\Delta\phi = 0\), gives the decomposition parameters \(R_\mathrm{inf}= \qty{4.80e-7}{\ohm}\) and
\(C_\mathrm{inf}= \qty{4.58e5}{\farad}\), and RC branches in Table 1.
The comparison of the CPE and the decomposition model of this example is shown in Figure 3. The left plot is a Nyquist plot displaying the imaginary part of the impedance versus the real part (blue crosses for the CPE and red circles for the decomposition). The right plot presents the phase angle as a function of frequency. The phase angle starts at approximately \(\qty{-0.85}{\radian}\) at low frequencies, exhibits slight undulations, and then rises sharply at higher frequencies, approaching \(\qty{-0.65}{\radian}\) at around \(\qty{100}{\Hz}\). This suggests that the decomposition model closely matches the behavior of the CPE within the frequency range of our interest.
Low Frequency Model
Figure 4 shows a Nyquist plot of the impedance spectrum for a lithium-ion battery cell. The plot is divided into three regions: the ohmic region at high frequencies ( \(> \qty{1}{\kHz}\) ), the charge transfer region at intermediate frequencies (\(\qty{1}{\Hz} \sim \qty{1}{\kHz}\)), and the diffusion region at low frequencies ( \(\qty{1}{\Hz}\) ). The ohmic region is characterized by high frequencies where the impedance is mainly due to the ohmic resistance of the cell, which includes the resistance of the electrolyte, electrode materials, and current collectors. The charge transfer region corresponds to the intermediate frequencies where the impedance is influenced by the charge transfer resistance at the electrode/electrolyte interfaces, and is typically observed as a semicircular arc in the Nyquist plot. At low frequencies, the impedance is dominated by the diffusion of lithium ions within the electrodes and represented by a linear increase in the spectrum.
If only the impedance spectrum that is in the upper half of the y-axis is considered, in the ohmic region we can only get the intercept of the x-axis, i.e., a real impedance, which is modelled using a resistor; in the charge transfer region, the arc in the impedance spectrum can be modelled using one or more parallel RC branches in series; and in the diffusion region the linear growth is modelled using a constant-phase element (CPE). Thus, the equivalent circuit model built from impedance spectrum and its low frequency approximation are as shown in Figure 5.
As discussed in Section 3.1.2, with given max phase ripple \(\Delta\phi\) and frequency band of interest, we can decompose CPE into generalized Randles circuit in Figure 2. Therefore, the low frequency equivalent circuit model is constructed in Figure 6, where \(R_\text{sum} = R' + R_{\text{inf}}\) and \(R' = R_{\text{0}}+ R_\text{I} +R_\text{II}\).
Simulink
As analyzed in Section 3.1, all the circuit topologies in this work are the same as generalized Randles model, thus the circuit model shown in Figure 7 is built in Simulink. The excitation signal is generated in the workspace and then passed into the model and fed into the system by a controlled current source, all the resistor and capacitor values are defined in the workspace and the port current (input signal) and voltage (output signal) are sampled and sent back to the workspace.
Model Identification
Vector fitting
Vector fitting is a robust and widely used numerical technique for system identification and parameter estimation, particularly in the field of electrical engineering and control systems. It is designed to approximate complex frequency responses by rational transfer functions, making it highly suitable for modelling and analysing systems based on measured impedance or admittance data [10]. The primary goal of this method is to convert a set of frequency domain data points into a rational function model:
\[ G(s) \approx \sum_{n=1}^{N} \frac{c_n}{s-p_n} + d + se \] where \(G(s)\) is the transfer function, \(c_n\) are residues, \(p_n\) are poles, \(d\) is the direct term, and \(e\) is the slope of the polynomial.
The process begins with an initial guess of the poles of the system. These poles are typically distributed along the imaginary axis to cover the frequency range of interest. Then, a least squares problem is constructed to iteratively adjust the poles and residues and minimize the difference between the measured data and the model’s response. [10], [11], [12]
Estimate Transfer Function Models in Time Domain
tfest is a MATLAB function designed to estimate transfer function models from time-domain or frequency-domain data. The transfer function model represents the relationship between the input and output of a system as a ratio of polynomials in the Laplace domain. The core principle behind tfest involves fitting a transfer function model to the given data by solving an optimization problem.
The transfer function model is represented as a ratio of polynomials: \[
G(s) = \frac{B(s)}{A(s)} = \frac{b_0 + b_1 s + b_2 s^2 + ... + b_n s^n}{a_0 + a_1 s + a_2 s^2 + ... + a_m s^m}
\] The input-output data is organized in an iddata object and the user specifies the order of the model. MATLAB formulates an optimization problem where the objective is to minimize the error between the measured output and the model output. The error \(e(t)\) is defined as: \[
e(t)=y_\text{measured}(t)-y_\text{model}(t)
\tag{5}\]
The parameters of the transfer function (coefficients of \(A(s)\) and \(B(s)\)) are adjusted to minimize a cost function, typically the sum of squared errors (SSE): \[ \mathrm{SSE} = \sum_t [e(t)]^2 \tag{6}\]
Then, numerical optimization techniques, such as the Levenberg-Marquardt algorithm or other nonlinear least squares methods, are used to solve this optimization problem.
Estimate Grey-box Models in Time Domain
greyest is a MATLAB function used to estimate grey-box models, which are based on physical principles but include unknown parameters that need to be estimated from data. These models can capture the underlying physics of the system while accommodating measured data. Firstly, a predefined model structure need to be given by the user before the estimation. Typically, the grey-box model is defined by a function that outputs the state-space matrices \(\boldsymbol{A}(\boldsymbol{\theta})\), \(\boldsymbol{B}(\boldsymbol{\theta})\), \(\boldsymbol{C}(\boldsymbol{\theta})\), \(\boldsymbol{D}(\boldsymbol{\theta})\) based on a set of parameters \(\boldsymbol{\theta} = \{\theta_1, \theta_2, ... , \theta_n\}\). This function represents the physical relationships and dynamics of the system. Then, the input-output data is organized in an iddata object and the user specifies the initial guess of the parameters set \(\boldsymbol{\theta}_0\). Similarly, greyest formulates the same optimization problem by Equation 5 and Equation 6 to minimize the error between the measured and modelled outputs. This approach has been explored in [13] and [14] and gives good results.
Simulation Results
Estimate Generalised Randles Model by Vector Fitting
The parameter values of generalised randles model are arbitrarily selected and shown in Figure 8.
The multi-sine signal given in Equation 7 is used to excite the circuit in Figure 8, where \(w_i\) is the amplitude, \(f_i\) is the frequency and \(\varphi_i\) is the phase for each sinusoidal wave. If the sine waves have equal amplitudes \(w_i\), it is recommended to use a non-zero Schroeder phase \(\varphi_i\) to reduce the crest factor by [6], [15].
\[ u(t) = \sum_{i=1}^{n} w_i\cdot \sin{(2\pi f_i t + \varphi_i)} \tag{7}\]
To identify the system, the excitation signal frequency band should cover the dynamic range of the system, i.e., the minimum frequency should be less than \(1/\tau_\text{max}\), the maximum frequency should be greater than \(1/\tau_\text{min}\), and the test time should be longer than \(\tau_\text{max}\), wherer \(1/\tau_\text{max}\) and \(1/\tau_\text{min}\) are maximum time constant and minimum time constant respectively. In addition, the sampling frequency \(f_\text{s}\) should be at least twice the maximum excitation signal frequency according to the Nyquist-Shannon theory. The generated multi-sine signal is used as current input \(u(t)=i(t)\) and simulated in Simulink to get the output voltage signal \(y(t)=v(t)\). To study the effect of noise, a series of zero-mean noise are added to input and output signals.
As a method to fit in frequency domain, frequency responses are needed. Here, fft in MATLAB is used to transform the input and output signal from time domain to frequency domain:
\[\begin{cases} U(f) = \mathscr{F}(u(t))\\ Y(f) = \mathscr{F}(y(t)) \end{cases} \] where \(\mathscr{F}\) is the Fourier operator, and thus the frequency response is given by: \[H(f) = \frac{Y(f)}{U(f)}\] It is worth noting that since noise has been added to \(u(t)\) and \(y(t)\), only the frequency response of the \(n\) with the largest magnitude should be considered when calculating \(H(f)\). This helps to suppress the effects of noise because vector fitting is an extremely noise sensitive method [16] and using a frequency response that contains noise will not give a correct estimation.
In this example, \(\tau_\text{max} = R_0C_\text{w} = \qty{30}{s}\), \(\tau_\text{min} = R_\text{I}C_\text{I} = \qty{0.06}{s}\). Therefore, \(n=4\), \(w_i= \qty{1}{A}\), \(\varphi_i = \qty{1}{\radian}\), \(f_i = \qty{0.02}{\Hz}, \qty{0.2}{\Hz}, \qty{2}{\Hz}, \qty{20}{\Hz}\) for \(i= 1,...,4\) and test duration \(t= \qty{200}{s}\) are selected to generate the multi-sine excitation signal. Figure 9 shows the input signal (time domain and frequency domain), output signal (time domain) and frequency response of the system excited by this signal.
Table 2 presents the results of 10 simulations conducted with both noise-free data and noisy data, estimated using vector fitting. The values in the table are expressed in per unit values, where the base value for each column represents the respective true value. Thus, an estimation result equal to 1 signifies a perfect estimation. Across all simulations, the estimated parameters are close to 1 with very small bias. When comparing noise-free and noisy results, the mean values remain relatively stable, though the standard deviation increases in the presence of noise. Additionally, analysing results at different sampling rates reveals that while the mean values of the parameters show little variation, a lower sampling rate results in an increased standard deviation.
Estimate CPE Decomposition Model
Selecting the CPE parameters \(\phi = 0.62\), Q = \(\num{1.5E4}\, \mathrm{S\cdot (j\omega)^{0.62}}\), specifying the frequency band of interest as \(\qty{0.01}{\Hz}\) to \(\qty{100}{\Hz}\), and \(\Delta\phi = 0\), gives the decomposition parameters in Table 3.
Multi-sine Exciation
From Table 3, \(\tau_\text{max} = R_1C_1 = \qty{15.916}{s}\), \(\tau_\text{min} = R_4C_4 = \qty{0.22}{s}\). Therefore, \(n=7\), \(w_i= \qty{1}{A}\), \(\varphi_i = \qty{1}{\radian}\), \(f_i = \qty{0.05}{\Hz}, \qty{0.1}{\Hz}, \qty{0.2}{\Hz}, \qty{0.5}{\Hz}, \qty{1}{\Hz}, \qty{2}{\Hz}, \qty{5}{\Hz}\) for \(i= 1,...,7\) and test duration \(t=\) \(\qty{200}{s}\) are selected to generate the multi-sine excitation signal. Input current signal and output voltage signal in time domain are given by Simulink model and noise is then added. The noise introduced is zero-mean, normally distributed, and has a specific standard deviation. Unlike in Section 4.1, the small impedance of the CPE decomposition results leads to a voltage response at the level of several \(\qty{1e-4}{V}\), therefore, the noise standard deviations here are given as relative values. Specifically, the standard deviation, denoted as \(\sigma\), represents the standard deviations of the current and voltage noise as a proportion of their respective maximum magnitudes respectively.
After obtaining the noise-free or noisy input and output signals \(u(t)\) and \(y(t)\), the MATLAB function tfest is used to estimate the continuous time transfer function:
sys = tfest(u,y,np)where np is the number of poles and sys is the identified transfer function. Suppose the system residues and corresponding poles are denoted as \(r_i\) and \(p_i\), direct term as \(k\), where \(i= 1,...,7\). Specify the following order:
\[ \begin{cases} p_i = 0 \qquad & \text{for}\quad i=n \\ |p_{i+1}| > |P_i| \qquad & \text{for}\quad i=1,2,...,n-1 \end{cases} \] from Equation 4, all components are reconstructed by: \[ \begin{cases} R_\text{inf} = k \\ C_\text{inf} = 1/r_n\\ C_i = 1/r_i \qquad & \text{for}\quad i=1,2,...,n-1 \\ R_i = -r_i/p_i \qquad & \text{for}\quad i=1,2,...,n-1 \end{cases} \]
Table 4 shows the results of 20 simulations at different sampling rates and different noise levels. Similarly, for comparing the estimation results of individual parameters, the data in the table are per unit values. Results with a per unit value less than 0.25 or greater than 4 are categorized as outliers. Comparing different parameters, estimation results of \(R_\text{inf}\) are accurate in all cases, while for the other parameters, although the mean values are acceptable, the presence of noise leads to an increase in the outliers and an increase in the standard deviation. By comparing the results for \(\qty{50}{Hz}\) and \(\qty{10}{Hz}\) sampling rates, the higher sampling rate reduces the number of outliers and concentrates the distribution of the results.
PRBS
As a commonly used excitation in system identification for parameter estimation, Pseudo-Random Binary Sequence (PRBS) is a type of deterministic signal that exhibits properties of randomness, making it highly suitable for exciting a system over a wide frequency range. This ensures that the system’s dynamic characteristics can be effectively captured during the identification process. In addition, PRBS can be exactly reproduced, facilitating repeatable experiments and validation, and the binary nature of PRBS provides a high signal-to-noise ratio, making it effective in noisy environments[2].
In MATLAB, PRBS is generated by the function:
prbs = frest.PRBS('Order',n,'NumPeriods',1,'Amplitude',Am,'Ts',Ts)where \(n\) is the order or the number of stage of the PRBS. If the number of periods is equal to 1, this will generate a sequence of length \(N\), \(N\) is defined by: \[N = 2^n-1\]
The power spectrum of a PRBS, as shown in Figure 10, features discrete frequencies with a magnitude envelope characterized by the sinc(x) function, reaching its first zero at the clock frequency \(f_\text{c}\) [2], [17]. Consequently, the PRBS is bandwidth-limited to a frequency below that of the clock. Furthermore, the frequency resolution is determined by \(N\) and \(f_\text{c}\). A rule of thumb to chose \(f_\text{c}\) is given in [2], [18]: \[ f_\text{c} = 2.5f_\text{max} \] where \(f_\text{max}\) is the highest frequency of interest.
In order to generate sufficiently long sequences, the order of the PRBS \(n\) should not be too small; and in order to reach too long test duration, a very high \(n\) should be avoided. It is suggested in [2] that \(n\) should in between 8 and 14. In addition, according to the Nyquist-Shannon theory, the sampling frequency \(f_\text{s}\) should be at least twice the clock frequency \(f_\text{c}\).
Here, \(f_\text{c} = \qty{5}{Hz}\), \(f_\text{s} = \qty{10}{Hz}\) and order \(n=12\) are selected to generate the PRBS signal to estimate the parameters of our CPE decomposition model in Table 2. The first 50 s of the input signal (current) and output signal (voltage) are shown in Figure 11.
Similarly, noise is added to the current and voltage signals in the same way as in Section 4.2.1 . The results of 10 simulations using the above PRBS for excitation and tfest for parameter estimation are listed in Table 5. In the noise-free case, the estimation results of all parameters are almost accurate ( these per unit values are very close to 1 ). As the noise level increases, the estimation results of \(C_\text{inf}\) and \(R_\text{inf}\) remain stable, while for other parameters, outliers appear and the standard deviation in each column increases.
Estimate Low Frequency Equivalent Circuit Model
As discussed in Section 3.1.3, CPE is decomposed using a predefined frequency range and a max phase ripple, and at low frequencies the equivalent circuit topology is shown in Figure 6, for which the model parameters are \[ \boldsymbol{\theta} = [R', Q, \phi] \]
where \(\phi\in(0,1)\), and \(\boldsymbol{A}(\boldsymbol{\theta})\), \(\boldsymbol{B}(\boldsymbol{\theta})\), \(\boldsymbol{C}(\boldsymbol{\theta})\), \(\boldsymbol{D}(\boldsymbol{\theta})\) are matrices depend on parameter vector \(\boldsymbol{\theta}\).
Use the same parameters in Table 3 to build CPE, i.e. \(\phi = 0.62\), Q = \(\num{1.5E4}\, \mathrm{S\cdot (j\omega)^{0.62}}\), the frequency band of interest from \(\qty{0.01}{\Hz}\) to \(\qty{100}{\Hz}\), \(\Delta\phi = 0\), select \(R'= \qty{1.615}{\mohm}\) to build the Simulink model. PRBS and the real current measured by BMS are subsequently used as input current to excite the Simulink model, respectively. Here, since the parameters we want to estimate is \(\boldsymbol{\theta}\) instead of all RC branches, for PRBS, the clock frequency of \(\qty{1}{\Hz}\), order=8, amplitude = \(\qty{40}{A}\) and test duration of \(\qty{254}{s}\) are used. The real current collected by BMS is shown in Figure 12, which has a sampling rate of \(\qty{1}{\Hz}\), and it can be observed from the frequency domain that the current is mainly distributed below \(\qty{0.1}{\Hz}\). This means that it is challenging to estimate the parameters using real current data.
After knowing the model structure, linear grey-box model is estimated by MATLAB. First, a grey-box model with identifiable parameters is built by idgrey:
init_sys = idgrey(odefun,parameters_initial,fcn_type)where odefun relates the model parameters \(\boldsymbol{\theta}\) to the state-space representation (\(\boldsymbol{A}(\boldsymbol{\theta})\), \(\boldsymbol{B}(\boldsymbol{\theta})\), \(\boldsymbol{C}(\boldsymbol{\theta})\), \(\boldsymbol{D}(\boldsymbol{\theta})\)), fcn_type indicates whether the model is continuous or discrete, and parameters_initial gives initial values of parameters \(\boldsymbol{\theta}_\text{init}\): \[\boldsymbol{\theta}_\text{init} = [R'_\text{init}, Q_\text{init}, \phi_\text{init}]\] With the initial model, the linear grey-box model sys is estimated using time-domain input and output data, sys is an identified model that has the same structure as the initial model init_sys.
sys = greyest(data,init_sys)If needed, for an identified model sys, the response can be easily simulated by:
y = sim(sys,udata)which returns the simulated response of sys using the input data udata.
The method is firstly used for noise-free data. Define the initial parameters: \(R'_\text{init} = \qty{10}{\mohm}\), \(Q_\text{init}= \num{1e4}\, \mathrm{S\cdot (j\omega)^{0.5}}\), \(\phi_\text{init} = {0.5}\) . The estimation results using PRBS and randomly selected real current data excitation are shown in Table 6. We can observe that, by using maximum likelihood estimation (MLE) greyest, all of the estimation results obtained are very accurate.
Then, by adding noise to input and output data, the possibility of applying this method in experiments is explored. Zero-mean normal distributed random sequences with standard deviation of \((\qty{0.1}{A})/3\) and \((\qty{1.5}{\mV})/3\) are superposed to current and voltage signals respectively. The average results of 10 times of simulation are shown in Table 7.
The effect of initial values selection on the estimation results (average results for 10 times of simulation) is shown in the Table 8, which uses the real currents from \(0\) to \(\qty{200}{s}\) as inputs, and the noise is added in the same way as in Table 7. Comparing different sets of initial parameters, the results are insensitive to variations in \(R'_\text{init}\) and \(\phi_\text{init}\), while too small \(Q_\text{init}\) can lead to wrong results. In practice, \(\phi_\text{init}:{0.5}\) is a reasonable guess, while \(Q_\text{init}\) may need to be adjusted based on the value of \(\hat{\phi}\). For example, if \(\hat{\phi}\) is close to 0 or 1, we may need to try another \(Q_\text{init}\) and re-estimate the model.
Model Validation
This work uses the current and voltage data collected by BMS to validate the low frequency equivalent circuit model, since the sampling frequency of the BMS is \(\qty{1}{\Hz}\). The data is reported while the BESS is providing frequency containment reserve (FCR) service. In addition, according to the BMS, the batteries are operating at a stable temperature of about \(\qty{47}{\degreeCelsius}\). Figure 13 shows the first 800 seconds of the current of a battery module and the voltage of each cell in that module, which consists of six cells connected in series. In the figure, a positive current indicates discharging and a negative current indicates charging, the current varies in the range of a few tens of amperes and the voltage fluctuates around \(\qty{3.7}{V}\) to \(\qty{3.8}{V}\).
Since the low frequency equivalent circuit model built in Section 3.1.3 does not consider any voltage sources, the first step in performing parameter estimation is to remove the open circuit voltage (OCV) of the cells from the voltage data. For the cells in this battery module, the red crosses in Figure 14 are known, and the OCV is a nonlinear function of the State of Charge (SoC), thus OCV can be obtained for different SoCs using piecewise linear interpolation.
Known that \(\mathrm{SoC_0} = 60\%\) and cell capacity \(C_n = \num{62.5}\, \mathrm{Ah}\), assume the cells in this module is well balanced, \(\mathrm{SoC_t}\) can be calculated by ampere-hour integral method: \[\mathrm{SoC_t} = \mathrm{SoC_0} - \frac{\int_0^t i( \tau ) \mathrm{d}\tau}{3600\cdot C_n}\]
Figure 15 shows the variation of the SoC in \(\qty{40000}{s}\) when providing FCR service. We can observe that the SoC fluctuates between 50% and 65%, with an average value of 59.25%. This implies that the battery is in a “quasi-steady state” over long time scales.
The method discussed in Section 4.3 can be then used. In init_sys, specify that the CPE is decomposed from \(\qty{0.01}{\Hz}\) to \(\qty{100}{\Hz}\) when creating odefun and perform the estimation using the current and voltage data from \(\qty{0}{s}\) to \(\qty{200}{s}\). However, after adjusting the initial parameters \(\boldsymbol{\theta}_\text{init}\), the parameters estimated are always \(\hat{R'} = \qty{0}{\ohm}\), \(\hat{Q} = \num{620}\, \mathrm{S\cdot (j\omega)^{0.071}}\) , and \(\hat{\phi} = 0.071\). By simulating the output of this identified system and superposing the removed OCV back to \(y(t)\), the comparison of real voltage and model output voltage is shown in Figure 16. The model tracks well with the real voltage during the \(\qty{200}{s}\) duration used for estimation but has large error thereafter.
Considering the estimation result \(\hat{\boldsymbol{\theta}}\), we can discover that \(\hat{\phi} \approx 0\), thus: \[\hat{Z} \approx \hat{R'} + \hat{Z}_\text{cpe} = \hat{R'} + \frac{1}{\hat{Q}} = \frac{1}{620\, \mathrm{S}} = 1.613\,\mathrm{m\Omega}\] which means the grey-box model considers the system to approximate a pure resistor.
The reason can be found by analysing the impedance spectrum of the CPE decomposition and the current data in Time Domain. Figure 17 shows the impedance spectrum of a CPE decomposed between \(\qty{0.01}{\Hz}\) and \(\qty{100}{\Hz}\), where the decomposition model fits the CPE very well in the selected range, while the region outside this frequency range has large errors. In contrast, Figure 12 shows that the current data is mainly DC and the frequency is distributed in the range below \(\qty{0.05}{\Hz}\). Therefore, the mismatch between the input current and the model produces incorrect estimation result.
Then, we choose to decompose the CPE in the range of \(\qty{0.1}{\mHz}\) to \(\qty{1}{\Hz}\) and re-establish the initial system init_sys. For parameter estimation, considering the “quasi steady state” of the battery in long time scales, instead of perform the parameter estimation in the whole duration, a series of windows are divided within the duration and moved in specific steps, and the parameter estimation is performed for each window subsequently. After that, the results of all the windows are averaged and given as the average parameter for the whole duration, which helps to enhance the stability of the estimation method.
For the selection of test duration, window length and step length, the test duration should not be too short to ensure “quasi steady state”, it should be at least longer than one hour; the window length should not be too short in order to capture the dynamics of the system, but to avoid changes in the state of the battery, the window should not be too long, generally should be in 20 minutes; the step length should be shorter than the length of the window so that different windows can overlap, generally take the length of a few minutes.
In this work, test duration of 0 to \(\qty{5000}{s}\), window length of \(\qty{1000}{s}\), and step of \(\qty{200}{s}\) are selected. The parameter estimation results for each window of cell 1 are listed in Table 9. In the estimation, the range of \(\phi\) is specified as 0.4 to 0.8, thus the windows with Index of 16 to 20 in the table enable the boundaries, which should be discarded when calculating the average parameter for the whole duration. The average estimated parameters of cell 1 for this duration are then calculated: \(\hat{\boldsymbol{\theta}}_1 = [\hat{R'}_1, \hat{Q}_1, \hat{\phi}_1\)], where \(\hat{R'}_1 = \qty{0.0012}{\ohm}\), \(\hat{Q}_1 = \num{7.983E4}\, \mathrm{S\cdot (j\omega)^{0.646}}\) , and \(\hat{\phi}_1 = 0.646\).
Fixing the parameters of cell 1 as \(\hat{\boldsymbol{\theta}}_1\), a identified system sys using data from \(0\) to \(\qty{5000}{s}\) is obtained. Thereafter the system input data (current) is extended to \(\qty{30000}{s}\) and the model output (voltage) is simulated, as shown in Figure 18. The blue line in the figure represents the real voltage, the purple line indicates the duration used for parameter estimation, and the orange line shows the output of our model, which closely follows the real voltage.
The root mean squared error (RMSE) of the model output is given by:
\[ \mathrm{RMSE} = \sqrt{\frac{1}{L}\sum_{t=1}^{L} \left( {V_t-\hat{V}_t}\right)^2} \times 100\% \]
where \(L\) is the total length of the time series, \(V_t\) and \(\hat{V}_t\) are real and model output data respectively. For Cell 1 in Figure 18, RMSE = \(\qty{4.6}{\mV}\) .
The same process is applied on each of the other 5 cells and the results are listed in Table 10. The RMSE of each cell is at the level of about \(\qty{4}{\mV}\) thus the model output is very close to the real voltage. Comparing the parameter estimation results of different cells, the three parameters of cell 2 to 5 are very close to each other, while there is a different between \(\hat{R'}\) of cell 1 and \(\hat{Q}\) and \(\hat{\phi}\) of cell 6 and the other cells.
Discussion
Directly use time-domain current and voltage data recorded during actual battery operation for parameter estimation is challenging due to the BMS’s low sampling frequency and it is usually pretty noisy. In this work, two kind of signals for excitation, multi-sine and PRBS, are used in generating the simulated data, and a comparison of Section 4.2.1 and Section 4.2.2 shows that, consistent with our expectations, estimation using PRBS has better results, which is more tolerant of noise, and maintains good estimation results for data with low sampling frequency.
The three models, generalised Randles model, CPE decomposition model and low frequency model, are incrementally more difficult to estimate when estimating with simulated data in author’s opinion. For CPE decomposition model, the decomposition method yields circuits that often have many RC branches and produce very small R values and very large C values. For low frequency model, the low sampling frequency limits the application of many estimation methods. Vector fitting works well in estimating generalised Randles model, but our simulation is performed under an ideal condition, in practice, to calculate the frequency response, a DFT is needed and thus introduces sampling errors (aliasing), leakage errors (windowing), and frequency errors (limited resolution). In addition, if the current and voltage data in time domain has a DC offset (may not fixed), how to remove this in frequency response is still a problem. Another method to estimate the parameters of the white-box model in time domain, transfer function estimation, or in MATLAB, is more robust than vector fitting in this work. The grey-box model based on the knowledge of EIS and CPE reduces the number of parameters to be estimated for the low frequency model to only 3, as verified in Section 5, yields satisfactory results.
A next direction of work could be to validate the model with more real data, to generalize the “quasi steady state” to a larger range, or to validate the CPE estimation results in the frequency domain.
Conclusion
In conclusion, this project explores the use of time-domain measurements to estimate the parameters of an equivalent circuit model of a lithium cell. 3 ECMs: generalized Randles model, CPE decomposition model and low frequency model are built and several estimation methods are then used. Thereafter, a method for directly estimating the \(R'\) and CPE parameters in a low frequency model using grey-box estimation is developed, which is finally validated using experimental data from a battery module, with a RMSE of around \(\qty{4}{\mV}\).
References
Citation
@online{liu2024,
author = {Liu, Wenyu},
title = {Li-Ion {Cell} {Equivalent} {Circuit} {Model} {Parameters}
{Estimation} from {Time-Domain} {Measurements}},
date = {2024-06-05},
url = {https://wenyuliu.ch/blog/2024/06/Li-Ion-Cell-Equivalent-Circuit-Model-Parameters-Estimation-From-Time-Domain-Measurements//},
langid = {en}
}