
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)12876-9
10.1016/j.heliyon.2024.e36845
e36845
Research Article
A two-step identification approach for an extended nonlinear double-capacitor model
Genario de Oliveira Jose Jr. genario.oliveira@tuwien.ac.at
a⁎
Aras Cisel b
Pallewar Pankaj c
Charkhgard Mohammad c
Sivaraman Thyagesh c
Hametner Christoph a
a Christian Doppler Laboratory for Innovative Control and Monitoring of Automotive Powertrain Systems, Vienna, Austria
b AVL Research and Engineering, Istanbul, Turkey
c AVL List GmbH, Graz, Austria
⁎ Corresponding author. genario.oliveira@tuwien.ac.at
28 8 2024
30 9 2024
28 8 2024
10 18 e368455 10 2023
9 8 2024
22 8 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Equivalent circuits are one of the most used models for Li-ion cells in the automotive area. However, it is a challenge to these models to be able to capture the cell discharge capacity under different loads, while still being accurate on both continuous charge and dynamic tests, fast to compute, and easy to parametrize from non-specialized data. To tackle this challenge, this paper proposes an extension of the nonlinear double capacitor model by increasing its order, parameter dependency with C-rate, and an identification procedure that exploits the pseudo-linear nature of the problem to find the parameter maps. An analogy between the parts of the circuit and the single particle model is also presented to reduce the search space of the identification algorithm and to enhance model interpretability. The performance of the proposed model extension is analyzed and compared to a state-of-the-art model on a challenging LiFePO4 dataset with different characteristics and validated on a realistic drive cycle, obtaining a mean absolute average error of around 20 mV for both training and validation tests.

Keywords

Battery management systems
Energy storage systems
Electrochemical systems
==== Body
pmcNomenclature

A Model system matrix

a Ratio between the surface area and volume of a sphere of radius R m−1

Ae Electrode surface area m2

Aid System matrix with SOC as an input

B Model input matrix

Bid Input matrix with SOC as an input

Ci Capacitance of an RC element F

Cs Surface lithium concentration molm−3

Cid Output matrix of Vs with SOC as an input

Cn Normalized C-rate with respect to an undisclosed capacity

CRCid Output matrix of the RC voltage sum with SOC as an input

CRC Output matrix of the RC voltage sum

cs,max Maximum lithium concentration molm−3

Csi Capacitance of the left subcircuit branch F

Cv Output matrix of Vs

D Diffusivity coefficient in the solid m2s−1

Did Feedthrough term of Vs with SOC as an input

DR0 Feedthrough term of y with SOC as an input

F Faraday's constant Cmol−1

I Applied current A

Ii Left subcircuit branch current A

Iagg Auxiliary current used to aid the identification A

iRCi Resistor current in the RC element A

JLi Lithium molar flux molm−2s−1

Li Linear spline function

pi Poles from the Padé approximation of Xs/Xavg

Q Cell capacity C

Qe Electrode capacity C

R Particle radius m

R0 Model series resistance Ω

Ri Resistance of an RC element Ω

Rsi Resistance of the left subcircuit branch Ω

Rss Model DC gain Ω

s Laplace transform variable

SOC State of charge

T Temperature K

t Time s

uid Original input vector augmented with SOC

Vauxi Auxiliary voltage state V

Vcell Cell terminal voltage V

Vs Voltage which links both circuits through the OCV V

x Model output, same as Vcell

x Model state vector

Xavg Average particle stoichiometry

xid Original state vector without SOC

Xs Surface particle stoichiometry

zi Zeros from the Padé approximation of Xs/Xavg

Greek

δ Electrode thickness m

ϵ Solid phase volume fraction m

λi Ratio between Ri and Rss

λOCV Regularization term on the OCV parameters

Φ Output Regressor matrix with fixed τ, τRCi and λi

ΦOCV Auxiliary matrix to compute the OCV as a matrix multiplication

ΦRss Auxiliary matrix to compute Rss as a matrix multiplication

τ Electrode solid diffusion time constant s

τRCi Time constant associated with the RC elements s

τsi Time constant associated with the left subcircuit branches s

θ Stacked parameter vector with fixed τ, τRCi and λi

θOCV OCV parameter vector

θRss Rss parameter vector

1 Introduction

Lithium-ion batteries play a decisive role in decarbonization strategies for the automotive industry and grid storage applications. This will only increase in importance over time, especially in light of new regulations to ban the selling of new combustion engine vehicles and reduce the emissions of greenhouse gases. Hence, maximizing the efficiency and lifetime of battery and hybrid electric vehicles is crucial to fulfilling these challenges. Understanding different cell models and their limitations in a battery management system (BMS) environment is paramount to ensure a safe and adequate operation of the electric vehicle subsystems, also having a significant impact on energy management strategies, state of power applications, and other critical functionalities in a BMS, such as most state of charge and state of health estimation techniques.

In the context of BMS applications, three goals will be the target for this work: Modeling the fast voltage collapse as the cell becomes discharged (discharge capacity) under different loads; Parameter identification from test data; and fast model computation. Modeling the discharge capacity is vital because if left unchecked, it may lead to unpredicted vehicle shutdown due to the voltage boundaries being violated. The model parameters should be easily identifiable from test data, which is often the case for simple models, but gets increasingly harder as the model complexity increases. Finally, given that the computational resources in a BMS are limited, the model should be easy to compute in an online fashion. While these challenges should be chemistry-independent, a LiFePO4 (LFP) dataset will be used for the methodology proposed in this work.

1.1 Literature review

