Identifiability and parameter estimation of the single particle lithium-ion battery model

TL;DR

This paper analyzes the identifiability and parameter estimation of the single particle model for lithium-ion batteries, revealing only six independent parameters and demonstrating practical estimation from impedance data with good predictive accuracy.

Contribution

It provides a comprehensive analysis of the model's identifiability and introduces a practical method for parameter estimation from experimental data.

Findings

Only six independent parameters in the model.

Model is structurally identifiable under certain conditions.

Estimated model predicts voltage with a maximum error of 20 mV.

Abstract

This paper investigates the identifiability and estimation of the parameters of the single particle model (SPM) for lithium-ion battery simulation. Identifiability is addressed both in principle and in practice. The approach begins by grouping parameters and partially non-dimensionalising the SPM to determine the maximum expected degrees of freedom in the problem. We discover that, excluding open circuit voltage, there are only six independent parameters. We then examine the structural identifiability by considering whether the transfer function of the linearised SPM is unique. It is found that the model is unique provided that the electrode open circuit voltage functions have a known non-zero gradient, the parameters are ordered, and the electrode kinetics are lumped into a single charge transfer resistance parameter. We then demonstrate the practical estimation of model parameters from measured frequency-domain experimental electrochemical impedance spectroscopy (EIS) data, and show additionally that the parametrised model provides good predictive capabilities in the time domain, exhibiting a maximum voltage error of 20 mV between model and experiment over a 10 minute dynamic discharge.

Tables

Table 1. Table I: Parameters used for generating synthetic electrochemical impedance data for an LCO cell reported in  [ 1 ] .

Parameter	Units	Anode	Cathode
Electrode thickness ${\delta }_i$	$\mu m$	$73.5$	$70.0$
Particle radius $R_i$	$\mu m$	$12.5$	$8.5$
Active material volume fraction ${ϵ}_i$	-	$0.4382$	$0.3000$
Li diffusivity in active material $D_i$	$m^2{s}^{-1}$	$5.5\text{×}{10}^{-14}$	$1.0\text{×}{10}^{-11}$
Reaction rate constant $k_i$	$m^{2.5}{\mathrm{mol}}^{-0.5}{s}^{-1}$	$1.764\text{×}{10}^{-11}$	$6.667\text{×}{10}^{-11}$
Max. active material concentration $c_i^{max}$	$\mathrm{mol}{m}^{-3}$	$\mathrm{30\hspace{0.17em}555}$	$\mathrm{51\hspace{0.17em}555}$
Electrode surface concentration $𝒜$	${\mathrm{cm}}^{-2}$	$982$	$982$
Electrolyte concentration $c_e$	$\mathrm{mol}{m}^{-3}$	$1000$

Table 2. Table II: RMS and maximum error between the measured and predicted impedance at 10   % times 10 percent 10\text{\,}\mathrm{\char 37\relax} , 50   % times 50 percent 50\text{\,}\mathrm{\char 37\relax} and 80   % times 80 percent 80\text{\,}\mathrm{\char 37\relax} DoD.

	Single DoD estimation		Combined DoDs estimation
DoD	RMS Error	Max Error	RMS Error	Max Error
$10\%$	$3.28m\Omega$	$5.60m\Omega$	$6.66m\Omega$	$16.73m\Omega$
$50\%$	$9.56m\Omega$	$21.36m\Omega$	$47.88m\Omega$	$95.16m\Omega$
$80\%$	$6.48m\Omega$	$11.32m\Omega$	$6.84m\Omega$	$12.82m\Omega$

Equations

$$\frac{\partial c_i}{\partial t}=\frac{D_i}{r_i^2}\frac{\partial }{\partial r_i}\left(r_i^2\frac{\partial c_i}{\partial r_i}\right),$$

$${\frac{\partial c_i}{\partial r_i}|}_{r_i=0}=0\text{and}{D}_i{\frac{\partial c_i}{\partial r_i}|}_{r=R_i}=-j_i,$$

$$c_i(0,r_i)=c_i^0,$$

$$j_{-}=\frac{+I}{a_{-}{\delta }_{-}\mathcal{F}\mathcal{A}}\text{and}{j}_{+}=\frac{-I}{a_{+}{\delta }_{+}\mathcal{F}\mathcal{A}},$$

$$V=U_{+}\left(x_{+}^s\right)-U_{-}\left(x_{-}^s\right)+{\eta }_{+}-{\eta }_{-}.$$

$$x_i^s=c_i^s/c_i^{max}.$$

$$j_i=\frac{i_{0,i}}{\mathcal{F}}\left[\exp \left(\frac{{\alpha }_i\mathcal{F}}{RT}{\eta }_i\right)-\exp \left(\frac{-\left(1-{\alpha }_i\right)\mathcal{F}}{RT}{\eta }_i\right)\right].$$

$$i_{0,i}=k_i\mathcal{F}\sqrt{c_e}\sqrt{c_i^s}\sqrt{c_i^{max}-c_i^s}.$$

$${\eta }_i=\frac{2RT}{\mathcal{F}}{\sinh }^{-1}\left(\frac{j_i\mathcal{F}}{2i_{0,i}}\right).$$

$$\frac{\partial {\bar{x}}_i}{\partial t}=\frac{D_i}{R_i^2}\frac{1}{{\bar{r}}^2}\frac{\partial }{\partial \bar{r}}\left({\bar{r}}^2\frac{\partial {\bar{x}}_i}{\partial \bar{r}}\right),$$

$${\frac{\partial {\bar{x}}_i}{\partial \bar{r}}|}_{\bar{r}=0}=0\text{and}{\frac{\partial {\bar{x}}_i}{\partial \bar{r}}|}_{\bar{r}=1}=\frac{-R_i}{D_i{c}_i^{max}}{j}_i,$$

$${\bar{x}}_i(0,\bar{r})=0.$$

$${\frac{\partial {\bar{x}}_{+}}{\partial \bar{r}}|}_{\bar{r}=1}$$

$${\frac{\partial {\bar{x}}_{-}}{\partial \bar{r}}|}_{\bar{r}=1}$$

$${\eta }_{+}$$

$${\eta }_{-}$$

$${\tau }_{+}^d$$

$${\tau }_{+}^k$$

$$Q_{+}^{th}$$

$$\frac{\partial {\bar{x}}_i}{\partial t}=\frac{1}{{\tau }_i^d}\frac{1}{{\bar{r}}^2}\frac{\partial }{\partial \bar{r}}\left({\bar{r}}^2\frac{\partial {\bar{x}}_i}{\partial \bar{r}}\right)$$

$${\frac{\partial {\bar{x}}_i}{\partial \bar{r}}|}_{\bar{r}=0}=0\text{and}{\frac{\partial {\bar{x}}_i}{\partial \bar{r}}|}_{\bar{r}=1}=-\frac{{\tau }_i^d}{3Q_i^{th}}I,$$

$${\bar{x}}_i(0,\bar{r})=0.$$

$$V=U_{+}\left(x_{+}^s\right)-U_{-}\left(x_{-}^s\right)+{\eta }_{+}-{\eta }_{-},$$

$${\eta }_{+}$$

$${\eta }_{-}$$

$$\mathbf{\theta }={\left[{\tau }_{+}^d\frac{{\tau }_{+}^d}{3Q_{+}^{th}}\frac{{\tau }_{+}^k}{3Q_{+}^{th}}{\tau }_{-}^d\frac{{\tau }_{-}^d}{3Q_{-}^{th}}\frac{{\tau }_{-}^k}{3Q_{-}^{th}}\right]}^T.$$

