Accelerating parameter estimation in Doyle–Fuller–Newman model for lithium-ion batteries

Abstract

Purpose

This paper aims to solve the parameter identification problem to estimate the parameters in electrochemical models of the lithium-ion battery.

Design/methodology/approach

The parameter estimation framework is applied to the Doyle-Fuller-Newman (DFN) model containing a total of 44 parameters. The DFN model is fit to experimental data obtained through the cycling of Li-ion cells. The parameter estimation is performed by minimizing the least-squares difference between the experimentally measured and numerically computed voltage curves. The minimization is performed using a state-of-the-art hybrid minimization algorithm.

Findings

The DFN model parameter estimation is performed within 14 h, which is a significant improvement over previous works. The mean absolute error for the converged parameters is less than 7 mV.

Originality/value

To the best of the authors’ knowledge, application of a hybrid optimization framework is new in the field of electrical modelling of lithium-ion cells. This approach saves much time in parameterization of models with a high number of parameters while achieving a high-quality fit.

Citation

Reddy, S.R., Scharrer, M.K., Pichler, F., Watzenig, D. and Dulikravich, G.S. (2019), "Accelerating parameter estimation in Doyle–Fuller–Newman model for lithium-ion batteries", COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, Vol. 38 No. 5, pp. 1533-1544. https://doi.org/10.1108/COMPEL-12-2018-0533

1. Introduction

The applications of lithium-ion batteries have drastically increased over the past decade. With the continuous implementation of the Li-ion (Li⁺) cells in household appliances, automotive, aerospace and defense industries, accurate modeling and simulation of them is paramount. Accurate analysis of the battery can sometimes require the internal state of the cell to be known. This internal state can include abstract quantities (e.g. state of charge [SoC] and state of health, i.e. the usable capacity, power capabilities among others) and physical quantities (e.g. potentials and concentrations). Some of these quantities can be measured through experimentation. In several cases, the material properties of the cell can also be of interest. These material properties sometimes cannot be measured directly and must be estimated, often non-intrusively. This gives rise to the traditional parameter estimation problem.

Parameter estimation techniques attempt to identify certain parameters in a model using only the model response. Parameter estimation techniques can be non-intrusive and non-destructive depending on whether the model response can be obtained non-intrusively and non-destructively. The parameter estimation problem in this work can be stated as follows: Given only the voltage, how can the material properties and model parameters of the lithium-ion cell model be estimated?

A significant effort has been dedicated to solving the parameter estimation problem in Li⁺ batteries. This literature review predominantly focuses on studies that perform the parameter estimation offline. Schmidt et al. (2010) successfully identified 33 parameters in their electrochemical model using a pattern search algorithm. They also used the Fisher information to determine the identifiable parameters. Speltino et al. (2009) performed parameter identification in a single-particle-model (SPM) to identify nine parameters. Santhanagopalan et al. (2007) used the Levenberg–Marquardt algorithm to identify five parameters in the Doyle–Fuller–Newman (DFN) and the single-particle-model under constant charge and discharge conditions. Scharrer et al. (2013) made use of a space-mapping parameter surrogate model to the DFN model to successfully identify three parameters. Their work made use of a Morris-One-At-A-Time sensitivity analysis to identify the three most sensitive parameters in the model.

Forman et al. (2012) performed parameter identification of 88 parameters using a genetic algorithm. To date, this is the latest attempt in estimating a significant number of parameters in the DFN model. Recently, Jin et al. (2018) also performed sensitivity analysis to identify the five most sensitive parameters. They then used Levenberg–Marquardt algorithm to estimate the values of these five parameters. A parallel genetic algorithm was used by Zhang et al. (2013) to identify 29 parameters in the pseudo-two-dimensional DFN model. They reported a computing time of 22.3 h to identify the 29 parameters. Uddin et al. (2016) estimated a total of three parameters in the DFN model using the differential evolution algorithm.

Previous works have reported solution times ranging from 22 h up to three weeks. This work drastically accelerates the parameter estimation of several parameters in the single particle and Newman models. In the work of Forman et al. (2012) it was stated that the parameter estimation took approximately three weeks.

In this work, the estimation of parameters in the DFN model took approximately 14 h. This great speed-up is due to the sophisticated minimization algorithm used to perform the parameter estimation. This work makes use of several minimization algorithms and continuously switches between them to accelerate convergence and avoid local minima.

2. Electrochemical model