There is a vast literature on lithium-ion cell modeling. The model considered state-of-the-art is the Doyle-Fuller-Newman (DFN) model, developed in [1]. While it performs extremely well, when it comes to voltage prediction, there are two major drawbacks to using it: a large computational burden and the parameter identification from test data can be extremely challenging and time-consuming. Reference [2] deals with the non-invasive identification of the DFN parameters and compares them to other parameters obtained with cell teardown. Subsequently, a sensitivity analysis of the full parameter set was done, and roughly three weeks were needed to find the model parameters in a parallel computing setup. In [3], another non-invasive identification procedure, exploiting the difference in parameter sensitivity, was developed using a Cuckoo search algorithm. Pre-identification experiments were conducted to determine the open circuit potentials (OCP), stoichiometry values, and the relationship between electrolyte conductivity and concentration. Additionally, the identified model performance yielded overall better results when compared to an invasive parametrization strategy. Another interesting approach to diminish the required time to identify the model parameters was developed in [4], where an LSTM network is used to search the model parameters feasibly, dramatically reducing the identification time to a few hours. A possible way to minimize the model computation time is shown in [5], which is analyzed to solve the model equations and two model reduction strategies.

The simplest electrochemical model commonly used is the Single-Particle Model (SPM). By ignoring the electrolyte and reducing all the electrode dynamics to a single particle, the computational tractability of the model is improved significantly while retaining most of the physical behavior expected from a Lithium-ion cell model, especially at low C-rates. Several works are available on identifying the SPM model, such as [6], where it is shown that the voltage response of the SPM is uniquely defined by six grouped parameters, assuming that the OCP functions are known. This is highly relevant from an identification point of view because it means that it is impossible to fully parametrize the SPM correctly from just voltage and current data without any prior assumptions. In [7], a single particle model was also developed, this time including the electrolyte dynamics, achieving good accuracy, a maximum of 1% voltage error compared against a 313th order computational fluid dynamics (CFD) model while having only 12 states. The fitting was conducted in the frequency domain. [8] modeled a commercial LFP cell with a P2D model coupled with a 2D axisymmetric heat transfer model, combining both optimization and cell teardown.

Parameter identification is often easier on equivalent-circuit models (ECM) approaches; the model is computed significantly faster at the cost of voltage prediction. Despite this, the model quality is often enough for the state of charge estimation, which is one of the primary purposes of a BMS. Probably the most commonly used ECM is the Enhanced-Self-Correcting Model shown in [9] since it captures the battery dynamics well and is simple. One of the shortcomings of this model is that it cannot predict the discharge capacity under different currents well. One step in this direction is the Nonlinear Double-Capacitor (NDC) model, shown in [10], in which a new variable mimicking the surface SOC is used in the Open Circuit Voltage (OCV) curve to model the diffusion effects, focusing on the low-frequency behavior of the model and using an internal resistance and RC elements to model the overpotentials.

In order to improve the model performance on ECMs, often the parameters are not constant but a function of SOC and sometimes C-rate. In this case, the parameter identification of these models can become quite challenging due to overfitting and generalization issues. [11] is an excellent example of how to conduct this parameter identification; linear spline functions represent the functional form of the parameter dependencies, the model parameters are identified using a genetic algorithm, and the validation is done for an LFP dataset. Another work conducted on an LFP dataset for different temperatures and SOCs was the characterization in the frequency domain from EIS data shown in [12], where it is seen how important it is to have a good initial guess for the parameters when the number of parameters of the model is large, which is the case when a functional dependence of the parameters with SOC and temperature is present. One assumption commonly made in literature is that the OCV curve is available. This is not always the case when identifying the parameter from experimental data. The local identifiability of an RC model with constant parameters and without prior information on the OCV curve is analyzed in [13]. The model parameters and the OCV-SOC relationship are identified using a nonlinear least squares algorithm, showing good agreement with the experimental discharge data. In [14], a simpler version of the model proposed here with non C-rate dependent parameters is coupled with a modified Preisach model in order to model the hysteresis and a nonlinear observer is developed, being an improvement over classical EKF and UKF approaches for this class of model. A comparative study of different equivalent circuit models on different chemistries, namely LFP, NMC, LMO, and NCA, is conducted in [15]. The models considered are the classic 1 or 2 RC models linked with a hysteresis term when necessary. They find that these ECMs perform better with dynamic current profiles but have worse performances on loads with lower dynamics, such as constant charge/discharge pulses. These non-dynamic loads are used in this work and are one of the reasons why C-rate dependent parameters were investigated.

A very similar approach to the one proposed in [10] and extended here is presented in [16], where the diffusion in the solid is also approximated inspired on the SPM model by considering an average SOC and a surface SOC. The major difference between these approaches is how these surface and average soc variables show up in the output equation of the models. An identification procedure is also developed showing good agreement with the test data.

The research gap identified by us, which we tackle in this paper, is finding an equivalent circuit model that is at the same time fast to compute, easy to identify from non-specialized data, and can correctly model the discharge capacities under different C-rates. One example of model extension is done in [17], where a hybrid data-driven approach is done. The accuracy of the model developed in [10] is improved by using training a neural network on the outputs of the original model, improving the accuracy over different C-rates.

1.2 Contribution