$$\{\begin{array}{ll}{\tau }_{+}^d& ={\theta }_1\\ {\tau }_{+}^k& =({\theta }_1{\theta }_3)/{\theta }_2\\ Q_{+}^{th}& ={\theta }_1/(3{\theta }_2)\\ {\tau }_{-}^d& ={\theta }_4\\ {\tau }_{-}^k& =({\theta }_4{\theta }_6)/{\theta }_5\\ Q_{-}^{th}& ={\theta }_4/(3{\theta }_5)\end{array}$$

$$H(s,\mathbf{\theta })=H(s,{\mathbf{\theta }}^{\ast })\text{for almost all }s,$$

$$\frac{\partial {\bar{u}}_i}{\partial t}=\frac{1}{{\tau }_i^d}\frac{{\partial }^2{\bar{u}}_i}{\partial {\bar{r}}^2},$$

$${\bar{u}}_i\left(\bar{r}=0\right)=0\text{and}{\frac{\partial {\bar{u}}_i}{\partial \bar{r}}|}_{\bar{r}=1}-{\bar{u}}_i(\bar{r}=1)=\frac{-{\tau }_i^d}{3Q_i^{th}}I.$$

$${\bar{u}}_i(0,\bar{r})=0.$$

$$\frac{d^2{\bar{U}}_i(s,\bar{r})}{d{\bar{r}}^2}-s{\tau }_i^d{\bar{U}}_i(s,\bar{r})=0$$

$${\lambda }^2-s{\tau }_i^d=0\Rightarrow \lambda =\pm \sqrt{s{\tau }_i^d},$$

Full text

Nomenclature

Roman

$\mathcal{A}$

Electrode surface area, $\text{\,}\mathrm{m}^{2}$

$a_{i}$

Specific active surface area, $\text{\,}\mathrm{m}^{-1}$

$c_{e}$

Electrolyte lithium concentration, $\text{\,}\mathrm{m}\mathrm{o}\mathrm{l}\,\mathrm{m}^{-3}$

$c_{i}$

Lithium concentration (active material), $\text{\,}\mathrm{m}\mathrm{o}\mathrm{l}\,\mathrm{m}^{-3}$

$c_{i}^{0}$

Initial lithium concentration, $\text{\,}\mathrm{m}\mathrm{o}\mathrm{l}\,\mathrm{m}^{-3}$

$c_{i}^{max}$

Maximum lithium concentration, $\text{\,}\mathrm{m}\mathrm{o}\mathrm{l}\,\mathrm{m}^{-3}$

$c_{i}^{s}$

Surface lithium concentration, $\text{\,}\mathrm{m}\mathrm{o}\mathrm{l}\,\mathrm{m}^{-3}$

$D_{i}$

Electrode active material diffusivity, $\text{\,}\mathrm{m}^{2}\,\mathrm{s}^{-1}$

$\mathcal{F}$

Faraday’s constant, $\text{\,}\mathrm{C}\,\mathrm{m}\mathrm{o}\mathrm{l}^{-1}$

$H^{0}$

SPM transfer function, $\text{\,}\mathrm{\SIUnitSymbolOhm}$

$H_{i}^{d}$

Diffusion model transfer function, $\text{\,}\mathrm{A}^{-1}$

$I$

Applied current, $\text{\,}\mathrm{A}$

$i_{0,i}$

Exchange current density, $\text{\,}\mathrm{A}\,\mathrm{m}^{-2}$

$j_{i}$

Reaction rate, $\text{\,}\mathrm{m}\mathrm{o}\mathrm{l}\,\mathrm{m}^{-2}\,\mathrm{s}^{-1}$

$k_{i}$

Reaction rate constant, $\text{\,}\mathrm{m}^{2.5}\,\mathrm{m}\mathrm{o}\mathrm{l}^{-0.5}\,\mathrm{s}^{-1}$

$L_{j}$

Single DoD loss function, $\text{\,}\mathrm{\SIUnitSymbolOhm}^{2}$

$L$

Combined DoDs loss function, $\text{\,}\mathrm{\SIUnitSymbolOhm}^{2}$

$N_{\omega}$

Number of data points per EIS measurement

$N_{DoD}$

Number of EIS datasets at different DoDs

$Q$

Charge/discharge capacity, $\text{\,}\mathrm{C}$

$Q_{i}^{th}$

Theoretical electrode capacity, $\text{\,}\mathrm{C}$.

$R$

Gas constant, $\text{\,}\mathrm{J}\,\mathrm{m}\mathrm{o}\mathrm{l}^{-1}\,\mathrm{K}^{-1}$

$R_{ct}^{0}$

Charge-transfer resistance, $\text{\,}\mathrm{\SIUnitSymbolOhm}$

$R_{i}$

Active material particle radius, $\text{\,}\mathrm{m}$

$r_{i}$

Particle radial coordinates, $\text{\,}\mathrm{m}$

$\bar{r}$

Particle dimensionless radial coordinates

$T$

Temperature, $\text{\,}\mathrm{K}$

$t$

Time, $\text{\,}\mathrm{s}$

$U_{i}$

Open-circuit voltage, $\text{\,}\mathrm{V}$

$\bar{u}_{i}$

Change of variable $\bar{r}\bar{x}_{i}$

$V$

Cell terminal voltage, $\text{\,}\mathrm{V}$

$\bar{V}^{0}$

Voltage deviation from equilibrium, $\text{\,}\mathrm{V}$

$V_{i}$

Volume of electrode active material, $\text{\,}\mathrm{m}^{3}$

$\mathbf{x}_{0}$

Linearisation reference point

$x_{i}$

Lithium stoichiometry (active material)

$x_{i}^{0}$

Initial lithium stoichiometry

$x_{i}^{s}$

Surface lithium stoichiometry

$\bar{x}_{i}$

Change of variable $x_{i}-x_{i}^{0}$

$Z_{j}$

Measured cell impedance, $\text{\,}\mathrm{\SIUnitSymbolOhm}$

Greek

$\alpha_{i}$

Anodic charge-transfer coefficient

$\alpha_{i}^{0}$

OCV gradient w.r.t stoichiometry, $\text{\,}\mathrm{V}$

$\beta_{i}^{0}$

OCV gradient w.r.t capacity, $\text{\,}\mathrm{V}\mathrm{/}\mathrm{(}\mathrm{A}\,\mathrm{s}\mathrm{)}$

$\delta_{i}$

Electrode thickness, $\text{\,}\mathrm{m}$

$\epsilon_{i}$

Solid-phase volume fraction

$\eta_{i}$

Overpotential, $\text{\,}\mathrm{V}$

$\theta_{j}$

Grouped parameter

$\boldsymbol{\theta}$

Vector of grouped parameters, $\boldsymbol{\theta}\in\mathbb{R}^{6}$

$\boldsymbol{\tilde{\theta}}$

Vector of identifiable parameters, $\boldsymbol{\tilde{\theta}}\in\mathbb{R}^{3}$

$\boldsymbol{\theta}_{id}$

Vector of identified parameters

$\tau_{i}^{d}$

Diffusion time constant, $\text{\,}\mathrm{s}$

$\tau_{i}^{k}$

Kinetic time constant, $\text{\,}\mathrm{s}$

$\omega$

Angular frequency, $\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}\,\mathrm{s}^{-1}$

I Introduction

In addition to their widespread use in consumer electronics, the importance of lithium-ion batteries is increasing in large-scale energy storage applications, such as automotive, off-grid and grid applications. Such large-scale battery energy storage systems are required to operate for a decade or more. Battery degradation, which results in capacity and power fade over time, becomes a particularly important issue in applications with such long service life. However, battery management system (BMS) technologies, currently used in commercial applications for monitoring and controlling batteries to ensure their safety, have very limited insights into battery degradation.

Two important functions of a conventional BMS are (i) to monitor the battery voltage, current and temperature and ensure that these remain within their region of safe operation, and (ii) to estimate unmeasurable metrics such as State-of-Charge (SoC) and State-of-Health (SoH). The SoC denotes the remaining energy in the cell compared to a fully charged cell, while the SoH is usually defined as the capacity or resistance of an aged cell compared to that of a pristine cell. A mathematical model of a battery is required to infer SoC and SoH from the available measurements of voltage, current and temperature. Conventional BMSs commonly employ electrical equivalent circuit models (ECMs) consisting of a voltage source and a network of resistors and capacitors whose current-voltage response mimics that of a battery. The parameters of such ECMs may vary with SoC, SoH, current and temperature and must therefore be identified from experimental data under a large range of operating conditions. However, the battery behaviour predicted by ECMs cannot credibly be extrapolated outside the range of validity over which they have been parametrised. Also, accounting for battery degradation is challenging due to the lack of physical significance of the model parameters [1].

Motivations

In order to improve the lifetime performance of battery packs, researchers are investigating the application of physics-based electrochemical models instead of ECMs in the BMS, for example to enable electrochemical control of fast charging [2, 3, 4]. Due to their physical basis, such first-principle models can provide insights into battery ageing because they can more easily be coupled to degradation models. Physics-based models could be solved online in next-generation BMSs to enable health-aware management and control algorithms that use internal electrochemical safety limits (e.g. local overpotential or lithium concentration) instead of the commonly employed terminal voltage and current safety limits [5, 6].