The dynamics of lithium-ion batteries is of a highly multi-physics nature. The physics of the processes in Li⁺ cells are governed by strongly coupled, highly non-linear system of partial differential equations.

This section presents the model used in parameter identification.

Although simplifications can be made to the mathematical model of electrochemistry in the Li⁺ battery, the simulation of such simplified processes is still computationally expensive. For this reason, an efficient implementation of the mathematical model is needed. Owing to the nature of materials inside a cell several simplifying assumptions have been made, often applied in the field of battery modelling, to enable computational simulation of the electrical and chemical processes inside a cell. One crucial assumption is made: all electrode particles are spheres of radius R_s,i, where i ∈ {a, c} denotes the anode and cathode domain. This results in a simplified one-dimensional diffusion equation, which implies uniform superficial current and constant isotropic diffusion inside. Thus, the entire equation system may be quickly solved in contrast to full 3D-simulations, while accurately describing the insertion process.

Figure 1 shows the model domain schematically, including the layered structure of a cell, as well as the sub-domains annotations.

The DFN model: Each electrode is represented by homogenously distributed spherical particles as the limiting factor, connected via electrolyte. The model equation system results in:

(1) ∂cs,i∂t−1r2∇⋅(r2Ds,i∇cs,i)=0inΩs,i×Ωiεℓ∂cℓ∂t−∇⋅(RTF2t+κ(cℓ)cℓ∇cℓ+t+κ(cℓ)F∇ϕℓ)=AijBV in Ωi−∇⋅(RTF(2t+−1)κ(cℓ)cℓ∇cℓ+κ(cℓ)∇ϕℓ)=FAijBV,in Ωi

where c_s denotes the concentration of lithium inside the solid electrode, r is the spheres’ radial dimension and D_s is the solid diffusion coefficient in electrode i, with Ω_s,i = (0, R_s,i). Owing to the homogenous distribution of the particles and the assumption of a small dimension orthogonal to the layered structure, a single one-dimensional cut through the electrolyte domain models the electrolyte geometry, i.e. Ω_ℓ = (0, L), Ω_a = (0, L_a) and Ω_c = (L – L_c, L) (Figure 1).

The constant inner surface A_i = 3ε_s/R_s,i arises as the constant particle surface to particle volume ratio, ε_ℓ and ε_s denote the active volume fraction in the liquid and solid phase, respectively, t₊ is the charge transfer constant and κ(c_ℓ) is the conductivity in the liquid phase.

Taking into account all electrolyte and electrode quantities allows setting the particle boundary condition to:

(2) Ds,i∂cs,i∂n→=jBV,i(ϕℓ,ϕs,UOCP(cs))  on Γs,ijBV,i=i0(cs)(cℓcℓ,0exp(αFRT(ϕs−ϕℓ−UOCP(cs)))−exp(−(1−α)FRT(ϕs−ϕℓ−UOCP(cs)))),

where ϕ_s denotes the electrode potential and U_OCP(c_s) is the open circuit potential of the electrode at a given lithium concentration c_s. The model of U_OCP(c_s) used in this work is based on the Redlich–Kister expansion as introduced by Karthikeyan et al. (2008):

(3) UOCP(ξ)=E^0+RTFln⁡(1−ξξ)+RTF∑k=0nAk((2ξ−1)k+1−2ξk(1−ξ)(2ξ−1)1−k),

where ξ = c_s/c_total is a measure for the lithiation state of an intercalation electrode. Using the definition of U_OCP(ξ), Pichler (2018) derived the exchange current density on the basis of activity functions derived from transition theory:

(4) i0(ξs):=kBVexp⁡(FRT((ξs−α)UOCP(ξs)−∫0ξsUOCP(x)dx))

We assume constant behavior of the electronic quantities in the solid domain. This permits us to state the conservation of electric charge and current as:

(5) ∫ΩiFAejBV,idx=Iapp

where I_app is the applied current to the cell and A_e denotes the electrode cross section area.

In addition, we capture the effect of electrolyte losses by an ohmic resistance R_I, such that we may state the cell voltage u_cell as the algebraic condition:

(6) ucell=IappRI+ηc+ηa,

where the cell electrodes’ overpotentials η_i = ϕ_s,i – ϕ_ℓ – U_OCP(c_s,i) are used to simplify the notation.

Although most of this equation system is standard in literature, the final form differs in the liquid activity represented by cℓcℓ,0, that is only multiplied with one of the two exponential branches, whereas in the literature it is most often multiplied with both. The presented version arises from the distinction of the equilibrium activity a_ℓ,0 and the non-equilibrium activity a_ℓ that is probably neglected or missed in other works.

