Thermal uncertainty analysis of a single particle model for a Lithium-Ion cell

In recent years, the sustained increase in the worldwide demand for lithium batteries goes side by side with the need for reliable methods to assess battery performance. Of particular importance is assessing Lithium-Ion cell’s thermal behavior given its role on hazard and aging of batteries. An essential task towards developing battery management systems is the estimation of physical parameters and their uncertainties in terms of both models and observations. Consequently, this paper analyzes the uncertainty in a thermal single particle model for a lithium-ion cell. In the first part of the manuscript, we explore the model adequacy by analyzing the forward and backward uncertainty propagation in the model in terms of diffusion coefficients, reaction rate constants, and observations of cell voltage. In the second part of the manuscript, we infer the cell’s reversible heat given the energy balance equation and the cell’s temperature observations. We argue that the methods proposed here may be extended to analyze other more general lithium-ion models.


Introduction
The worldwide demand for Lithium-Ion batteries has continuously increased during the last decades. The development of models to simulate battery behavior has gone side by Dedicated  side. Broadly speaking, first principle models of batteries are based on multiphysics and are fundamental from a theoretical point of view [1][2][3]. Unfortunately, these models are complex and computationally expensive to simulate, which is a limitation in applications, i.e., design and analysis, where the intense simulation of battery behavior under multiple regimes is required. Because of these arguments, the worldwide community has developed reducedorder models. These types of models, many of them phenomenological, reproduce battery dynamics successfully.
Computer models arising from the discretization of battery models will be useful as long as there is confidence in their predictions, mainly to assist decision making [4]. Since reducedorder models are typically easier to implement and faster to simulate, quantitative methods are required to evaluate reduced-order models' ability to predict quantities of interest as input parameters vary across regimes. An appropriate approach to address this question is the forward and backward uncertainty propagation analysis. In this regard, it is important to analyze the model adequacy to predict the lithium-ion cell thermal behavior given its role on battery performance, safety, aging, wear and tear.
To address these questions, in this paper, we consider a thermal single-particle model (TSPM), studied by Guo et al. [5], that assumes non-constant temperature along time and couples lithium concentration and thermal balance equations. It is a reduced-order model whose numerical simulation cost allows it to combine with Monte Carlo simulations.
First, we analyze the forward uncertainty propagation of the transport and kinetic coefficients-i.e., its identifiability-through global sensitivity analysis. To accomplish this task, we analyze each coefficient's contributions to the likelihood's variance. We define the likelihood as the conditional probability of voltage observations given these model coefficients in the present setting. We find that four parameters (the diffusion coefficient and the reaction rate activation energies, of both electrodes) with the same scale are identifiable from voltage measurements.
On the other hand, to carry out backward uncertainty analysis, we elicit prior distribution models for the activation energy parameters and use them, with the likelihood, to analyze the conditional probability of these coefficients given voltage observations. Of note, voltage observations are easier to obtain compared to temperature observations, see [6]. We apply this strategy to Guo's reduced-order model to find a conditional probability of the activation energy parameters given voltage observations. We remark that this is not a trivial task since "...the task of parameter identification is challenging. [...] many parameters are weakly identifiable from current and voltage measurements....", see [7]. In this context, we must highlight that we make inferences with a computer model whose solutions mimic the continuous model solutions with a specific order of error in practice. In this paper, we use the estimate introduced in [8] for the solution of the inverse problem in terms of the computer model error. Based on voltage observations, this estimation provides evidence to decide on the grid thickness used to solve the computer model.
As a second line of evidence of the TSPM's predictive capacity, we infer the reversible heat of the lithium-ion cell using the energy balance equation and temperature observations during a discharge cycle. A central part of our method is estimating the temperature derivative solving a linear inverse problem with the Bayesian paradigm. We use spatial statistics methods to propose a prior model of the temperature derivative, which we regard as a distributed parameter, using a Gaussian Markov random field. We accomplish this task in two steps. First, we pose a uniformly spaced prior model, and secondly, we subsample a central subinterval where the pointwise variance is smaller. This strategy allows us to reduce the variance of the posterior model and reduce the number of degrees of freedom. This approach is based on well known model reduction methods, see for example Spantini et al. [9] or Bui et al. [10].