However, this is a challenging research area from both the modelling and numerical solution perspective, and also the application perspective. One key challenge is whether such models are observable and can be used reliably for state estimation. Recent promising results indicate that electrochemical models could indeed be used in this way [7, 8, 9, 10, 11, 12, 1, 13]. Another important condition for the use of an electrochemical model in a BMS is that the parameters of the model must be credible, not just to recreate the measured voltage accurately, but so that internal states, such as overpotential, are also accurate. In general, parameters in these models are obtained from the literature (for example [14, 1]), and whilst this is useful for initial design studies of cell performance, for use in a BMS, the model must be parametrised specifically for the cells used in the pack. Relatively few authors have tackled this parameter estimation problem [15, 16, 17, 18, 19] although it has been investigated in the context of ECMs of battery dynamics [20, 21, 22, 23]. The present work is therefore motivated by the need to obtain meaningful model parameters for using electrochemical battery models for control and estimation purposes.

Contributions

In this paper we investigate the parameter estimation challenge in a simplified version of the ‘Doyle-Fuller-Newman’ pseudo-two-dimensional (P2D) lithium-ion battery model [24], the so-called single particle model (SPM) [25, 26], first to see whether groups of parameters are identifiable in principle, and then to investigate whether in practice they can be identified from simulated and experimental frequency-domain Electrochemical Impedance Spectroscopy (EIS) data [27].

The main contribution of our work is to show the importance of considering the structural and practical identifiability of model parameters before undertaking any parameter estimation. Attempting the identification of all model parameters simultaneously without considering their identifiability is likely to lead to ill-conditioned optimisation problems. An electrochemical model identified in such a way may indeed fit the battery voltage and current response well, but the estimated parameters will have no physical significance and therefore the internal states cannot be reliably used for control. We also demonstrate parameter estimation of an electrochemical battery model from frequency-domain data obtained from commercial cells, and this is a relatively new and unexplored research topic [28, 29].

Our first key result is that the parameter identifiability of any lithium-ion battery model, whether ECM or first-principle electrochemical model, is largely conditional on the slope of each electrodes’ open-circuit voltage (OCV) as a function of SoC. A flat OCV curve overshadows all dynamics and results in parameter unidentifiability. A second key message from this work is that physics-based models for lithium-ion batteries are over-parametrised and therefore one must try to identify the minimum number of identifiable grouped parameters prior to attempting any parameter estimation.

II The single-particle model

The single particle model (SPM) was first introduced by [25] and later extended to lithium-ion batteries by [26]. The SPM is an electrochemical battery model describing the cell thermodynamics, the diffusion of lithium in the active material of the electrodes, and the interfacial kinetics of the lithium intercalation/de-intercalation reaction at the electrode/electrolyte interface. In contrast with the more complex P2D model [24], electrolyte dynamics are neglected in the SPM and therefore the reaction rate is assumed uniform for each electrode. We chose the SPM in this paper because it is the ‘simplest’ electrochemical lithium-ion model and a good starting point for parameter identification studies. We anticipate that our methodology could be extended to other electrochemical models suitable for higher C-rates. Importantly, we do not anticipate that our main conclusions will change with respect to more complex models: the electrode OCV functions will remain key in paramaterisation, and the over-parametrisation issue will become more challenging. Although the validity range of the SPM may vary depending on the physico-chemical properties of the cell (e.g. electrolyte conductivity, porosity), it is commonly accepted in the literature that the SPM is valid at currents below 1C to 2C111In the battery literature the electrical current is commonly expressed as a C-rate, which is a value of current normalised by the cell capacity so that a current amplitude of 1C will discharge the cell in one hour., where concentration gradients in the electrolyte are negligible [5, 15, 30]. For many applications, including electric vehicles and grid storage, average C-rates are lower than 1C and this limitation is not particularly restrictive. In fact, we demonstrate later in this paper that reasonable voltage errors may be obtained with the linearised SPM on data including peaks in current of up to 6C. Finally, for the purpose of parameter estimation, in our work we used frequency-domain linear response EIS experimental data. Linear EIS data are recorded by applying a small AC current signal (with zero DC offset) to the cell with an amplitude much lower than 1C (typically around C/10) to ensure that the cell voltage response remains linear. It is therefore reasonable to use the SPM for modelling the battery impedance from linear EIS data with such low amplitude currents.

II-A Spherical particle diffusion model

In the SPM, the diffusion of lithium in the active material of each battery electrode $i$ (anode or cathode) is governed by the Fickian diffusion equation in spherical coordinates

[TABLE]

where $r_{i}$ denotes the radial coordinates, $c_{i}$ the lithium concentration profile and $D_{i}$ the lithium diffusion coefficient (assumed uniform and constant) in electrode $i$. The subscript $i$ either takes the value $+$ or $-$ to refer to the cathode or anode domain respectively. The diffusion equation is subject to Neumann boundary conditions at the particle centre $r_{i}=0$ and surface $r_{i}=R_{i}$,

[TABLE]

and the initial condition equation,

[TABLE]

where $c_{i}^{0}$ is the initial concentration profile in each particle. The molar flux $j_{i}$ in (2) expresses the rate of the lithium intercalation/de-intercalation reaction. In contrast to the Doyle–Fuller–Newman model [24] where algebraic equations must be solved to obtain the molar flux for each particle, the molar fluxes in the SPM can be directly related to the applied battery current $I$ through the relationships,

[TABLE]

where $a_{i}=3\epsilon_{i}/R_{i}$ is the specific active surface area in electrode $i$, $\epsilon_{i}$ is the volume fraction of active material in an electrode, $\mathcal{F}$ is Faraday’s constant, $\mathcal{A}$ is the electrode surface area which is assumed equal for both electrodes, and $\delta_{i}$ is the thickness of the electrode $i$. By convention, the battery discharge current is positive and the charging current is negative. The anode is the negative electrode upon discharge and the cathode is the positive electrode. Although strictly speaking these swap upon charging, we will in this paper use the widely accepted abuse of terminology that the ‘anode’ is the negative electrode and the ‘cathode’ is the positive electrode. The above molar flux expressions assume an electrochemical reaction involving a single electron transfer, which is the case for the lithium intercalation/de-intercalation reaction.

II-B Voltage measurement equation

The initial-boundary value problems (1-3) for each electrode constitute the dynamic part of the single-particle model. The battery terminal voltage $V$ is given by the nonlinear measurement equation,

[TABLE]

The anode and cathode OCV $U_{-}$ and $U_{+}$ respectively are empirical nonlinear functions of the surface stoichiometry $x_{i}^{s}$ of each particle defined according to,

[TABLE]

The overpotential $\eta_{i}$ is the voltage drop due to the departure from equilibrium potential associated with the intercalation/de-­intercalation reaction in each electrode. The relationship between the reaction rate $j_{i}$ and the overpotential $\eta_{i}$ is given by the Butler­-Volmer kinetics equation [24]

[TABLE]

The exchange current density $i_{0,i}$ depends on the reaction rate constant $k_{i}$ and the reactants and products concentrations, i.e. the electrolyte concentration $c_{e}$ (assumed constant in this model), and the lithium solid-phase concentration $c_{i}^{s}$ at the surface of the particle, through the relation

[TABLE]

By assuming that the anodic charge transfer coefficients in (7) are $\alpha_{i}=0.5$, i.e. anodic and cathodic charge transfer coefficients are equal, the overpotential $\eta_{i}$ can be expressed as a function of the reaction rate $j_{i}$ by

[TABLE]

II-C Identification of grouped parameters

First-principle models, such as the SPM, are usually over-parametrised in the sense that only a subset of parameters can be estimated from measured input-output data. A first necessary step towards credible parameter estimation is to reformulate the model in terms of the minimum number of parameter groups. In order to identify such groups of parameters, we first introduce the dimensionless radial coordinates $\bar{r}=r_{i}/R_{i}$, and the stoichiometry $x_{i}=c_{i}/c_{i}^{max}$, instead of the concentration in the governing equations. In addition, the change of variable $\bar{x}_{i}=x_{i}-x_{i}^{0}$ is introduced, with $x_{i}^{0}=c_{i}^{0}/c_{i}^{max}$ being the initial stoichiometry in each particle, in order to set the initial condition to zero in the governing equations. Note that the time independent variable $t$ is kept dimensional for simplicity because the diffusion time constants are different in each electrode. Introducing these dimensionless variables, the initial-boundary value problem (1-3) can be written as

[TABLE]

subject to the boundary conditions

[TABLE]

with initial conditions

[TABLE]

The boundary condition at the surface of the particle in each electrode can also be written in terms of the input current $I$ instead of the fluxes $j_{i}$ using relations (4) according to

[TABLE]

Similarly using relation (4) the expression for the overpotential in each electrode (9) becomes

[TABLE]