This work aims to achieve several objectives: firstly, modeling the voltage response under various constant C-rates, which is a limitation of commonly used ECMs; secondly, implementing a parameter identification strategy for the proposed modified model. Throughout these tasks, a focus is to ensure a swift computation time for the model. To address these goals, we introduce an enhanced version of the NDC model initially developed in [10]. This extended model incorporates C-rate dependent parameters, allowing for an accurate representation of discharge capacity under diverse loads. We establish a novel analogy between the electrode transfer function of the SPM and the extended model, resulting in a reduction of the number of parameters requiring identification. Additionally, we propose a two-step identification approach to concurrently determine the model parameters and the open circuit voltage. It will be shown that the proposed extended approach fulfills the research gap identified by us. The proposed model extension and identification approach are validated on an experimental LFP dataset, where the model accuracy, both in terms of output voltage and discharge capacity at different C-rates, is compared against a state-of-the-art approach model. A comparison of the identification times is also done. Notably, the ECM structure of the model ensures efficient online computation.

1.3 Organization

This work is organized as follows: Section 2 details the model structure and briefly explains the cell and experimental datasets. Section 3 describes the identification procedure and validation steps. Section 4 discusses the results obtained for the dataset presented in 2.

2 Extended equivalent circuit model

The ECM structure is developed in this section, together with a method to reduce the number of effective parameters of part of the circuit, based on an analogy between the transfer function of a part of the model and the SPM solid diffusion PDE in an electrode. A time-varying state-space model representation is also derived. The dataset that is used for parameter identification and validation is also presented at the end of this section.

2.1 Model structure

The equivalent circuit model considered here is shown in Fig. 1, an extended version of the nonlinear double capacitor model, shown in [10], with more parallel RC branches in the left circuit and current assumed positive on charge. In this model, the link between both circuits is the voltage signal Vs, which is interpreted here as a surface SOC.Figure 1 Generic nonlinear multiple capacitor circuit.

Figure 1

The derivation of a transfer function between I and Vs is done by equating the Laplace transform of the voltage drops in the left circuit for the parallel branches, so for an arbitrary branch k, the following equations are considered:(1) (Rs0+1sCs0)I0=(Rsk+1sCsk)Ik

(2) Ik=CskCs0sRs0Cs0+1sRskCsk+1.

Where (1) comes from the fact that the voltage drop in the parallel branches is the same. Using Kirchhoff's Current Law and the capacitor's constitutive voltage relation in the Laplace domain, the following is obtained:(3) I(s)=I0(s)+I1(s)+...+In(s),

(4) Vs(s)=I0(s)sCs0.

Combining Eqns. (2), (3) and (4), one obtains:(5) I(s)sCs0Vs(s)=(1+Cs1Cs0sRs0Cs0+1sRs1Cs1+1+...CsmCs0sRs0Cs0+1sRsmCsm+1).

The capacity is assumed to be known and the SOC is defined as:(6) SOC(s)=I(s)sQ,

where Q is the cell capacity. By substituting Eqn. (6) into Eqn. (5), the following equation is obtained:(7) Vs(s)SOC(s)=1Cs0sRs0Cs0+1+Cs1sRs1Cs1+1+...CsmsRsmCsm+1×QsRs0Cs0+1.

The transfer function above has a DC gain of one if the sum of the capacitor values equals the capacity of the cell, which is assumed here.

Additionally, in the next sections, it will be shown how to choose the values for the k-th elements Rsk and Csk of the proposed circuit in order to impose similar dynamics between Vs(s) and SOC(s), and an electrode average and surface stoichiometries in the context of the SPM.

2.2 Single particle model diffusion transfer function

Here the transfer function between the electrode average and surface stoichiometry will be derived from the SPM equations for the diffusion in the electrode. The transfer function between the surface lithium concentration Cs and the lithium flux JLi(s), taken from [18] and already assuming spherical diffusion in the solid particle, is written below:(8) Cs(s)JLi(s)=RDe2RsD−11+RsD+e2RsD(RsD−1),

where R is the particle radius and D is the diffusivity coefficient in the solid. Further simplification of the equation above is done by assuming that both the electrodes can be well approximated by one spherical particle, spatially independent variables, and uniform distribution of reactions throughout the cell thickness. This results in:(9) JLi(s)=I(s)aδFAe,

(10) a=3ϵR,

where a is the ratio between the surface area and volume of a sphere of radius R, ϵ is the solid phase volume fraction, δ the electrode thickness, F is the Faraday constant and Ae is the electrode surface area. Additionally, the electrode stoichiometries are defined as:(11) Xs(s)=Cs(s)cs,max,

(12) Xavg(s)=I(s)sQe,

where Xs(s) is the electrode surface stoichiometry, Xavg(s) the average electrode stoichiometry, cs,max the electrode's maximum lithium concentration and the electrode capacity Qe=cs,maxϵδFAe. Substituting the electrode diffusion time constant τ=R2D, as in [6], and Eqns. (9)-(12) into Eqn. (8), one obtains:(13) Xs(s)Xavg(s)=τ3s(e2sτ−1)1+sτ+e2sτ(sτ−1).

The transfer function between Xs(s) and Xavg(s), shown in Eqn. (13), is fully parametrized by one lumped parameter, the diffusion time constant τ.

2.2.1 Transfer function matching

The main idea behind the left circuit of the NDC-like approaches is to simplify the solid electrode diffusion effects from both electrodes, as a bulk SOC and a surface SOC, which appears in the model output through the OCV curve, instead of modeling each electrode individually, using the OCP functions together with the surface stoichiometries for each electrode. Because of this, it is proposed here to use Eqn. (13) as the basis for identifying the parameters from Eqn. (7) by approximating the stoichiometries of an average electrode Xs(s) and Xavg(s) by the surface SOC, represented by Vs(s), and the bulk SOC, simply referred to as SOC(s).