Related work
There are many research efforts aimed at analyzing cycling performance of lithium ion batteries taking into account thermal effects [5,[11][12][13][14]. Marcicki et al. [15] considered a lumped energy balance equation and Kalman filtering to estimate heat generation during cycling for Li-ion cells. Bizeray thesis [16] contributes towards state and parameter estimation of thermal-electrochemical lithium-ion battery models, noteworthy are the performance estimation analyses based on extended Kalman filtering. Tagade et al. [17] carried out a synthetic analysis of experimental design and bayesian calibration of an electrochemical termal model. Chen et al. [18] proposed a 3D model aimed at exploring battery thermal behavior taking into account the battery shape. Raijmakers et al. [6] provided a review on heat generation principles, as well as methods and challenges in battery temperature measurement methods.

Contributions and limitations
We explore the adequacy of a thermal single particle model for a Lithium-Ion cell. Our forward uncertainty propagation analysis reveals a relation between forward uncertainty propagation and practical parameter identifiability when we consider the likelihood function as a quantity of interest (QoI). Moreover, in the backward uncertainty analysis we find a probability distribution for these parameters. On the other hand, we estimate the reversible heat using this model for a lithium-ion cell with quantified uncertainty. We use methods of spatial statistics to elicit a prior statistical model for the derivative of the temperature. The corresponding inference process accounts for dimension reduction exploiting the fact that the conditional probability of the temperature derivative is a Gaussian distribution. We point out that our analysis can be applied to other models for lithium-ion cells. We show the importance of prior model elicitation to estimate the reversible heat.
The remainder of this manuscript is organized as follows: Sect. 2 presents the analyzed model and the theoretical analysis tools. Section 3 makes a discussion of our results. Finally, Sect. 4 summarizes our findings and their generality.

Methods
For the sake of making this paper self-contained, in this section, we gather the theoretical methods used in our in silico analysis.

Thermal Li-ion single cell model
We consider the thermal single particle model for a Lithium-Ion cell described by Guo et al. [5], which reads as follows. Assuming each electrode can be represented by a single intercalation spherical particle, the mass balance of lithium ions in an intercalation particle of electrode active material is described by Fick's second law in a spherical coordinate system where c s,j is the concentration of lithium ions in the intercalation particle, t is time, r is the radius direction coordinate, D s,j is the solid phase diffusion coefficient, which is a function of temperature, and the subscript j = p/n represents the positive/negative electrode. We also have the initial conditions and the boundary conditions where r = R j is the particle radius, and with I the total current, F = 96,487 C/mol the Faraday constant and S j the total electroactive surface area of electrode j.
Assuming that the spatial temperature distribution in the cell can be neglected, the cell temperature T is a function of time only, satisfying Here M and C p are the mass and the specific heat capacity of cell, respectively. In addition, U j is the Open Circuit Potencial (OCP) and x j,surf = c s,j| r=R j c s,j,max is the State of Charge (SoC). The overpotential η j is defined as being k j the reaction rate constant and c e the (assumed constant in this model) electrolyte concentration in solution phase. Finally, R cell is the electrolyte resistance and q is the rate of heat transfer between cell and surroundings. The initial cell temperature is assumed to be ambient temperature In this context, it is possible to explicit the cell voltage V cell as The transport and kinetic parameters involved in the model equations are functions of temperature. The solid phase diffusion coefficients D s,j and the reaction rate constants k j depend on the temperature through Arrhenius' correlation where E a di,j and E a re,j are the diffusion coefficient activation energy and the reaction rate activation energy, respectively, of electrode j. Of note, the dependence of coefficients D s,j on temperature implies that Eqs. (1) and (4) are coupled.

Computer model
Broadly speaking, battery model analysis is concerned with studying the solution of model (1)-(6) under different parameter regimes. Let us denote the activation energy parameter vector and let V cell (t) be the cell voltage at time t. Equations (1)-(6) define a forward mapping M from these model parameters to cell voltage The forward mapping V cell = M(θ ) has no analytical expression. Thus, to explore problem (1)-(6), we rely on a computer model where h accounts for a discretization size and N is the dimension of the corresponding finite dimensional space. In this paper, we propose M h implementing second-order centered finite differences in space and Crank-Nicolson scheme to march in time to deal with equation (1), as well as a two stage Runge-Kutta second-order method (namely, Heun's method) to approximate the solution of thermal equilibrium Eq. (4).
Practical validation of a computer model is necessary if we are to use model M h as a surrogate of model M to explore backward and forward uncertainty propagation. In the absence of an analytical expression of V cell , in order to provide a numerical convergence rate of the computer model, we compute where the norms ||·|| ∞ are computed on the corresponding coarser mesh. Of note p h ≈ 1.99, for h = 2 − j for j = 2, 3, 4, . . . This value of p h suggests a quadratic order of convergence.

Forward and backward uncertainty analysis
for a given set of parameters θ with additive noise, e.g.
In the synthetic examples shown in Sect. 3, the standard deviation σ of the noise is set such that the signal to noise ratio (SNR) satisfies to mimic the precision of instruments used to measure voltage. We have made the hypothesis that there is no interference, see Johnson [26]. If we assume observations are independent then we obtain a model of the conditional probability of v given θ In Sect. 3 we analyze the uncertainty in the forward model V cell = M(θ ) through global sensitivity analysis. To accomplish this task, we explore the practical parameter identifiability by taking the likelihood (8) as a quantity of interest in the Sobol sensitivity analysis.
According to Salltelli et al. [27], the first order Sobol sensitivity analysis for a quantity of interest Y = Y (θ ) is defined by where var(Y ) is the variance of the quantity of interest Y , var θ i (·) denotes the variance of argument (·) taken over θ i , and E θ ∼i (·) denotes the expectation of argument (·) taken over all factors but θ i . We use (9) to explore the contribution of the variance in each model parameter to the variance in the likelihood (8). In Sect. 3 we argue that all activation energy parameters θ are practically identifiable given voltage observations. On the other hand, in the backward uncertainty analysis, we care about the uncertainty in the parameters θ given observations of the voltage. Thus, we shall use the Bayesian paradigm to model the conditional probability of the vector of model parameters θ given a vector of codes into a probability density function available information about the parameters that is independent from the voltage observations. If π (θ ) is a prior distribution of the model parameters, then the uncertainty of the model parameters given observations of the voltage is given by the posterior distribution where Z (v) is the normalization constant, often unknown. There are no analytical methods to explore the posterior distribution (11), e.g., to obtain its estimators, for example, the maximum a posteriori θ M AP = max θ π |V (θ |v), the posterior conditional mean θ C M = E[π |V (θ |v)] or posterior median θ median . Thus, in Sect. 3 we shall use a stochastic collocation method called Markov Chain Monte Carlo (MCMC) to approximate numerically the posterior distribution (11).

Elicitation of the prior distribution for the activation energy parameters
To elicit a prior distribution (10) we shall assume E a di,j , E a re,j ∈ [0, a × 10 5 ], j=p/n. We set a = 2 for E a di,p and a = 1 for the rest of the activation energy parameters. Also, we use the fact that the parameter correlation structure is implicitly set into the forward model V cell = M(θ ) to pose the prior model as the product of a prior model for each parameter. Consequently, we model E a di,j /(a × 10 5 ), E a re,j /(a × 10 5 ) ∼ Beta(α, β), where instrumental parameters α, β are chosen following E[log(θ i )] = ψ(α) − ψ(α + β), var(log(θ i )) = ψ 1 (α) − ψ 1 (α + β), being ψ and ψ 1 are the first and second logarithmic derivatives of the Gamma function.

Bayesian approach
We shall use the equation of energy balance (4) to estimate the reversible heat given temperature observations T (t i ) at times t i for i = 1, . . . , m. Following the notation of Marcicki [15], we assume the irreversible part of the generated heat to be q irr = I η p − η n + I R cell , while Newton's law of cooling holds for the heat rejected to the environment q ∞ = h (T − T amb ). Furthermore, we shall assume that the cell mass M and the cell specific heat capacity C p are known. Consequently, if an approximation of the temperature derivative dT dt is available, then we can use the energy balance (4) to estimate the reversible heat q rev by In order to use Eq. (12), we propose a Bayesian approach to estimate the temperature numerical derivative given noisy observations of T . If we regard numerical differentiation as an inverse problem, then the forward problem is a convolution with a Heaviside function H where f is the derivative of T in the zero noise limit and γ 2 T is the noise variance. In the results presented in Sect. 3 we set the signal to noise ratio We shall use Eq. (13) to pose both, a computer model to solve the direct problem and a likelihood model. We discretize the integral in Eq. (13) using the composite medium point rule obtaining where we have assumed uniform timing in the observations t i = i s for i = 0, . . . , m and s j = 1 2 (t j−1 + t j ) for j = 1, . . . , m.
, the discretized inverse problem (14) can be writen as where A = (a i j ) m i, j=1 is the lower triangular matrix defined by Assuming temperature observations are independent, the likelihood is given by the conditional probability where obs = γ 2 T I m and I m ∈ R m×m are the covariance and identity matrices respectively. We shall solve the inverse problem in two stages. First, we pose a Gaussian prior model for the probability distribution of the temperature derivative with zero mean and covariance matrix pr . Here, we define the covariance matrix in terms of the precision matrix −1 pr = γ −2 L T L, where γ 2 is the marginal variance and matrix L ∈ R m×m is a model proposed by Bui [29] as follows. We let For the end points f 1 and f m we postulate Then γ −2 L T L = −1 pr is the precision matrix where The marginal variance γ 2 is chosen to be a function of the number of degrees of freedom, i.e. γ 2 = γ 0 m 3 for some γ 0 > 0 and m is the number of degrees of freedom. Equations (16) and (17), together with Bayes rule (see [30]), define a Gaussian model for the conditional probability of the derivative of the temperature f given a vector of temperature measurements T with explicit expressions for the posterior covariance matrix (19) and the posterior mean vector The second stage of the inverse problem is carried out through dimension reduction in the number of degrees of freedom.

Dimension reduction
It is known that in the Bayesian approach to inverse problems, the update from prior distribution (17) to posterior distribution (18) occurs along a low-dimensional subspace of the space where the distributed parameter f lies. There are research efforts aimed at understanding, computing, and exploiting those data-informed subspaces in the parameter space, see [9,31,32], for both linear and nonlinear inverse problems. Among the data-informed subspace analysis methods, parameter dimension reduction is of paramount importance. Indeed, there is an interplay between smoothness of the forward model (13), the dimension of the parameter space, and the signal to noise ratio that determines the uncertainty in the posterior distribution (18), and this dimension is a feature we can change to reduce the uncertainty in the posterior distribution. In the context of methods for parameter space model reduction, there are approaches based on optimization [10] and approaches based on marginalization [32] among other techniques. It is safe to say that although there are no general parameter dimension reduction methods, most methods aim to enhance the predictive capability of the posterior distribution as a statistical model, e.g., to quantify the uncertainty in the solution of the inverse problem better.
In the present setting, we reduce the dimension of the distributed parameter f , inferred as shown in the Eqs. (18), (20), by computing the pointwise posterior variance and undersampling f at points with smaller variance. To accomplish this task, we use matrix A in Eq. (15) and a linear interpolation f = Bg, where B ∈ R m×n , with n < m, is a linear transformation aimed at collocating the degrees of freedom used to represent f . In the present setting, we define B such that Bg it leaves unchanged r ∈ N degrees of freedom near the endpoints, and undersamples the remainder degrees of freedom every other point.
where I r ∈ R r ×r is the identity matrix, 0 are zero matrices and Global sensitivity analysis is a measure of practical identifiability. Complementary to the analysis of the posterior distribution shown in Fig. 2, the marginal posterior distribution is well defined for parameters with a significant Sobol index given a data set is the interpolation matrix, with n = (m + 1)/2 if m is odd m/2 i fm is even .
This collocation strategy reduces both the posterior variance, and the number of degrees of freedom in the distributed parameter g, as we will see in Sect. 3. The forward model (15) is transformed into Next, we pose a posterior model on g following the derivation of the likelihood (16), as well as the prior model (17) by means of the corresponding marginal varianceγ and matrixL ∈ R n×n to define the covariance matrix in order to form the precision matrixˆ −1 pr =γ −2L TL . Thus obtaining whereˆ post = (AB) T −1 obs AB +ˆ −1 and

Forward and backward uncertainty analysis
To establish the adequacy of the thermal Li-Ion cell model (1)-(6), we start from synthetic simulated data created with 41,000 time steps. Specifically, we solve the direct problem with the values specified in [5] for an initial temperature of 298 K and a discharge rate of 1C. The time required to reach the cutoff voltage is 4100 s. We have used the forward and backward uncertainty analysis described in Sect. 2.3. First, by carrying out a Sobol sensitivity analysis we examine the contribution of each one of the model parameters θ , defined by equation (7), to the likelihood (8), which is the conditional probability of a set of voltage observations given these parameters. Table 1 offers evidence that all activation energy parameters have a meaningful contribution to the variance of the likelihood. We argue that this is a measure of practical identifiability.  Table 1, the posterior distribution is narrower for parameters with larger Sobol sensitivity index.
Next, we use MCMC to sample iteratively the posterior distribution (11) by means of the t-walk from Christen and Fox [33]. Two standard tests provide numerical evidence of the MCMC convergence, see Roberts and Rosenthal [34]. Namely, the trace plot of − log(π V | (θ |v)π (θ )) versus the iteration number shown in Fig. 1, and the computed integrated autocorrelation time (IAT) divided by the number of parameters, which in our case is given by I AT /4 = 13.4.  The posterior distribution gives useful information for the identification of model parameters. In particular, by means of the MCMC samples is possible to compute estimators of the posterior distribution (for example, the median θ median shown in Fig. 2) with a quantified uncertainty.
Also, this MCMC samples allows us to analyze the relation between forward uncertainty analysis and practical parameter identifiability. Indeed, from Table 1 and Fig. 2 we see that posterior variance is smaller for those parameters with larger sensitivity. This outcome is expected given that carry out a global sensitivity analysis (Sobol) of the likelihood over the whole parameter support. Consequently, a measure of practical identifiability given a data set, is that the Sobol sensitivity index is roughly the same order of magnitude for all parameters. This result should be general given a data set, a likelihood and a parameter support. Figure 3a and b show, for temperature and voltage respectively, the true value as well as the values predicted by the posterior median and by the 5-95% percentiles of the posterior distribution. We note that predictions are achieved with very low forward uncertainty.
The inference results shown in Figs. 2 and 3 are obtained by simulating the direct problem with 4100 time steps, while synthetic data, as described in Sect. 2.3, was created using 41000 time steps. Indeed, numerical simulations of the direct problem in computer MacBook Pro (15-inch, 2018), 2.6 GHz Intel Core i7 of six cores, with 16 GB 2400 MHz DDR4 take 1.04 s with 4100 points and 5.92 s with 41000 points. Namely, simulating the posterior distribution with the coarse model saves more than 80% computation time. We argue that this computing efficiency is of paramount importance if we are to sample the posterior distribution of the activation energy parameters θ via MCMC.

Estimation of the reversible heat
Information about the estimation of the temperature derivative for a coarse-grained computer model (4100 time steps) can be found in Fig. 4a. Figure 4b shows the temperature distribution predicted by multiplying by matrix A, the discrete direct operator defined in Sect.  obtained after computing the pointwise posterior variance and undersample the temperature derivative at the points with smaller variance. For this task, we allocate 10 % of the total initial number of degrees of freedom uniformly placed in a neighborhood of the endpoints, as described in Sect. 2.5.2. Of note, the dimension reduction improves the estimation of the reversible heat, and reduces the posterior variance.
The same analysis is presented in Figs. 6 and 7 for a fine-grained computer model (41000 time steps). This simulation renders roughly the same result as the coarse-grained one. Namely, we obtain roughly the same estimate of the reversible heat and the same uncertainty and predictions of both temperature and voltage. Consequently, it may be preferable to use the coarse-grained model if we have to analyze multiple datasets.
We argue that our inference is robust to the numerical discretization, and subsampling the distributed parameters in the region where the pointwise variance is smaller reduces the variance in the posterior distribution.

Conclusions
In this paper, we apply forward and backward uncertainty analysis to assess the adequacy of a thermal single-particle Li-Ion cell model. We have used parameter sensitivity analysis as diagnostics preceding parameter estimation. We offer numerical evidence that forward uncertainty analysis may serve as a practical identifiability test given a set of voltage observations. Afterward, we use error control of the posterior distribution concerning the numerical approximation of the direct problem to reduce the computational burden in approximating such a posterior distribution in order to identify the activation energy parameters via Markov Chain Monte Carlo. In the second part of the manuscript, we infer the reversible heat of a Li-Ion cell given observations of the temperature and an energy balance equation. To accomplish this task, we exploit the linearity of the energy balance equation to construct the solution of the inference problem as a conditional probability given by a Gaussian distribution with explicit expressions for the mean and the covariance, Eqs. (18)- (20). Furthermore, we exploit the fact that the solution grid is different from the data grid to apply a variance reduction method and obtain a posterior distribution with reduced uncertainty, and a smaller number of degrees of freedom, Eqs. (21)- (22), and Algorithm 1.
Code and data availability For the sake of reproducibility, all computations shown in this section were carried out using Anaconda Python [35] fixing the seed of random number generation as follows: numpy.random.seed(2021). Code can be obtained from the corresponding author upon request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.