with the voltage measurement (5) expression unchanged. Six physically meaningful groups of parameters naturally arise in these equations. Indeed one can identify three groups of parameters for each electrode: a diffusion time constant $\tau^{d}_{i}$, kinetics time constant $\tau^{k}_{i}$ and the maximum theoretical electrode capacity $Q^{th}_{i}$. The expressions for these six physically meaningful groups of parameters are defined as follows

[TABLE]

Note that for convenience the cathode theoretical capacity $Q_{+}^{th}$ is defined negative by convention and without loss of generality, to yield the same model structure in both the anode and cathode. Substituting these grouped parameters into the governing equations results in the following diffusion equation for both the cathode and anode:

[TABLE]

subject to the boundary conditions

[TABLE]

with initial condition

[TABLE]

The voltage measurement equation remains

[TABLE]

with the cathode and anode overpotentials given by

[TABLE]

Inspecting (20) through (25), one can identify six groups of parameters $\theta_{j}$ that fully parametrise the SPM given by the parameter vector $\boldsymbol{\theta}\in\mathbb{R}^{6}$:

[TABLE]

Unlike the parameters defined in (17-19), these parameters $\theta_{j}$ appear only once and not as a product of each other in the model equations. Moreover, one can show that there is a one-to-one mapping between the six parameters $\theta_{j}$ and the six parameters (17)–(18)–(19) defined by:

[TABLE]

Therefore, assuming that the initial electrode stoichiometries $x_{i}^{0}$ (i.e. state-of-charge) are known and that the OCVs are known functions of the surface stoichiometry for each electrode, the six parameters contained in $\boldsymbol{\theta}$, or equivalently the six parameters defined in (17–19), are sufficient to fully parametrise the SPM.

III Structural identifiability

Six grouped parameters have been identified as sufficient to fully parametrise the SPM. However, this does not imply that these six parameters can be identified from the battery current-voltage response. Several approaches can be used to investigate the parameter identifiability of a dynamical model from input-output data. In this section, the so-called structural identifiability [31] of the SPM is discussed. Structural identifiability investigates the mathematical identifiability of the model irrespective of the identification data considered (which are also assumed noise-free). A definition of structural identifiability for linear time-invariant dynamic models that can be cast into a transfer function $H(s,\boldsymbol{\theta})$ parametrised by a vector of parameters $\boldsymbol{\theta}$ is given as follows [32, 33]:

Definition 1**.**

Consider a model structure $\mathcal{M}$ with the transfer function $H(s,\boldsymbol{\theta})$ parametrised by $\boldsymbol{\theta}\in\mathcal{D}\subset\mathbb{R}^{n}$ where $n$ denotes the number of parameters of the model. The identifiability equation for $\mathcal{M}$ is given by

[TABLE]

where $\boldsymbol{\theta},\boldsymbol{\theta}^{\ast}\in\mathcal{D}$. The model structure $\mathcal{M}$ is said to be

• •

globally identifiable* if (28) has a unique solution in $\mathcal{D}$,*

• •

locally identifiable* if (28) has a finite number of solutions in $\mathcal{D}$,*

• •

unidentifiable* if (28) has a infinite number of solutions in $\mathcal{D}$.*

The SPM is not a linear dynamical model because of its voltage measurement equation (23) incorporating nonlinear Butler-Volmer kinetics (24,25), and nonlinear OCV functions. Therefore Definition 1 cannot directly be used to investigate the SPM structural identifiability unless we first linearise the model by assuming a small perturbation around a fixed depth-of-discharge (DoD) point. Since we assume the nonlinear OCV functions are measured and known a priori (for example using the techniques described in [34, 35]), linearisation is a valid approach to obtain the parameters of the diffusion sub-models, which are linear, and to obtain a linearised approximation of the kinetics. The subject of identifiability and parameter estimation of the nonlinear kinetics term in the output equation would be an interesting topic for further research.

III-A Diffusion model transcendental transfer function

Definition 1 requires the model transfer function; we therefore first derive this for the spherical particle diffusion model and then combine it in subsequent sections with the transfer function of the linearised voltage measurement equation to derive the transfer function for the linearised SPM. The transfer function for the spherical diffusion model has previously been derived in the literature for determining the diffusion coefficient of various electrode materials in a laboratory experimental setup [36, 37] and studying the cell frequency response analytically [38], but is included here for completeness. In order to simplify the derivation, the change of variable $\bar{u}_{i}=\bar{r}\bar{x}_{i}$ is introduced into the initial-boundary value problem (20-22) describing the diffusion in a spherical particle. Under this change of variable (20) can be equivalently written

[TABLE]

and the boundary conditions (21) become

[TABLE]

Note that the homogeneous Neumann boundary condition at the centre of the particle can be reduced to a simpler homogeneous Dirichlet boundary condition so that $\lim_{r_{i}\rightarrow 0}c_{i}(r_{i})$ remains finite. By introducing this change of variable the initial condition (22) conveniently becomes

[TABLE]

Since the initial-boundary value diffusion problem (29)–(30)–(31) is linear, an equivalent transfer function can be determined without loss of generality. Taking the Laplace transform of (29) yields

[TABLE]

where $s$ is the frequency-domain Laplace variable. The characteristic equation for this differential equation is

[TABLE]

and its general solution is therefore

[TABLE]

with $A_{i}(s)$, $B_{i}(s)$ two constants (with respect to $\bar{r}$) to be determined using the boundary conditions (30). Substituting (34) into the particle centre boundary condition (30) at $\bar{r}=0$ yields

[TABLE]

And substituting (34) and (35) into the surface boundary conditions  (30) at $\bar{r}=1$ yields

[TABLE]

Substituting (36)–(35) into (34) yields the general solution

[TABLE]

The variable of interest is the surface stoichiometry $\bar{X}^{s}_{i}(s)$, rather than $\bar{U}_{i}(s,\bar{r})$, since this is the variable involved in the voltage measurement equation. By substituting $\bar{u}_{i}=\bar{r}\bar{x}_{i}$ into (37), evaluating at $\bar{r}=1$ and dividing by the input current $I(s)$, the transfer function $H_{i}^{d}(s)$ from current to surface stoichiometry for the spherical diffusion model is given by

[TABLE]

The cathode and anode diffusion transfer functions expressed in terms of the parameter vector $\boldsymbol{\theta}$ are therefore respectively

[TABLE]

and

[TABLE]

III-B Linearisation of the voltage measurement equation

The voltage measurement equation is a nonlinear function of the anode and cathode surface stoichiometry $x_{-}^{s}$ and $x_{+}^{s}$ respectively, and the input current $I$. Assuming that the input current amplitude $I$ remains small and the battery is operated close to its initial DoD, i.e. the stoichiometry in both electrodes remains close to the initial value $x_{i}^{0}$, the voltage equation can be linearised using a first-order Taylor series approximation about the reference point $\mathbf{x_{0}}=\left(x_{+}^{0},x_{-}^{0},I_{0}=0\right)$ according to

[TABLE]

where $\bar{V}^{0}=V-V(\mathbf{x_{0}})$ denotes the deviation of the voltage from the equilibrium voltage at the linearisation point,

[TABLE]

The partial derivative of the voltage $V$ with respect to the input current $I$ evaluated at the reference point is given by

[TABLE]

which simplifies to

[TABLE]

when substituting for $I_{0}=0$. The partial derivative with respect to the surface stoichiometry in the cathode and anode are given respectively by

[TABLE]

and

[TABLE]

Substituting for $I_{0}=0$ yields the simple expressions

[TABLE]

where $\alpha_{i}^{0}$ denotes the gradients of the OCV functions in each electrode, with respect to stoichiometry, at the linearisation point. Substituting the partial derivatives in (41) by their respective expression given by (44), (47) and (48) therefore yields the linearised voltage measurement equation

[TABLE]

where $R_{ct}^{0}(\boldsymbol{\theta})$ denotes the charge transfer resistance arising from the intercalation/de-intercalation reaction kinetics and defined according to

[TABLE]

III-C Transfer function of the linearised SPM

Taking the Laplace transform of the linearised voltage equation (49) and dividing by the input current $I$ yields the transfer function $H^{0}(s,\boldsymbol{\theta})=\bar{V}^{0}(s)/I(s)$ of the linearised SPM about the equilibrium point $\mathbf{x_{0}}=\left(x_{+}^{0},x_{-}^{0},I_{0}=0\right)$:

[TABLE]

Substituting for the anode and cathode diffusion transfer function $H_{+}^{d}$ and $H_{-}^{d}$ with their respective expressions (39) and (40) gives

[TABLE]