In order to match the transfer functions defined in Eqns. (7) and (13), a Padé approximation, as in [18], will be done in Eqn. (13). The resulting transfer functions for different orders are shown in Table 1. For m RC branches in the circuit shown in Fig. 1, there are 2m+2 parameters that need to be determined by comparing the coefficients with the Padé approximant of order m of Eqn. (13) a total of 2m+1 equations are obtained; thus, without loss of generality, the following condition is imposed to get a unique solution:(14) Rs0=0.

An example of this is shown next in more detail.Table 1 Padé approximations of Xs(s)Xavg(s).

Table 1Order	Transfer Function	
1	2τ21s+1τ35s+1	
2	τ2495s2+4τ35s+1τ23465s2+3τ55s+1	
3	4τ3225225s3+2τ2585s2+2τ15s+1τ3675675s3+2τ22275s2+τ15s+1	

2.2.2 RC parameter determination for m=2

Introducing the new variable τsk=RskCsk, Eqn. (7) becomes:(15) Vs(s)SOC(s)=1Cs0sτs0+1+Cs1sτs1+1+Cs2sτs2+1Qsτs0+1,

expanding Eqn. (15) and comparing the coefficients with the second order Padé approximant for the numerator:(16) τs1τs2=τ2/495,

(17) τs1+τs2=4τ/35,

and for the denominator:(18) Cs0τs1τs2+Cs1τs0τs2+Cs2τs0τs1=Qτ23465,

(19) Cs0(τs1+τs2)+Cs1(τs0+τs2)......+Cs2(τs0+τs1)=3Qτ55,

(20) Cs0+Cs1+Cs2=Q.

The time constants τs1 and τs1 are determined by Eqns. (16)-(17). Once these two time constants are determined, four unknowns remain, namely Cs0, Cs1, Cs2 and τs0 and three Eqns. (18)-(20). These need to be augmented by Eqn. (14), or τs0=0, in order to fully determine the parameters, resulting in a linear system to be solved.

2.3 State-space model

The state-space representation for the model depicted in Fig. 1, with orders m=2 and n=2 will be parametrized here. The state vector for this model is defined as:(21) x=[iRC1iRC2Vsaux1Vsaux2SOC]T,

where iRC1 and iRC2 are the resistor currents in the right RC branches, SOC is the state of charge, Vsaux1 and Vsaux2 are auxiliary states not directly related to circuit elements and SOC is the state of charge, defined as:(22) SOC˙=1QI(t),

where Q is the cell capacity in A.s, I(t) the input current in A, which is also the system input and assumed positive on charge. In order to implement the second-order Padé approximation for the left NDC subcircuit, written as:(23) Vs(s)SOC(s)=τ2495s2+4τ35s+1τ23465s2+3τ55s+1,

the poles and zeros of this transfer function are explicitly computed as(24) [p1(t),p2(t)]=12τ(t)(−189±21861),

(25) [z1(t),z2(t)]=17τ(t)(−198±14949).

The second-order approximation, taken from Table 1 and shown in Eqn. (23), was chosen here in order to have some trade-off between complexity and accuracy, achieving adequate results. Other orders might be considered as well e.g. taking into account the frequency band in which the approximation should be valid to select the order. Note that in order to capture the C-rate dependent effects, the parameter τ is considered time-varying and a function of current. More details will be given in the next section. Defining the matrices A(t) and B as(26) A(t)=[−1τRC100000−1τRC200000−1p10−1p1001p2(1−z1z2)−1p21p2z1z200000],

(27) B=[1τRC11τRC2001Q],

where the variables τRC1 and τRC2 are the time constants associated with the RC elements in the left circuit and are defined as τRCi=RiCi. To avoid a notation overload, the variables p1(t), p2(t), z1(t) and z2(t) are written without the time-dependency. Finally, the state update equation is written as(28) x˙=A(t)x+Bu,

where the state vector is defined in Eqn. (21). The matrices A and B are defined respectively in Eqns. (26) and (27). The SOC state dynamics, defined in Eqn. (22), is also included here. The system output y is the cell voltage Vcell and it depends on the auxiliary output Vs. This auxiliary output plays the role similar to a surface SOC, and is computed from the state vector as(29) Vs=Cv(t)x,

(30) Cv(t)=[00p1p2(1−z1z2)(1−p1p2)p1p2z1z2].

The output voltage from the RC elements, uRC is computed using another auxiliary matrix(31) CRC(t)=[R1(t)R2(t)000].

Finally, the output equation for the state-space representation, Vcell, is written as(32) y=OCV(Cv(t)x)+CRC(t)x+R0(t)u.

Substituting Eqns. (29)-(31) into Eqn. (32)(33) y=OCV(Cv(t)x)+CRC(t)x+R0(t)u.

The resistances R0(t), R1(t), and R2(t) are also considered to be time-varying, depending on the current. The discrete-time equivalent of the system described here was simulated using a backward discretization scheme.

2.4 Parameter scheduling

The state-space model developed in the previous section depends on the input current. A possible physical explanation for τ is that the assumption of spherical diffusion, in Eq. (13), does not hold entirely for the LFP cells tested, with the τ here being used more as an effective parameter. Additionally, there is an assumption that the lithium transport happens entirely because of diffusion, ignoring the role of additional transport phenomena such as migration.

Modeling the direct resistance R0 and RC resistances R1, R2 as functions of the current also significantly improves the model fit. This is motivated by the nonlinear dependency of the overpotentials, which are linked with the resistance terms, with the current. Given that the overpotentials also depend on the surface stoichiometry and remembering the analogy made between Vs and Xs, it would also be feasible to consider a second dependency of the resistances R0, R1 and R2 with Vs. This was not done independently for each variable here to avoid additional complexity and simplify the identification problem, but will be considered later. This also applies to the RC time constants. Since the OCV is not assumed to be known, it will also be identified jointly with the data, and all the functional dependencies will be represented as linear splines, as done in [19], for dimensions one, and as multiplicative linear splines when the dimension is two, resulting in(34) f(w)=a0+∑i=1nwaiLi(w,wi),