3. Experimental setup

The measured voltage data obtained through cycling a Panasonic NCR18650B commercial cell is used in this work.

For long term behavior, the cell is first charged at C/3 rate (C-rate = 3.35 A) until the voltage reaches 4.113 V followed by a constant voltage charge at 4.113 V until the current tapered down to 160 mA (≈C/20 rate), then discharged at C/3 rate until 3.498 V again followed by a constant voltage discharge at 3.498 V for 40 min or until the current dropped to 160 mA, respectively. This is repeated three times, afterwards two full capacity estimation cycles according to the data sheet are executed: the cell is charged at C/2 rate until 4.2 V, constant voltage charged at 4.2 V until the current dropped to C/50 rate and discharged at 1 C rate until 2.5 V. The cell is then charged to 3.498 V again and discharged to a specific SoC level for a total of seven cycles (85 per cent, 75 per cent, 65 per cent, 55 per cent, 45 per cent, 35 per cent and 25 per cent). At each level a set of current pulses are applied such that the short term dynamic behavior of the cell is reflected as much as possible in the voltage. The pulse sequence subsequently applies C/5, 1.25 C and 1.35 C pulses in charge (+) and discharge (−) direction for 10 s, followed by 15 min rest after each pulse. The pulse sequence ends with a combined 5 s-pulse sequence of +C/5, +C/5, −C/5, −C/5, −1.35 C, +1 C with 5 s rest in-between and a discrete discharge/charge stair profile of 0.2 C, 0.35 C, 0.5 C, 0.75 C, 1.25 C for 10 s per level.

Figure 2 shows the voltage and current measured throughout this time of roughly 2.5 days. All tests were carried out using an Arbin BT-2000 battery testing system and Memmert incubator with Peltier cooling (model IPP600) for maintaining the temperature at 25°C by forced air cooling.

4. Parameter identification

4.1 Framework

The traditional approach to solve the parameter identification problem involves minimizing the difference between the measured response and predicted response. If the cell voltage curve obtained through experimental measurement is V_E(t), where the time t varies between the start of the curve at t = 0 and its end at t = T, and the cell voltage curve obtained by solving the mathematical model for a given parameter set π is V(π,t), then the correct parameter set π can be estimated by solving the optimization problem given by:

(7) πest=arg⁡min⁡π∈Π∫0T(VE(t)−V(π,t))2dt.

The minimization algorithm subsequently updates the parameter set π to minimize the error norm. It should be mentioned that each computation of that error norm requires the solution of the mathematical model using the given parameter set π. In the case of an infeasible parameter set π_inf, that will lead to an inadmissible system state (e.g. c_ℓ(t_crash) ≤ 0) at some time t_crash, the simulation result V(π_inf, t) will be set equal to zero for t ≥ t_crash for practical application of the integration in the range of (0, T ).

An algorithm that can efficiently minimize the error with a few model evaluations is very appealing. This minimization algorithm must be robust and should be able to avoid local minima. For this reason, a newly developed hybrid optimizer is used to solve the above optimization problem.

4.2 Hybrid optimizer

Owing to the large computational time required to solve the mathematical model and because of the non-linearity of the cost-function space, an efficient and robust minimization technique is needed. The minimization technique in this work is a single objective hybrid optimizer (SOHO). The SOHO algorithm features three individual algorithms. The three algorithms are the single objective variants of the NSGA-III (Deb and Jain, 2014), NSDE-R (Reddy and Dulikravich, 2019) and MOEA-DD (Li et al., 2015). It is well known from the no free lunch theorem that no single algorithm is superior over another for an entire problem set. This means that the superiority of one algorithm over another for a problem set is paid for by the loss of its superiority over another problem set. This drives the need to couple several optimization algorithms to increase their robustness over a larger set of problems.

The SOHO is initialized with one of the three previously stated algorithms. Each algorithm operates till convergence. If stagnation is detected, an alternative algorithm is selected randomly from the remaining two. This random selection of algorithm adds a stochastic nature to search process and avoids user bias. All runs in this work were initialized with the NSDE-R algorithm. Figure 3 shows the three algorithms and the switching scheme. This hybridization allows SOHO to avoid local minima and increases the convergence rate to the global minimum.