One can see that only the difference between the parameters $\theta_{3}$ and $\theta_{6}$ describing the cathode and anode kinetics appears in the charge transfer resistance term $R_{ct}(\boldsymbol{\theta})$ given by (50). As a result, there are an infinite number of pairs $(\theta_{3},\theta_{6})$ that will yield the same transfer function and only the lumped parameter $R_{ct}$ can be estimated using the linearised model at a given DoD. We have therefore reduced the parameter space to five parameters by combining the cathode and anode kinetics into the charge-transfer resistance term $R_{ct}$.

This assumes that the OCV slopes $\alpha_{+}^{0}$ and $\alpha_{-}^{0}$ are known parameters, but these are not directly measurable, since in practice, OCV can only be measured with respect to capacity (units Ah or Coulombs), not against stoichiometry (which is non-dimensional). By recognising that the theoretical capacity of an electrode is given by $Q_{i}^{th}=c_{i}^{max}\mathcal{F}V_{i}$ with $V_{i}=\epsilon A\delta_{i}$ (the volume of active material in the electrode), one can easily show that a variation of stoichiometry $\delta x_{i}$ is proportional to a variation of charge/discharge capacity $\delta Q$ according to

[TABLE]

Therefore the derivative of the OCV with respect to $x_{i}$ can be related to its derivative with respect to charge/discharge capacity $Q$, denoted $\beta_{i}^{0}$, according to

[TABLE]

Substituting $\alpha_{i}^{0}$ for $\beta_{i}^{0}$ using (54) in (III-C) and recalling the expressions of $Q_{i}^{th}$ in terms of $\theta_{j}$ in (27) yields the transfer function

[TABLE]

where the function $f$ defined as

[TABLE]

is only parametrised by the electrode diffusion time constant $\theta_{1}$ or $\theta_{4}$ in the cathode or anode respectively. The coefficients $\beta+^{0}$ and $\beta_{-}^{0}$ are the measurable and assumed known OCV gradients with respect to capacity in the cathode and anode respectively. This requires access to half-cell or reference electrode cell data, which may be obtained from commercially available cells as we have recently demonstrated [34, 35]. In the absence of individual electrode OCV data, it is not possible to parametrise the SPM directly, although it may be reasonable to use a 2-electrode OCV measurement combined with available OCV data from the literature for the graphite negative electrode to infer the positive electrode OCV function.

We have thus reduced the parameter estimation problem of the linearised SPM to three independent parameters: the cathode diffusion time constant $\theta_{1}=\tau_{+}^{d}$, the anode diffusion time constant $\theta_{4}=\tau_{-}^{d}$ and the charge-transfer resistance $R_{ct}$. We will subsequently define the vector $\boldsymbol{\tilde{\theta}}\in\mathbb{R}^{3}$ of identifiable parameters as

[TABLE]

The transfer function of the linearised SPM expressed in terms of the three parameters defined in $\boldsymbol{\tilde{\theta}}$ is given by

[TABLE]

Fig. 1 shows a Nyquist plot of the SPM frequency response predicted by (58) at $25\text{\,}\mathrm{\char 37\relax}$ DoD for the LCO cell parameters given in Table I. At low frequencies, the response features a $45\text{\,}\mathrm{\SIUnitSymbolDegree}$ slope which is representative of semi-infinite linear diffusion. This $45\text{\,}\mathrm{\SIUnitSymbolDegree}$ slope is commonly referred to as the ‘diffusion tail’ in the battery electrochemical impedance spectroscopy (EIS) community and is usually modelled using a so-called Warburg constant-phase element [27, 39]. The cell frequency response tends towards a capacitive behaviour with decreasing frequency as exhibited by the vertical asymptote on the Nyquist plot. This capacitive behaviour arises from the variation of the average particle concentration at low frequencies, as opposed to the surface concentration variation only, due to the finite-length nature of the diffusion problem. Although the low frequencies giving rise to this capacitive behaviour are usually avoided in the lithium-ion battery EIS literature, we argue that the departure from the semi-infinite diffusion $45\text{\,}\mathrm{\SIUnitSymbolDegree}$ tail actually provides a more informative dataset for the parameter estimation of the SPM. Fig. 1 also shows the effect of varying the diffusion time constants $\tau_{+}^{d}$ and $\tau_{-}^{d}$ in the cathode and anode respectively by $\pm$50\text{,}\mathrm{\char 37\relax}$$ from their nominal values. Generally, a larger diffusion time constant (slower diffusion) yields higher cell impedance as expected. Moreover a smaller diffusion time constant (faster diffusion) results in the capacitive behaviour occurring at higher frequencies; this is because faster diffusion results in easier variations of the average particle concentration. Finally, the frequency response is much more sensitive to the cathode diffusion than the anode diffusion at this DoD because of the flat anode OCV curve here, as shown by the low value of the anode OCV slope $\beta_{-}^{0}$ compared to $\beta_{+}^{0}$. This behaviour will be explained in more detail in the subsequent sections.

III-D Structural identifiability analysis

Definition 1 can now be applied to the SPM transfer function (58). To check the structural identifiability of the linearised model, we need to show that

[TABLE]

for almost all $s$, implies the parameters equality,

[TABLE]

Because the charge transfer resistance is the only additive parameter independent of $s$ on both sides of the equation, we clearly have $R_{ct}^{0}=R_{ct}^{0\ast}$, and the structural identifiability equation reduces to

[TABLE]

for almost all $s$. Since $f$ is a non-trivial function of $s$, this equality holds in the general case for most values of $s$, if and only if $\tau_{+}^{d}=\tau_{+}^{d\ast}$ and $\tau_{-}^{d}=\tau_{i}^{d\ast}$, and so the linearised SPM is structurally identifiable. There are however a few cases where other solutions exist:

• •

If $\beta_{+}^{0}=0$ (resp. $\beta_{-}^{0}=0$), then any pair $\left(\tau_{+}^{d},\tau_{+}^{d\ast}\right)$ (resp. $\left(\tau_{-}^{d},\tau_{-}^{d\ast}\right)$) satisfy the identifiability equation and the linearised SPM becomes unidentifiable. This makes sense because a ‘flat’ OCV function hides any diffusion dynamics effect of that electrode as shown in Fig. 1 for the anode diffusion when $\mid\beta_{-}^{0}\mid\ll\mid\beta_{+}^{0}\mid$ .

• •

If the magnitudes of the OCV function in each electrode are equal $\beta_{+}^{0}=-\beta_{-}^{0}=\beta^{0}$, then interchanging the diffusion time constants $\tau_{+}^{d}=\tau_{-}^{d\ast}$ and $\tau_{-}^{d}=\tau_{+}^{d\ast}$ satisfies the structural identifiability equation, and the linearised SPM is structurally identifiable provided the diffusion time constants are ordered.

In conclusion, the linearised SPM is structurally identifiable in the general case. In practice, the fact that a ‘flat’ OCV function results in unidentifiable parameters is not problematic since identification may be performed using data at several DoDs, ensuring the OCV functions have a significant slope in each electrode.

IV Frequency-domain parameter estimation

This section discusses the implementation of the estimation algorithm for the transfer function of the linearised SPM from frequency-domain EIS experimental data. The experimental methodology is discussed in section V. We define the vector $\boldsymbol{\theta}_{id}$ of parameters identified by the estimation algorithm in the general sense, where $\boldsymbol{\theta}_{id}$ could be equal to the full vector of parameters $\boldsymbol{\theta}$ given by (26) or a subset of parameters such as $\boldsymbol{\tilde{\theta}}$ given by (57). The estimation is performed by finding the best fit in the least-squares sense, i.e. by minimising the sum of squared errors between both the real and imaginary parts of the experimental and predicted impedances over a number of frequencies [32]. First, the loss function for the estimation at a single DoD is presented and extended to simultaneous estimation at several DoDs. Then, we show how all combined purely resistive parameters, including charge-transfer resistance $R_{ct}$, can be estimated separately by simple linear regression at a given DoD.

IV-A Single DoD parameter estimation

The measured impedance of a cell at a given frequency $\omega_{i}$ is denoted by the complex number $Z_{j}\left(\omega_{i}\right)=Z_{j}^{\prime}\left(\omega_{i}\right)+iZ_{j}^{\prime\prime}\left(\omega_{i}\right)$, where $Z^{\prime}$ and $Z^{\prime\prime}$ denote the real and imaginary part respectively, and the subscript $j$ denotes the DoD. Similarly, the impedance predicted by the SPM will be denoted $H_{j}\left(\omega_{i},\boldsymbol{\theta}_{id}\right)$ in the subsequent sections.