(35) g(w,d)=a0+∑i=1nw∑j=1ndaijLi(w,wi)Li(d,di),

where ai and aij represent the weights defining the linear splines and the breakpoint vectors wi and di. The function Li(.) is defined as(36) Li(w,wi)={w−wi,if w≥wi0,otherwise.

2.5 The dataset

Test data from a LiFePo4 cell was used to identify the model parameters and validate the results. There are four different testing profiles that are used for identification, and one test profile is used for validation. Due to confidentiality reasons, we are not able to disclose some key cell-specific or test-specific details, such as nominal capacities, voltage levels, C-rates, and test durations. Thus, when applicable, these values are shown normalized or without the label ticks in the figures. The first testing profile, shown later, is a pseudo-OCV test at a low C-rate. The second test, shown partially in Fig. 2, is a series of continuous charging and discharging profiles at different C-rates. The remaining two tests are more dynamic, with shorter pulses and relaxation profiles at higher C-rates. The normalized current profile, SOC channel, and voltage are also shown partially for these tests in Figure 3, Figure 4. The OCV is assumed to be unknown and will be identified jointly with the model parameters and the SOC is obtained with coulomb counting. Fig. 5 shows the validation profile, which is a drive cycle and have a current profile with levels different than the training set, passing as previously unseen data, in which the generalization capabilities of the model will be explored.Figure 2 Quasi-static Charging Test.

Figure 2

Figure 3 Dynamic Test I.

Figure 3

Figure 4 Dynamic Test II.

Figure 4

Figure 5 Drive Cycle Data.

Figure 5

3 Model identification procedure

In this section, the model identification procedure to find the model parameters from experimental data will be outlined. There are three steps that are fundamental to ensure the success of this approach, namely initialization, optimization, and model refinement.

Each of these steps will be explained in detail in the next subsections. Since the SOC, and the capacity, are assumed to be known for the parameter identification, the system represented by Eqns. (28) and (33) are rearranged to eliminate the SOC dynamics given by the Coulomb Counting equation and introduce SOC as an input, yielding(37) x˙id=Aid(t)xid+Biduid,

(38) y=OCV(Cid(t)xid+Did(t)uid)+...CRCid(t)xid+DR0(t)uid,

with uid=[ISOC]T, xid=[iRC1iRC2Vsaux1Vsaux2]T and the matrices from the equations above:(39) Aid=[−1τRC10000−1τRC20000−1p10001p2(1−z1z2)−1p2],

(40) Bid=[1τRC101τRC200−1p101p2z1z2],

(41) Cid=[00p1p2(1−z1z2)(1−p1p2)],

(42) CRCid=[R1R200],

(43) Did=[0p1p2z1z2],

(44) DR0=[R00].

The Eqns. (37)-(44) denotes the state-space system which is used for parameter identification.

3.1 OCV initial guess

The OCV relationship is represented with a linear spline as in Eqn. (34). The parameters and breakpoints are found by approximating the OCV curve extracted from a pseudo-OCV test, as in Fig. 6. The breakpoints for the OCV are then fixed and re-optimized only during the last stage. While this OCV initialization gives a pretty good estimation of the OCV, some areas need to be estimated jointly with the other cell parameters, given its dependency of Vs in the model, instead of SOC. Additionally, all the test data was scanned for resting times of more than ten minutes and they are used as a constraint for the OCV curve, with the OCV needing to be smaller than a rest from charge and greater than a rest from discharge at any specific SOC, given that no hysteresis representation is being considered here. These points are shown in Fig. 7. This procedure makes it harder for the identification algorithm to use the OCV curve as an effective parameter when Vs is between 0 and 1, given that the constraints are quite tight in some areas.Figure 6 Pseudo-OCV data.

Figure 6

Figure 7 Rest points extracted from data.

Figure 7

3.2 Two-step parameter optimization

The identification of the model parameters is made in two different steps. A global optimization, which optimizes the parameters τRC1, τRC2, the functional relationship τ(I) and the parameters λ1 and λ2, which will be defined later, and a local optimization which optimizes the resistance maps and the OCV. This is done because once the states are fully determined, the reduced local problem becomes linear in the parameters and is solved with a quadratic programming procedure.

Using a matrix notation for Eqn. (34), and using the variable θOCV for the parameter vector, one can rewrite it as(45) OCV=ΦOCV(Vs)θOCV.

Similarly, the same pseudo-linear structure for the OCV, shown in (45) and in Eqns. (34) and (36), is present for the resistance terms. However, to decrease overfitting issues and improve the model validation, a small modification was done introducing the following auxiliary variables λ1, λ2 and the steady-state gain Rss(I,Vs), which is represented as in Eqn. (35).(46) R0=(1−λ1−λ2)Rss,

(47) R1=λ1Rss,

(48) R2=λ2Rss,

(49) 0≤λ1+λ2≤1,

(50) λ1,2≥0

where a fixed ratio between the resistances is used in order to simplify the identification problem and to avoid identifying three different parameter maps instead of one, the relationship between R0, R1, R2 and Rss is shown in Eqns. (46)-(50). While the price to be paid is that this correlates R0, R1, and R2, identifying three parameter maps at once is quite challenging and the steady-state gain of the model is preserved.

Writing Rss also in matrix form, one obtains(51) Rss=ΦRss(Vs,I)θRss.

For the second step, the output Eqn. (38) is written as a function of the states in matrix form as(52) y=Φ(Vs,I)θ,

(53) Φ(Vs,I)=[ΦOCV(Vs)ΦRss(Vs,I)Iagg],

(54) θ=[θOCVθRss],

(55) Iagg=(1−λ1−λ2)I+λ1iRC1+λ2iRC2.

The cell terminal voltage in the model is now in a format suitable for estimation, as written in Eqn. (52). The definition of the variables Φ(Vs,I), θ and Iagg shown respectively in Eqns. (53)-(55). The identification algorithm proposed is then to minimize a loss function V defined as:(56) minτ(I),τRCi,λi⁡V(τ(I),τRCi,λi)s.t.τ(I)≥0,τRCi≥0,λi≥0,∑λi≤1,

where the loss function is defined as the solution for the convex subproblem(57) V(τ(I),τRCi,λi)=minθ⁡SSEw(y−Φ(Vs,I)θ)+...λOCV||dOCV(θ)dVs||2s.t.dOCV(θ)dVs≥0,

λOCV is a variable that penalizes the non-smoothness of the identified OCV map and was determined heuristically to avoid overfitting the OCV curve. SSEw is the sum of squared errors of the model using a weighting vector that gives a higher weight to the points close to the discharge limit and the charge limit. The problem defined in Eqn. (57) is convex, hence, the problem can be easily rewritten into a quadratic program for a given set of (τ(I),τRCi,λi), guaranteeing local optimally for the parameter vector θ, which is found by using the command ‘quadprog’ in MATLAB. The convergence of the proposed approach depends on the convergence of the global problem defined in Eqn. (56), which is sensitive to the parameter initialization, which is discussed next.

3.2.1 Parameter initialization

In order to solve the problem posed in Eqn. (56), a divide-and-conquer approach is used. First, only the continuous charge and discharge profiles at different C-rates are used for the identification of a reduced version of the model, obtained by not considering the RC elements in the main circuit, thus, eliminating two states and simplifying the problem. This can be done because, under a constant charge or discharge for a long period of time, the RC elements are well approximated by a steady-state gain Rss as a function of the current. No dependency of Rss with Vs is considered at this stage. After that, the identified values of τ(I) and Rss(I) are used to initialize the full problem with all the data sets. The parameter initialization strategy used here was tailored to the experimental dataset available and should be modified to consider the specific characteristics of each different testing data. Despite that, if the dataset has similar characteristics, with charge/discharge data and some dynamic tests, the modifications required should be minimal.

3.3 Model refinement

Finally, the last step of the identification procedure is to optimize the breakpoints for the current and Vs that are an integral part of finding the breakpoints that define the linear spline functions used for the OCV and Rss. Two different steps are considered; the first is to delete the breakpoints that have less impact on the loss function defined in Eq. (57) until the model quality is still considered acceptable by the user. This is shown for the OCV in Fig. 8. The second step is a breakpoint addition phase, where a new breakpoint is added iteratively, followed by a gradient-based optimization of Eq. (57) where the model parameters are fixed and all breakpoints are optimized at once. This procedure of breakpoint addition is shown in Fig. 9. At the end of the addition step, the number of SOC breakpoints was fixed to 21, where the loss function, normalized by the smallest value, stops decreasing significantly in Fig. 9.Figure 8 Linear spline breakpoint removal step.

Figure 8

Figure 9 Linear spline breakpoint addition step.

Figure 9

3.4 Model limitations

Given that the proposed model extension and identification step is innately data-driven, there are some key limitations that this approach is subject to. The biggest one is that it is hard for this type of model to extrapolate outside the test settings well, e.g. If the test C-rates are 1 and 3, extrapolating the model parameters to simulate the behavior on a C-rate of 5 is unreliable, on the other hand, simulating for a C-rate of 2 usually would be fine. Directly related to this, is the fact that if the model is required to be accurate on diverse C-rates that are key for the model purpose, it is necessary that the testing data used for the parameter identification cover these regions.

4 Results and discussion

After the model parameters are identified using the procedure described in the previous section, the performance on the different training datasets, together with the drive cycle validation and said parameters, are displayed here. Whenever relevant, the results will be compared versus a state-of-the-art model, RC2, which consists of an OCV curve, an internal resistance, and two RC elements. All the resistances for the RC2 model are SOC-dependent and the OCV curves of both models are the same.

4.1 Quasi-static charging tests

The simulation of the dataset with continuous constant discharge and charge at different C-rates is shown in Fig. 10, with also the rest intervals and quasi-rest voltages. Fig. 11 shows the data and model voltage response on the continuous discharge regions versus the depth of discharge in Ah, together with the R2 values for each curve. It is evident from these figures that the model captures the correct voltage behavior, especially at the extreme points, with some deviations being present during the charge/discharge processes depending on the C-rate. The relative error between the predicted discharge capacities by the models and the one computed at different C-rates by coulomb counting is depicted in Table 2. These errors are minimal for the proposed approach and indicate that the proposed model achieves the initial goal of modeling the fast voltage drops that happen at different C-rates. The errors for the RC2 model are bigger, but might be acceptable depending on the intended application. A smaller error on the discharge capacity might be achieved for the RC2 model if an additional C-rate dependency on its parameters is considered. However, since there is already a SOC dependency, the challenge in identifying a joint SOC and C-rate dependency is much greater than any of them individually. The very slow diffusion dynamics, which are visible in the long rest periods in Fig. 10 are slightly under-modeled, especially when combined with low C-rates. This was expected to be the case since the beginning, since both electrode solid diffusion dynamics are approximated by an equivalent electrode transfer function. This can be seen from the different diffusion dynamics at higher and lower SOC values, each affected more by one of the electrodes. In the case of LFP, this is not a bad approximation, since the LFP electrode is much more nonlinear in the cell operating region than the graphite one. More importantly, and also depicted in the identified OCV shown in Fig. 14, the rest points match the OCV with a small error, so even if the exact dynamics are not correct, the model relaxes approximately to the correct rest voltage. The OCV constraints, shown in the same figure, play a key role in enforcing that this happens. It is also crucial that the initial guess of the OCV meets these constraints, so the optimizer can start from a feasible point. Additionally, the last steps during the model refinement, where only the required breakpoints are used, play a key role in keeping the number of model parameters in check. The identified dependency of both Rss and τ with C-rate are shown in Figure 12, Figure 13. These figures show that the identified parameters change significantly with low C-rates and do not change much otherwise. The NDC time constants have a more pronounced impact and are more identifiable from the lower C-rate values, since for the high C-rates the continuous charge tests are not present. The time constants of the two RC elements were identified as τRC1=250s and τRC2=2600s.Figure 10 Quasi-static charging data and prediction.

Figure 10

Figure 11 Discharge at different C-rates.

Figure 11

Table 2 Discharge Capacity Relative Error.

Table 2Model	1Cn	5Cn	10Cn	17Cn	25Cn	40Cn	
2NDC-RC2	0.16%	0.12%	0.05%	-0.01%	-0.01%	-0.05%	
RC2	5.53%	4.83%	5.04%	4.27%	4.67%	8.47%	

Figure 12 Identified Rss, λ1, λ2 = 0.164.

Figure 12

Figure 13 Identified τ.

Figure 13

Figure 14 Identified Open-Circuit Voltage.

Figure 14

4.2 Dynamic tests

The two remaining dynamic tests in the training set are shown in Figure 15, Figure 16 together with the model performance on these tests, with an additional zoom in the areas with faster dynamics, showing that the voltage response in these areas, which is directly linked with the RC elements, matches the dataset well. These are the tests where the proposed modifications have the least impact, due to the short duration of the current pulses, compared with the other tests, the concentration gradient between surface and bulk doesn't increase too much, which means that the impact of the left subcircuit is low in the voltage response. Hence, the standard RC2 model is sufficient to model these without any problem.Figure 15 Dynamic data I and prediction.

Figure 15

Figure 16 Dynamic data II and prediction.

Figure 16

4.3 Validation and performance

The validation on the cycle test, depicted in Fig. 17, is highly dynamic. While the current levels are exciting only a region of the current-dependent parameters, they are not explicitly contained in the training set, which displays the interpolating model behavior of these parameters. There is a slight offset present in the step responses in this dataset which might be due to the unmodeled hysteresis behavior in this LFP cell. On the other hand, the effectiveness of the joint OCV identification is shown, since the rest voltages predicted by the model match the drive cycle rest periods quite well. Table 3 displays the mean absolute average error in mV, considering a sampling time of 1s, for the entire dataset for both models. The error level of the validation profile is of the same magnitude as the training set, indicating that the training procedure was successful and the data was not overfitted, at least in the region excited by the drive cycle. It is evident by looking at Table 3 that the inclusion of the subcircuit does not improve significantly the performance of the model on more dynamic cycles, or even in the validation drive cycle. However, the inclusion of the NDC elements greatly improves the model accuracy on the Quasi-static charging data and Pseudo-OCV tests. So there is a complexity trade-off that should be considered when choosing between these models for a real application. The last model refinement steps are also done in a fashion that eliminates redundant degrees of freedom from the model. Additionally, an extra model refinement step was performed at the end where the transfer function poles and zeros are not constrained anymore by Eqns. (25) and (24). This did not result in any significant improvement in the model quality. The whole identification procedure was done in MATLAB and took around 3 hours on a laptop with an i7-1260 processor, 32 GB of RAM, and SSD. As is the case for most equivalent-circuit models, the proposed model is extremely fast to be computed and definitely real-time capable, simulating around 250 h of cell data in roughly 0.8 seconds, while the RC2 model would take 0.6 seconds, using the same computer. Considering a 1 s sampling time, the model has a real-time factor of 1125000. This means that it can be used online without any problem and also for demanding optimizations.Figure 17 Drive cycle data and prediction.

Figure 17

Table 3 Mean Absolute Average Error in mV for each test.

Table 3Model	Pseudo-OCV	Quasi-static	Dyn. I	Dyn. II	Drive Cycle	
2NDC-RC2	11.9	14.7	24.7	21.1	19.8	
RC2	34.7	34.1	25.1	22.0	20.5	

5 Conclusion

The goal of this paper, which reflects our identified research gaps, was: to obtain an equivalent circuit model that could predict the cell discharge capacity under different C-rates; model simultaneously the dynamic tests and continuous charge tests; a straightforward parameter identification strategy from non-specialized data; and a fast computation time.

The proposed approach extends the model presented in [10] by introducing C-rate dependent parameters, identifying the time constants associated with the left subcircuit in a physics-motivated comparison with the SPM, which constrains the degrees of freedom from the parameter optimization. An identification strategy was presented to find the model parameter maps jointly with the OCV.

This was validated by comparing the performance of the proposed method with an OCV-R0-RC2 model, which is a state-of-the-art model in several applications, on a challenging LFP dataset.

The extended model performance on continuous load tests and prediction of the discharge capacity over different C-rates is greatly improved over the baseline. On the other hand, it did not improve the voltage error metrics significantly for the more dynamic tests, where the impact of the diffusion effects is reduced in comparison to the more static tests. By incorporating C-rate dependent parameters, it predicts discharge capacity and fast voltage collapse, particularly in low SOC regions, crucial for avoiding unexpected vehicle shutdowns in a BMS. These results are achieved despite known limitations, such as the absence of a hysteresis model and a good solid diffusion approximation for both high and low SOC values at the same time. All of this is achieved while still maintaining fast computation times and a straightforward parameter identification strategy.

Hence, if only performance in the more dynamic cycles is important for the model application, the extra complexity of the proposed model does not pay off. However, if estimating the SoE at low SOC values, predicting the discharge capacities, and optimizing the continuous charging modeling are relevant for a given application, there is a significant performance increase by adopting the model extension shown here.

Possible research directions to extend this work might be exploring distinct dynamics for each electrode, which would involve the identification of both open circuit potentials or refining the solid diffusion model approximation based on the SPM.

CRediT authorship contribution statement

Jose Genario de Oliveira: Writing – review & editing, Writing – original draft, Visualization, Validation, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Cisel Aras: Writing – review & editing, Visualization, Supervision, Resources, Data curation. Pankaj Pallewar: Project administration, Investigation. Mohammad Charkhgard: Writing – review & editing, Validation, Supervision, Investigation. Thyagesh Sivaraman: Resources, Project administration, Funding acquisition. Christoph Hametner: Writing – review & editing, Visualization, Validation, Supervision, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement

The financial support by the Austrian 10.13039/501100012416 Federal Ministry for Digital and Economic Affairs ; the 10.13039/100010132 National Foundation for Research, Technology and Development ; and the 10.13039/501100006012 Christian Doppler Research Association are gratefully acknowledged.
==== Refs
References

1 Doyle M. Fuller T.F. Newman J. Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell J. Electrochem. Soc. 140 6 1993 1526
2 Khalik Z. Donkers M. Sturm J. Bergveld H. Parameter estimation of the Doyle–Fuller–Newman model for lithium-ion batteries by parameter normalization, grouping, and sensitivity analysis J. Power Sources 499 2021 229901
3 Li W. Demir I. Cao D. Jöst D. Ringbeck F. Junker M. Sauer D.U. Data-driven systematic parameter identification of an electrochemical model for lithium-ion batteries with artificial intelligence Energy Storage Mater. 44 2022 557 570
4 Xu L. Lin X. Xie Y. Hu X. Enabling high-fidelity electrochemical p2d modeling of lithium-ion batteries via fast and non-destructive parameter identification Energy Storage Mater. 45 2022 952 968
5 Xia L. Najafi E. Bergveld H. Donkers M. A computationally efficient implementation of an electrochemistry-based model for lithium-ion batteries 20th IFAC World Congress IFAC-PapersOnLine 50 1 2017 2169 2174
6 Bizeray A.M. Kim J. Duncan S.R. Howey D.A. Identifiability and parameter estimation of the single particle lithium-ion battery model IEEE Trans. Control Syst. Technol. 27 5 2019 1862 1877
7 Smith K.A. Rahn C.D. Wang C.-Y. Control oriented 1d electrochemical model of lithium ion battery Energy Convers. Manag. 48 9 2007 2565 2578
8 Muñoz P. Humana R. Falagüerra T. Correa G. Parameter optimization of an electrochemical and thermal model for a lithium-ion commercial battery J. Energy Storage 32 2020 101803
9 Plett G. Battery Management Systems, Volume I: Battery Modeling 2015 Artech House
10 Tian N. Fang H. Chen J. Wang Y. Nonlinear double-capacitor model for rechargeable batteries: modeling, identification, and validation IEEE Trans. Control Syst. Technol. 29 1 2021 370 384
11 Hu Y. Yurkovich S. Guezennec Y. Yurkovich B. Electro-thermal battery model identification for automotive applications J. Power Sources 196 1 2011 449 457
12 Birkl C. Howey D. Model identification and parameter estimation for lifepo4 batteries IET Hybrid and Electric Vehicles Conference 2013 (HEVC 2013) 2013 1 6
13 Tian N. Wang Y. Chen J. Fang H. On parameter identification of an equivalent circuit model for lithium-ion batteries 2017 IEEE Conference on Control Technology and Applications (CCTA) 2017 187 192
14 Movahedi H. Tian N. Fang H. Rajamani R. Hysteresis compensation and nonlinear observer design for state-of-charge estimation using a nonlinear double-capacitor li-ion battery model IEEE/ASME Trans. Mechatron. 27 1 2022 594 604
15 Li C. Cui N. Cui Z. Wang C. Zhang C. Novel equivalent circuit model for high-energy lithium-ion batteries considering the effect of nonlinear solid-phase diffusion J. Power Sources 523 2022 230993
16 Tran M.-K. DaCosta A. Mevawalla A. Panchal S. Fowler M. Comparative study of equivalent circuit models performance in four common lithium-ion batteries: lfp, nmc, lmo, nca Batteries 7 3 2021
17 Tu H. Moura S. Wang Y. Fang H. Integrating physics-based modeling with machine learning for lithium-ion batteries Appl. Energy 329 2023 120289
18 Forman J.C. Bashash S. Stein J.L. Fathy H.K. Reduction of an electrochemistry-based li-ion battery model via quasi-linearization and Padé approximation J. Electrochem. Soc. 158 2011
19 Hu Y. Yurkovich S. Linear parameter varying battery model identification using subspace methods J. Power Sources 196 5 2011 2913 2923