The NSGA-III uses simulated binary crossover (Deb and Agrawal, 1995) and polynomial mutation (Deb, 2001) to perform the recombination. The parents to be mated are selected randomly from the entire population set. The NSDE-R uses the “rand/1/bin” (Robič and Filipič, 2005) mutation to perform the recombination where the parents to be mated are randomly selected from population set of unique members. The MOEA-DD also uses the same recombination operators as the NSGA-III algorithm but selects its parents randomly from the N best members. More details on these algorithms may be found in Deb (2001).

4.3 Parameter set

The DFN-model used in this work is defined using 44 parameters. Table AI in the Appendix shows the parameters to be identified for the model. The parameters to be estimated are: the separator resistance, along with the particle radii, diffusion coefficients, reaction rates and active mass of both the cathode and the anode, electrode area, separator porosity and the tortuosity of the cathode, anode and separator. A total of 15 terms (n = 15) in the Redlich–Kister expansion are used to define the OCP curve for each the anode and the cathode. The first RK coefficient for the anode is always set to zero because of linear dependency to α_{0, C}. Thus, the total number of RK coefficients is 29 for both the anode and the cathode.

4.4 Problem setup

As previously mentioned, the parameter estimation problem is solved by minimizing the L²-norm of the difference between the calculated and measured voltage curves. The calculated curve is obtained by solving the mathematical model while the measured voltage curve is obtained experimentally. It should be mentioned that the so-called “inverse crime” (Wirgin, 2004) is avoided in this work because the two voltage curves are obtained using different methods and because of the inherent measurement errors present in the experimentally obtained voltage curve.

The minimization is performed using the SOHO algorithm. The SOHO algorithm will search for the parameters, within a user-specific bound, that best minimizes this L²-norm. The lower and upper bounds for each of the parameters to be estimated in the model is given in Table AI. It should be mentioned that the bounds on each variable are conservative and larger than usual. This is to mimic the lack of prior knowledge about the parameters. The initial values of each parameter were randomly selected using the SOBOL’s algorithm (Sobol, 1967).

Owing to the large allowable range for most parameters, the optimization algorithm should first efficiently search a large parameter space but then must focus its search on a smaller region where there is a greater chance of finding the global minimum.

In this respect, Forman et al. (2012) divided the optimization problem into a global optimization run followed by local optimization run. Solving two separate optimization problems greatly increases the computational cost and time. This work makes efficient use of the recombination operators to solve both optimization problems in a single run. The crossover (Deb and Agrawal, 1995) and mutation distribution indices (Deb, 2001), η_c and η_m, control the proximity of the new candidate solution to its parents. A higher value of each index leads to a solution that is closer to its parents. The trade-off between global and local search is controlled by adapting the distribution indices as a function of generations. Each distribution index linearly increased from a value of 1 to 50 as function of generations. This leads to a more global search at the beginning which then gradually becomes a local search. The SOHO algorithm was run for a total of 1,000 generations, although in all cases, the minimum was found in less than 500 generations.

The SOHO algorithm is parallelized in a master-slave arrangement. The master node performs all optimization computation (recombination, selection, etc.) while each slave node solves the mathematical model. A total of 100 parallel runs (i.e. 100 slave nodes) are used throughout this work.

It should be mentioned that the solution of the DFN model was terminated if either the time step became less than 10⁻⁶ or if the maximum allowable working time was exceeded. The maximum allowable time was set as twice the average computing time. This greatly reduces computing time as infeasible parameter combinations runs are not evaluated. These termination criteria add additional degree of non-linearity and discontinuity to the cost function space. It also adds several “flat” regions where the gradient is zero. For this reason, a gradient based method will find it very difficult to converge to the correct values of the model parameters. The SOHO algorithm is not affected by any of these function space modifications.

5. Results

Previous results show that a single particle model is able to accurately model the battery response. For certain cases, e.g. very high currents, the single particle model may not be able to accurately represent the Li⁺ cell dynamics and a more complete model, such as the DFN-model, may be required.

The DFN model was defined using a total of 44 parameters. Figure 4 shows the estimated voltage and measured voltage using the 44 converged parameters.

It can be seen that the results of the DFN model are similar to those measured. The charging and discharge peaks coincide well for the entire time range.

Figure 5a shows the convergence history for the DFN model estimation problem.

Here, the residual is seen to sharply decrease within the first 10,000 evaluations. Figures 5b, 5c and 5d show error distributions at three different locations along the convergence history. It can be seen that even parameter sets in the early regions of the convergence history (Case 1) have majority of the errors within 1 per cent, with a significant number of them centered close to zero. In the Case 3, a large fraction of the errors is within 0.25 per cent.