The optimal parameter estimate in the least-squares (maximum likelihood) sense is given by the estimator $\boldsymbol{\hat{\theta}}_{id}$ which minimises the loss function $L_{j}\left(\boldsymbol{\theta}_{id}\right)$ according to

[TABLE]

where the loss function is the sum over $N_{\omega}$ chosen frequencies of the squared magnitude of the error between predicted and measured complex impedance according to

[TABLE]

This is also equal to the sum of the squared real and imaginary part errors according to

[TABLE]

IV-B Combined DoDs parameter estimation

As mentioned in section III-D and shown in section V, the estimation of the parameter $\boldsymbol{\tilde{\theta}}$ based on the impedance data at a single DoD may be unidentifiable in practice due to the flatness of one of the electrode OCV functions. Therefore, the parameter estimation must be performed against impedance measured at several DoDs. As before, the optimal parameter estimate is given by the argument $\boldsymbol{\hat{\theta}}_{id}$ minimising the loss function according to

[TABLE]

where the loss function $L\left(\boldsymbol{\theta}_{id}\right)$ in this case is simply defined as the sum of the loss functions $L_{j}$ defined in (64) over $N_{DoD}$ levels of DoD, where $j$ denotes the DoD, according to:

[TABLE]

IV-C Estimation of purely resistive terms including charge-transfer resistance

Although the simultaneous parameter estimation at several DoDs could in theory improve the identifiability of the estimation problem (assuming the diffusion parameters are not a function of DoD), it also adds $N_{DoD}$ charge transfer resistances $R_{ct}^{0}$ to identify because of the dependency of $R_{ct}^{0}$ on stoichiometry (or equivalently DoD) in (50). The charge-transfer resistance $R_{ct}^{0}$ is an additive purely resistive term in the transfer function describing the cell impedance. Therefore, its only effect is to shift the impedance response along the real-axis on the Nyquist plot. In order to reduce the parameter space to be explored, a purely resistive term denoted $R_{0}$, combining the higher-frequency charge-transfer resistance $R_{ct}^{0}$, the ohmic contact resistance of the cell and any resistive effect due to passivation layers, can be estimated separately by simple linear regression and subtracted from the experimental frequency-response data. Recalling from section III-C that the semi-infinite diffusion-driven impedance at lower frequencies features a characteristic $45\text{\,}\mathrm{\SIUnitSymbolDegree}$ slope, the charge-transfer resistance can be estimated by fitting a $\mathrm{4}\mathrm{5}$° straight line to the experimental response at lower frequencies, and extrapolating to find the intercept of this with the real axis on the Nyquist diagram. Although the SPM could be extended to account for double-layer capacitance effects, this is beyond the scope of the present study. Therefore, the semi-circle at higher frequencies on the experimental EIS Nyquist plot, which arises from the parallel combination of charge-transfer resistance and double-layer capacitance, is ignored. Any high-frequency effects are considered as purely resistive and accounted for in the term $R_{0}$ estimated by the linear regression described below. We assume that any experimental high-frequency data points where charge-transfer processes are dominant (charge-transfer semicircle), have been removed; only data points at frequencies where diffusion processes are dominant remain. The linear regression for estimating the purely resistive term is illustrated in Fig. 5 for EIS data at $10\text{\,}\mathrm{\char 37\relax}$ DoD.

Denoting $x_{i}=Z^{\prime}\left(\omega_{i}\right)$ and $y_{i}=-Z^{\prime\prime}\left(\omega_{i}\right)$, the real part and negative imaginary part respectively of the experimental impedance at frequency $\omega_{i}$, the linear regression problem is

[TABLE]

The slope is set to $45\text{\,}\mathrm{\SIUnitSymbolDegree}$ to represent the diffusion behaviour as discussed previously. Therefore $\beta_{1}=1$ is known and the regression problem consists of estimating the real intercept $\beta_{0}$ ($=-R_{0}$) only. Defining the variable $\bar{y}_{i}=y_{i}-\beta_{1}x_{i}$, the regression problem can be written:

[TABLE]

This therefore constitutes an over-determined system of $N_{\omega}$ equations with the only unknown $\beta_{0}$, which can be written

[TABLE]

where $\bar{\mathbf{Y}}=\left[\bar{y}_{1}\!,\ldots\!,\bar{y}_{N_{\omega}}\right]^{T}\in\mathbb{R}^{N_{\omega}}$ is the vector containing the values of $\bar{y}_{i}$ at the $N_{\omega}$ frequencies, and $\mathbf{1}\in\mathbb{R}^{N_{\omega}}$ is a vector of all ones of length $N_{\omega}$. The value of $\beta_{0}$ can easily be determined in the least-squares sense by computing

[TABLE]

where the superscript $\mathbf{1}^{\dagger}$ denotes the pseudo-inverse of the vector $\mathbf{1}$ which is simply equal to a row vector with all components equal to $1/N_{\omega}$. This linear regression must be performed using only the experimental data points in the frequency range where low-frequency pseudo-capacitive effects arising from OCV variations are negligible. Therefore, very low frequency impedance data points must be discarded in the experimental dataset to perform a meaningful regression. In practice this was achieved by performing the regression first on all data points, and then on an iteratively reduced dataset where the lowest frequency impedance data point was removed until the goodness-of-fit was satisfactory. The coefficient of determination $R^{2}$ was used as an indicator to assess the goodness-of-fit and a threshold value of $0.98$ was chosen to stop the iterations.

V Results and discussion

V-A Parameter estimation using synthetic data

The parameter estimation algorithm was first tested against synthetic data generated using the linearised SPM with a known set of parameters for an LCO cell from the literature [1], which are summarised in Table I.

Fig. 2 shows contour plots of the natural logarithm of the loss function $L_{j}\left(\boldsymbol{\tilde{\theta}}\right)$, (64), at several DoDs for the synthetic data. It is assumed that the charge-transfer resistance $R_{ct}$ is known for each DoD since, as previously explained, it can be identified separately by straightforward linear regression. Therefore the two parameters to be identified are the cathode and anode diffusion time constants $\tau_{+}^{d}$ and $\tau_{-}^{d}$. The coordinates of the minimum values of the natural logarithm of (64) on these contour plots represent the optimal parameter estimates in the least-squares sense as defined in (62). The two cases for which the linearised SPM is structurally unidentifiable, as discussed in section III-D, can be witnessed by examining Fig. 2. Firstly, when the absolute value of the OCV functions are approximately equal in the cathode and anode, $\beta_{+}^{0}\approx-\beta_{-}^{0}$, such as $75\text{\,}\mathrm{\char 37\relax}$ DoD, the loss function clearly shows two minima: one minimum corresponds to the actual parameters, while the other minimum happens when the anode and cathode diffusion time constants are interchanged. Moreover, even if the absolute values of the electrode OCV slopes $\beta_{-}^{0}$ and $\beta_{+}^{0}$ are distinct (e.g. $5\text{\,}\mathrm{\char 37\relax}25\text{\,}\mathrm{\char 37\relax}95\text{\,}\mathrm{\char 37\relax}$ on Fig. 2), the loss function still features two local minima, although barely visible due to the narrow loss function along one axis. Therefore, an ambiguity persists and the anode and cathode dynamics cannot be distinguished. A possible explanation for this ambiguity is that the measured voltage is the difference between the anode and cathode OCV and not the absolute voltage of the electrode with respect to a reference electrode. Secondly, in the case where one OCV slope is an order of magnitude lower in one electrode compared to the other, the diffusion time constant of this electrode becomes unidentifiable, e.g. $5\text{\,}\mathrm{\char 37\relax}25\text{\,}\mathrm{\char 37\relax}95\text{\,}\mathrm{\char 37\relax}$ DoD in Fig. 2. For instance at $95\text{\,}\mathrm{\char 37\relax}$ DoD, the absolute value of the cathode OCV slope $\beta_{+}^{0}$ is small compared to the anode slope, and the minimum of the loss function is elongated along the $\tau_{+}^{d}$ axis, i.e. provided the value of $\tau_{-}^{d}$ is correct, almost any value for $\tau_{+}^{d}$ will result in a minimum value of the loss function.

These results suggest that the anode and cathode diffusion time constants cannot both be identified from frequency-data at a single DoD. In the best case, where the OCV slopes are equal for the chosen DoD, the values $\tau_{+}^{d}$ and $\tau_{-}^{d}$ can be determined but cannot be assigned to a specific electrode. Fig. 3 shows the natural logarithm of the loss function $L\left(\boldsymbol{\tilde{\theta}}\right)$ combining the EIS data at $5\text{\,}\mathrm{\char 37\relax}25\text{\,}\mathrm{\char 37\relax}75\text{\,}\mathrm{\char 37\relax}95\text{\,}\mathrm{\char 37\relax}$ DoD as defined in (66). This loss function shows a single global minimum corresponding to the actual set of parameters used to generate the synthetic data. The two time constants of the linearised SPM can therefore only be estimated unambiguously by considering EIS data at several DoDs. More specifically, complementary DoD values, where each electrode in turn presents a large OCV slope while the other is negligible, should be combined to ensure better identifiability. For instance, it is clear from Fig. 2 that combining EIS data at $5\text{\,}\mathrm{\char 37\relax}$ and $95\text{\,}\mathrm{\char 37\relax}$ DoD will result in a single minimum at the intersection of the lines of minimum loss function values for individual DoDs.

V-B Parameter estimation using experimental data

This section discusses the results of the parameter estimation using EIS experimental data measured from a commercial Kokam NMC cell. First, the experimental procedure used to measure the cell impedance and the OCV for both electrodes is discussed. The estimation of the charge-transfer resistance using linear regression is then briefly considered. Finally, the practical identifiability of the model diffusion time constants from the experimental impedance data is examined.

V-B1 Experimental setup

Electrochemical impedance spectroscopy (EIS) [27] is a widely used technique to measure the frequency response of electrochemical systems, including lithium-ion batteries. During a galvanostatic EIS (GEIS) experiment, the battery is excited with a small sinusoidal current and the associated voltage response is recorded. This is repeated at a number of frequencies. Potentiostatic EIS is a similar technique, where the voltage is used as the input sinusoidal signal and the recorded output is the current. The impedance of the cell at each frequency is calculated from the gain and phase shift between the input and output signals. EIS experimental data were obtained for a Kokam SLPB533459H4 $740$ mAh NMC cell at several DoDs using a BioLogic SP-150 potentiostat. The cell was kept at a constant $20\text{\,}\mathrm{\SIUnitSymbolCelsius}$ temperature in a Vötsch VT4002 thermal chamber in order to minimise the possible detrimental effects of temperature variations on the results of the parameter identification approach. The cell impedance was measured using single-sine GEIS, averaged over two periods, with a peak-to-peak AC current amplitude of $100\text{\,}\mathrm{mA}$, and $0\text{\,}\mathrm{mA}$ DC current. The impedance was measured at logarithmically-spaced frequencies with six frequencies per decade ranging from $5\text{\,}\mathrm{kHz}$ to $200\text{\,}\mathrm{\SIUnitSymbolMicro Hz}$. A resting time of 3 hours was applied before any EIS measurement to ensure that the cell was close to equilibrium. The EIS measurement using this experimental setup took approximately 9 hours for each DoD because of the low frequencies required. Clearly this is quite time consuming and ultimately one should aim at using the minimum number of low-frequency data and/or DoD linearisation points to accelerate the model parametrisation.

Fig. 4 shows the measured EIS data as a Nyquist plot for several DoDs ranging from $10\text{\,}\mathrm{\char 37\relax}$ to $90\text{\,}\mathrm{\char 37\relax}$. These impedance results are typical of lithium-ion batteries, which have smaller impedance at high frequencies (bottom left), arising from ohmic resistance and fast charge-transfer processes, and larger impedance at lower frequencies (top right), arising from slower diffusion processes and SoC variations. At high frequencies, the response features a depressed semicircle, which is characteristic of charge-transfer processes. As discussed in section IV-C, this part of the plot is discarded because charge-transfer processes are approximated by the resistor $R_{ct}^{0}$ only in (58). As expected, in the medium range of frequencies where semi-infinite diffusion is predominant, the frequency response shows the familiar $45\text{\,}\mathrm{\SIUnitSymbolDegree}$ diffusion slope. At very low frequencies the response tends towards a capacitive vertical line due to the variations of the average particle stoichiometry.

Fig. 7d shows the OCV of the Kokam SLPB533459H4 cell as a function of discharge capacity $Q$ and DoD measured experimentally. In order to measure the anode and cathode potential separately, a minimally invasive reference electrode consisting of a lithium-coated copper wire was inserted into the commercial pouch cell using the technique presented in [35]. The OCV of both electrodes were then measured using the galvanostatic intermittent titration technique (GITT) [34] by discharging the cell in $14.8\text{\,}\mathrm{mA}\text{\,}\mathrm{h}$ increments at C/10 and measuring the cell OCV after a one-hour voltage relaxation time to obtain 50 measurements points.

V-B2 Estimation of high-frequency purely resistive terms

The purely resistive term $R_{0}$ accounting for high-frequency effects, including charge-transfer, passivation layer and contact resistance, which shifts the diffusion-driven impedance to the right on the experimental EIS Nyquist plot, is estimated using the linear regression method discussed in section IV-C. Fig. 5 shows the results of the linear regression for the EIS data collected at $10\text{\,}\mathrm{\char 37\relax}$ DoD. As discussed previously, high frequency data points corresponding to charge-transfer resistance and double-layer capacitance are discarded and only data points at frequencies where diffusion effects are predominant are considered. The initial linear regression is performed on all low frequency data points but yields very poor fitting ($R^{2}=0.525$) due to the the pseudo-capacitive behaviour at very low frequencies arising from DoD variations. By iteratively removing the lower frequency data points, a good linear fit of the $45\text{\,}\mathrm{\SIUnitSymbolDegree}$ tail can be achieved ($R^{2}=0.995$). The purely resistive term $R_{0}$ is then approximated by finding the intercept of the linear regession curve with the real axis. The estimated values of $R_{0}$ as a function of DoD are reported on Fig. 6. Note that the cell ohmic resistance included in $R_{0}$ is approximately equal to $25\text{\,}\mathrm{m\SIUnitSymbolOhm}$ according to the intercept of the experimentally measured impedance with the real-axis of the Nyquist plot at high frequencies, Fig. 4. As shown, the resistive term $R_{0}$ including charge-transfer and contact resistance increases with DoD, which is consistent with previously reported Kokam NMC cell EIS data in the literature [40].

V-B3 Parameter estimation performance

Fig. 7 presents practical identifiability and parameter estimation results based on experimental impedance data at individual DoDs. Three characteristic values of DoD where chosen, $10\text{\,}\mathrm{\char 37\relax}50\text{\,}\mathrm{\char 37\relax}80\text{\,}\mathrm{\char 37\relax}$, in order to demonstrate the identifiability issues discussed in section III-D. At $10\text{\,}\mathrm{\char 37\relax}$ DoD, only the cathode OCV shows a significant slope, at $50\text{\,}\mathrm{\char 37\relax}$ DoD both the cathode and anode OCV are flat, and at $80\text{\,}\mathrm{\char 37\relax}$ DoD both electrodes feature a significant OCV slope. The experimental and fitted impedance responses are shown on Fig. 7e, f and g for the chosen values of DoD. The corresponding root-mean square (RMS) and maximum (Max) error between the predicted and measured impedance response for this single DoD parameter estimation algorithm are presented in Table II. A satisfactory fit is obtained in all cases with a root-mean square error of $3.28\text{\,}\mathrm{m\SIUnitSymbolOhm}9.56\text{\,}\mathrm{m\SIUnitSymbolOhm}6.48\text{\,}\mathrm{m\SIUnitSymbolOhm}$ for $10\text{\,}\mathrm{\char 37\relax}50\text{\,}\mathrm{\char 37\relax}80\text{\,}\mathrm{\char 37\relax}$ DoD respectively.

However, the uncertainties on the parameter estimates are large as exhibited by the contour plots of the loss function in Fig. 7a, b and c. At $10\text{\,}\mathrm{\char 37\relax}$ DoD the loss function is very sensitive to the cathode diffusion dynamics while anode dynamics have very little impact because of the very small slope of the anode OCV. Therefore the estimated anode diffusion time constant $\tau_{-}^{d}$ is highly uncertain and cannot be estimated from this data set alone. In contrast, at $80\text{\,}\mathrm{\char 37\relax}$ DoD the loss function is sensitive to both the anode and cathode dynamics because of the similar magnitude of the OCV slope in both electrodes. However, the minimum of the loss function is still relatively elongated along the anode parameter axis, which suggests that the uncertainty on the anode parameter remains larger than the cathode one. This confirms that impedance data at several DoDs must be combined to yield less uncertain parameter estimates.