Detailed Results may be found in the Appendix. The error statistic and the convergence information of the three selected cases of the DFN model are shown in Table AII. Even though the computing time of the DFN-model is of considerable magnitude, the SOHO algorithm is able to estimate the parameters shown in Table AIII in the DFN model in less than one day.This is a significant improvement in convergence time over the previous studies, which took approximately three weeks to obtain converged results. It should be mentioned that the computing time for the Newman model used in Forman et al. (2012), i.e. 63 s, is similar to the DFN model used in this work (30 s on average).

6. Conclusion

This work efficiently solves the parameter identification problem to match the voltage obtained using the Doyle–Fuller–Newman model and from experimental measurements. A total of 44 parameters are used to define the DFN model. The minimization of the error between computed and measured voltage was performed using an efficient single objective hybrid optimizer. The parameters of the DFN model were identified within one day. After identification, the model had a mean absolute error and root mean squared difference below 10 mV.

Figures

Figure 1.

Schematic view of the battery cell sub-domains

Figure 2.

The measured voltage (top) and the applied current (bottom) used in experimental study

Figure 3.

Algorithms and switching between the algorithms in the SOHO suite

Figure 4.

Measured and estimated voltage response obtained using the converged Newman model

Figure 5.

Parameter estimation of the DFN model showing (a) the convergence history of the SOHO algorithm; (b) error probability distribution of Case 1; (c) error probability distribution of Case 2; and (d) error probability distribution of Case 3

Table AI.

Parameters and admissible ranges in the DFN model

Variable descriptor	Symbol	Min	Max
Electrode area	AE	0	2
Cathode tortuosity	τC	0	1
Anode tortuosity	τA	0	1
Separator tortuosity	τS	0	1
Separator porosity	εs	0	1
Separator resistance	RI	0	1
Anode initial SoC	ξa,0	0	1
Cathode particle radius	rC	1.0E-8	1.0E-5
Anode particle radius	rA	1.0E-8	1.0E-5
Cathode diffusion coefficient	DC	0	1
Anode diffusion coefficient	DA	0	1
Cathode reaction rate	kC	−20	100
Anode reaction rate	kA	−20	100
Cathode active mass	mC	0	0.051
Anode active mass	mA	0	0.033
Cathode kth RK coefficients	Ak, C	−8	8
Anode kth RK coefficients	Ak, A	−8	8

Table AII.

Error statistics of the three selected cases

	Case 1	Case 2	Case 3
Evaluations to convergence	8,700	21,100	83,500
Approximate time to convergence (s)	5,455	13,330	52,354
Mean absolute error (mV))	18.91	15.82	6.47
Relative to measurement (%	0.492	0.326	0.276
Root mean squared difference (mV)	24.49	10.95	7.16
Normalized to measurement mean (%)	0.58	0.399	0.331

Table AIII.

Converged values of the parameters in the Newman model

Variable descriptor	Symbol	Value
Electrode area	AE	1.87
Cathode tortuosity	τC	8.29E-3
Anode tortuosity	τA	5.86E-1
Separator tortuosity	τS	3.20E-1
Separator porosity	εs	9.69E-1
Separator resistance	RI	4.68E-2
Anode initial SoC	ξa,0	7.17E-1
Cathode particle radius	rC	2.36E-6
Anode particle radius	rA	4.30E-6
Cathode diffusion coefficient	DC	9.98E-4
Anode diffusion coefficient	DA	5.00E-5
Cathode reaction rate	kC	−2.10
Anode reaction rate	kA	−9.25
Cathode active mass	mC	5.09E-2
Anode active mass	mA	3.26E-2
Cathode Redlich–Kister coefficients	A1, C	3.77
	A2, C	−2.9E-1
	A3, C	3.50
	A4, C	4.35E-1
	A5, C	4.24
	A6, C	3.64
	A7, C	3.41
	A8, C	−2.46
	A9, C	4.59
	A10, C	1.33
	A11, C	1.43
	A12, C	4.59
	A13, C	−1.9E-2
	A14, C	2.54
	A15, C	−4.94
Anode Redlich–Kister coefficients	A1, A	−4.85
	A2, A	2.76
	A3, A	−4.24
	A4, A	7.35E-1
	A5, A	2.29
	A6, A	−4.30
	A7, A	−3.89
	A8, A	−6.45E-1
	A9, A	3.66
	A10, A	2.50E-1
	A11, A	4.22
	A12, A	−5.86E-1
	A13, A	−8.20E-1
	A14, A	4.98