At $50\text{\,}\mathrm{\char 37\relax}$ DoD, where the OCV is flat for both electrodes, the loss function is completely insensitive to the anode dynamics and somewhat more sensitive to the cathode dynamics. However, the contour plot of the loss function is actually shallow along the cathode parameter as well, implying that the uncertainty on both parameters is very large. Furthermore, the minimum of the loss function happens for a cathode diffusion time constant centred around $60\,000\text{\,}\mathrm{s}$ as shown in Fig. 7b. This is approximately an order of magnitude higher than the value of $\tau_{+}^{d}$ estimated using the other data sets at $10\text{\,}\mathrm{\char 37\relax}80\text{\,}\mathrm{\char 37\relax}$ DoD. Throughout this work, the implicit assumption was made that both diffusion time constants do not vary with DoD. However, Fig. 7b suggests that lithium diffusivity in the electrode active material could drastically change with DoD. Experimental measurements of lithium diffusivities in the active material of a Kokam NMC cell are reported in [41] and confirm this. More especially, both the anode and cathode diffusivities are an order of magnitude lower around $50\text{\,}\mathrm{\char 37\relax}$ DoD compared to high DoDs [41]. Another source of uncertainty that has been ignored in the model is the error related to the OCV measurement and the calculation of the OCV slopes $\beta_{+}^{0}$ and $\beta_{-}^{0}$. At $50\text{\,}\mathrm{\char 37\relax}$ DoD, the OCV slopes in both electrodes are very small and therefore the relative error on measurement is large, especially considering that numerical differentiation is highly sensitive to noise. Finally, because the OCV slopes are flat, both the anode and cathode diffusion dynamics have very limited impact on the cell terminal voltage at $50\text{\,}\mathrm{\char 37\relax}$ DoD. Therefore, the EIS data at $50\text{\,}\mathrm{\char 37\relax}$ DoD will be discarded from the parameter estimation using combined DoDs discussed subsequently.

Fig. 8 shows the contour plot of the cost function $L\left(\boldsymbol{\tilde{\theta}}\right)$ combining the experimental impedance data measured at $10\text{\,}\mathrm{\char 37\relax}80\text{\,}\mathrm{\char 37\relax}$ DoDs. Similarly to the results with synthetic data, Fig. 3, the loss function presents a single minimum, elongated along the anode parameter axis suggesting a larger uncertainty on the anode parameter compared to the cathode one. Although this suggests a ‘well-behaved’ estimation problem with a single parameter estimate, these results must be considered with caution. Indeed, the two combined datasets at $10\text{\,}\mathrm{\char 37\relax}80\text{\,}\mathrm{\char 37\relax}$ DoDs used to perform the parameter estimation were chosen somewhat arbitrarily. Although these parameter estimates yield a good impedance fit for these DoDs, it results in very approximate fitted impedance at other DoDs, e.g. $50\text{\,}\mathrm{\char 37\relax}$ DoD. As shown in Table II, the impedance prediction error is much higher at $50\text{\,}\mathrm{\char 37\relax}$ DoD with $47.88\text{\,}\mathrm{m\SIUnitSymbolOhm}$ RMS error compared to $6.66\text{\,}\mathrm{m\SIUnitSymbolOhm}$ and $6.84\text{\,}\mathrm{m\SIUnitSymbolOhm}$ at $10\text{\,}\mathrm{\char 37\relax}$ and $80\text{\,}\mathrm{\char 37\relax}$ DoDs respectively. Firstly, as previously mentioned the OCV slopes $\beta_{i}^{0}$ in both electrodes are assumed fixed and perfectly known for each DoD. The sensitivity of the predicted impedance response to the OCV slopes is not negligible. The computation of the OCV slope is highly sensitive to measurement noise because the measured voltage variations are relatively small between each DoD, on the order of a few millivolts. Additionally, the measured OCV might not reflect closely enough the actual cell OCV. The OCV measurements were performed on a different cell (of the same type and size) which was modified by inserting a reference electrode. Although this technique is minimally invasive, it nonetheless may slightly modify the cell behaviour. Finally, electrode OCV is affected by temperature variation and hysteresis effects [34], which introduce further uncertainties. Secondly, our parameter estimation algorithm combining experimental data at several DoDs assumes that the diffusion time constants in both electrodes remain constant with DoD. Although solid-phase diffusivities are usually assumed constant with respect to active material lithiation in the lithium-ion battery electrochemical modelling community, e.g. [1, 14, 15, 42, 43], this assumption may introduce additional inaccuracies on the parameter estimate, especially for NMC cells [41]. All of these challenges provide fruitful avenues for further exploration.

V-C Model validation in the time domain

To demonstrate the performance of the linearised SPM using our frequency domain parameter estimation approach, we also include here a brief comparison between the modelled and measured voltage response in the time domain. The measurements were made on the same Kokam cell type as previously described, using the same equipment that was used for EIS measurements. The current demand was specified based on a dynamic loading profile at $1\text{\,}\mathrm{Hz}$ derived from an electric vehicle drive cycle, see [44]. The maximum current was $4.4\text{\,}\mathrm{A}$, which corresponds to a fairly high C-rate (almost 6C). The same current input signal used in the experimental measurements was provided to the linearised SPM, and the terminal voltage was calculated by solving the model in the time domain. This requires solution of (29-31, 49, 54), which requires knowledge of $Q_{i}^{th}$, the theoretical capacity of each electrode, which relates variations in stoichiometry $x$ to variations in measured discharge capacity. Unfortunately $Q_{i}^{th}$ is unknown, but by introducing two changes of variable, $\hat{u}_{i}=Q_{i}^{th}\bar{u}_{i}$ and $\hat{x}_{i}^{s}=Q_{i}^{th}\bar{x}_{i}^{s}$, this term can be eliminated from the linearised SPM, similarly to section III-C. This approach assumes the model is linearised; to run the SPM with non-linear OCV would require knowledge of $Q_{i}^{th}$. The time-domain linearised SPM was solved in MATLAB® using Chebyshev orthogonal collocation as in our previous work [1, 45].

The comparison between linearised SPM and measured data is shown in Fig. 9 over a time period of almost 10 minutes. The DoD at the beginning of the period was $10\text{\,}\mathrm{\char 37\relax}$ and at the end of the period was $20\text{\,}\mathrm{\char 37\relax}$. The OCV curves are almost linear over this range (see Fig. 7) and the variation in $R_{0}$ is small (Fig. 6). Initial results using this approach showed a consistent increasingly negative offset in voltage predictions over the time period examined. This is most likely due to the assumed OCV slope at this DoD being inaccurate, and was adjusted by multiplying the linearised cathode OCV slope by a factor of 0.78. Fig. 9 shows both the initial simulated and the adjusted simulated data in comparison to the measured data. The voltage errors between adjusted simulated vs. measured data are maximum $20\text{\,}\mathrm{mV}$ and RMS $10\text{\,}\mathrm{mV}$. These are satisfyingly reasonable given the fact the model was parametrised completely separately with frequency domain data, that different cells were used for OCV measurements compared to EIS and time domain data, and that the kinetics were linearised despite current peaks of up to 6C.

VI Summary and conclusion

The analysis of parameter identifiability from experimental data is crucial prior to any attempt at estimating the parameters of battery first-principle models. We have demonstrated that the lithium-ion battery single particle model is over-parametrised and that only six subgroups of parameters are necessary for full parametrisation. Assuming that the OCV vs. discharge capacity curves for each electrode are known, we have shown that the estimation of the linearised model parameters at a given DoD can only identify three parameters among these six subgroups of parameters, namely the cathode diffusion time constant $\tau_{+}^{d}$, the anode diffusion time constant $\tau_{-}^{d}$ and a charge-transfer resistance $R_{ct}$. Finally, we have shown the crucial role of the slope of electrode open-circuit voltage curves on the parameter uncertainty. This has been investigated through analysis of the model, as well as through fitting both time and frequency domain data. A flat electrode OCV curve results in high uncertainties on the electrode diffusion parameter. Therefore impedance data at a single DoD cannot, in general, result in accurate parameter estimation for both electrodes, and experimental data at multiple DoDs must be considered. Specifically, complimentary DoDs where in turn the anode and cathode OCV gradients are significant must be chosen to yield satisfactory identifiability. Future work could investigate the effects of temperature on parameter estimation, and the evolution of cell parameters as a battery ages.

Acknowledgment

Financial support is gratefully acknowledged from EPSRC UK (EP/K002252/1) and Samsung Electronics Co. Ltd. We thank Dr Adam Mahdi (University of Oxford) for his valuable advice on model structural identifiability analysis, and Dr Christoph Birkl and Dr Shi Zhao (University of Oxford) for providing the experimental OCV and dynamic test data.