On-board monitoring of 2-D spatially-resolved temperatures in cylindrical lithium-ion batteries: Part II. State estimation via impedance-based temperature sensing

Impedance-based temperature detection (ITD) is a promising approach for rapid estimation of internal cell temperature based on the correlation between temperature and electrochemical impedance. Previously, ITD was used as part of an Extended Kalman Filter (EKF) state-estimator in conjunction with a thermal model to enable estimation of the 1-D temperature distribution of a cylindrical lithium-ion battery. Here, we extend this method to enable estimation of the 2-D temperature field of a battery with temperature gradients in both the radial and axial directions. An EKF using a parameterised 2-D spectral-Galerkin model with ITD measurement input (the imaginary part of the impedance at 215 Hz) is shown to accurately predict the core temperature and multiple surface temperatures of a 32113 LiFePO$_4$ cell, using current excitation profiles based on an Artemis HEV drive cycle. The method is validated experimentally on a cell fitted with a heat sink and asymmetrically cooled via forced air convection. A novel approach to impedance-temperature calibration is also presented, which uses data from a single drive cycle, rather than measurements at multiple uniform cell temperatures as in previous studies. This greatly reduces the time required for calibration, since it overcomes the need for repeated cell thermal equalization.


Introduction
Monitoring the temperature of Li-ion batteries during operation is critical for safety and control purposes. The conventional approach to temperature estimation is to use numerical electrical-thermal models coupled with online measurements of the cell surface temperature and/or the temperature of the heat transfer medium [1] (Figure 1)a). 1
Using this approach in conjunction with state estimation techniques such as Kalman filtering, the cell internal temperature may be estimated with high accuracy [2,3,4,5]. However, large battery packs may contain several thousand cells [6], and so the requirement for surface temperature sensors on every cell represents substantial instrumentation cost. As a result, EV manufacturers often use fewer temperature sensors than are required to achieve full observability of the pack [7]. Moreover, rapid fluctuations in internal temperature may not be registered by surface mounted temperature sensors, regardless of the sampling frequency. This may mean thermal runaway cannot be detected, since associated timescales are often shorter than those associated with heat conduction through the cell [8]. Consequently, surface mounted temperature sensors even when used with a thermal model may be insufficient to track internal temperature or predict thermal runaway. One approach to overcome these problems is to embed flexible thin film micro-temperature sensors within the cell to enable in-situ internal temperature measurement [9,10,11,12,13]. Whilst this has some obvious advantages, the additional manufacturing and instrumentation requirements would significantly increase the cost and complexity of the system.
An alternative approach to temperature estimation uses electrochemical impedance spectroscopy (EIS) measure- ments at one or several frequencies to directly infer the internal cell temperature [14,15,16,17,18,19,20,21,22] ( Figure 1)b). This exploits the fact that impedance is a function of the cell internal temperature. For brevity, we refer to the use of impedance to infer information about the internal temperature as impedance-temperature detection (ITD). ITD has promise for practical applications, since methods capable of measuring EIS spectra using existing power electronics in a vehicle or other applications have been developed [23]. Most ITD studies assume the internal temperature is uniform [14,15], or explicitly acknowledge that the impedance measurement alone only predicts the average internal temperature [16]. Our recent work showed that ITD could be used in conjunction with a thermal model in a 'hybrid' ITD-based state estimation scheme, to enable estimation of the temperature distribution of the cell [19] (Figure 1c). The method is therefore analogous to the conventional state-estimation method, but the measurement input consists of ITD rather than a surface temperature measurement. However, that work relied on the assumption of 1-D conditions, which are unlikely to arise in all in situations.
In the present paper, we extend the ITD approach to problems involving 2-D thermal dynamics, using the loworder thermal model presented in Part I .

Impedance-based temperature sensing
The electrochemical impedance Z(ω) = Z (ω) + jZ (ω) of lithium-ion cells is a function of temperature, state of charge (SoC), and state of health (SoH). Within an appropriate frequency range, however, the dependence on SoC and SoH is negligible and the impedance can thus be used to infer information about the cell temperature. Previous ITD studies have used as a temperature-dependent parameter (TDP) the real part of the impedance at a specific frequency [19,16], the imaginary part of the impedance at a specific frequency [21], the phase shift at a specific frequency [14,15], and the intercept frequency [18]. The issue of which TDP and excitation frequency are most suitable for temperature inference is still an open question, with various arguments in favour of each [22,24]. In the present study, we use the imaginary part of the impedance at f = 215 Hz as the TDP, since this was found to give superior results to the real part for the studied cell. However, in principle the method is not limited to this option and could also be applied using any TDP.
The principle of operation of ITD is to relate the cell impedance to the mean 2 cell temperature by a polynomial fit: where a 1 , a 2 and a 3 are the constant polynomial coefficients. In our previous work the coefficients of the polynomial fit were identified via offline impedance measurements at multiple uniform temperatures [19]. However, here we introduce an improved calibration technique, which uses data from a single drive cycle (see Section 6.3). This greatly reduces the time required for calibration, since it overcomes the need for cell thermal equalization at multiple temperatures.
Finally, the principle of operation of the hybrid ITDbased state-estimation method is to use ITD as the measurement input to a state-estimator in conjunction with the 2-D thermal model presented in the Part I of this paper [26].

Overview of hybrid method
An overview of the process is shown in Figure 2 and described below.
Offline: The cell dimensions (jelly-roll inner radius, r in ; outer radius, r out ; and height, H) are measured and therefore known a-priori; the thermal properties (density, ρ; specific heat capacity, c p ; and radial conductivity, k r ) are available from the literature [27], and therefore known apriori. The remaining thermal properties (axial conductivity, k z ; and the convection coefficients on the left end, h l , right end, h r , and curved surface, h t ) are then estimated using a recursive least squares fitting algorithm (see Section 6.2). Specifically, the thermal model is simulated using experimental current and voltage data (V exp , I exp ), from a parameterization drive cycle and the temperatures at four locations (T 1 , T 2 , T 3 , and T 4 ) are compared with their corresponding thermocouple measurements. The parameterised model is then used to fit the impedance coefficients (see Section 6.3). Specifically, the predicted mean  temperature (T ) from the model is paired with the corresponding measured impedance value (Z ) at each sample. A second order polynomial fit is then applied to the resulting impedance-temperature data, thus identifying coefficients (a 1 , a 2 , a 3 ). Finally, the known and identified thermal parameters and the identified impedance coefficients are provided to the online thermal model.
Online: The thermal model (Section 4.2) uses online measurements of the voltage and current (V , I) to predict the heat generation, 2-D temperature distribution, and cell imaginary impedance at each time step. This is compared to the measured imaginary impedance which is used to update the state estimate of the model via an Extended Kalman Filter algorithm (EKF) (Section 4.3).

Problem definition
The model consists of the transient 2-D energy conservation equation in cylindrical coordinates. Heat generation is assumed to be uniform in space but time-dependant. The multi-layer structure of the battery is treated as a homogeneous solid with anisotropic thermal conductivity in the radial and axial directions. The temperature variation in the azimuthal (ϕ) direction is neglected. Convective heat transfer is assumed to occur at the outside surfaces, and there is zero convection at the inner annulus of the jelly roll (i.e. at the central mandrel). The convection coefficient of the ambient air may be different at each surface, although its free-stream temperature is assumed to be uniform throughout the chamber. A schematic of the model is shown in Figure 3.
The model is governed by the following 2-D boundary value problem [28]: where t is time and r and z are the position coordinates in the radial and axial directions respectively. The functions T (r, z, t) and q(t) are the temperature distribution and volumetric heat generation rate, respectively. The parameters ρ and c p are the density and specific heat capacity respectively, and k r and k z are the anisotropic thermal conductivities in the r and z directions. The boundary conditions are given by: where T ∞ is the free-stream air temperature of the chamber, and {h σ ; σ = t, r and l} are the convection coeffi-cients at the top, right and left surfaces. Note that the convection coefficient at the bottom surface is set to zero (since this corresponds to the inner radius of the cell jelly roll, which is not exposed to cooling). h t corresponds to the curved surface of the cell, whilst h l and h r correspond to left and right ends of the cell. Note that the placement of the heat sink (see later) results in an increased value of h l . As in our previous work [19], we consider only ohmic heat generation, given by: The entropic heat is neglected because (i) the net reversible heat would be close to zero when the cell is operating in HEV mode and (ii) the HEV drive cycles employed in this study operate the cell within a small range of SoC (47 − 63 %) and hence the entropic heat is small [1]. The heat generation is assumed to be uniformly distributed throughout the cell volume, hence the volumetric heat generation is given by: where V b is the cell volume.

Spectral-Galerkin model
The above problem can be simulated using the spectral-Galerkin approach described in Part I of this paper. The model can be expressed in state space form: where the states and the state matrices are as defined in Part I, and the model input is u = [q(t), 1] T . Four temperature outputs are chosen corresponding to the positions T 1 = T (r = r in , z = H/2), T 2 = T (r = r out , z = 0), T 3 = T (r = r out , z = H/2) and T 4 = T (r = r out , z = H) (see Figure 3). These outputs were chosen to match the thermocouple measurements in the experimental setup as described in the following section.
Lastly, the cell imaginary impedance at the selected frequency is also calculated as an output. This is obtained from the mean temperature using eq. 1. T is itself a function of the model states, x as described in the companion paper. Hence an expression for the impedance as a function of the cell states is obtained, The computed impedance is used as the measurement input in the state estimation algorithm described next.

State estimation
The state estimation consists of an extended Kalman filter (EKF), for estimating the temperatures at each of the four thermocouple locations, with the cell impedance as measurement input. The EKF is equivalent to that employed in [19] for the 1-D case, except that here it is applied to the 2-D SG model.
We firstly modify eq. (6) by rewriting it as an explicit state-space model in discrete time: whereĀ andB are system matrices in the discrete-time domain, given bȳ where ∆t is the sampling time of 1 s. We then set the impedance as the model output where f (x k ) is the non-linear function relating the state vector to the impedance measurement (i.e. eq. 8), and v k and n k are the process and measurement noise, respectively. Their corresponding covariance matrices are R v and R n . Note that, although the impedance is the model output for the purpose of the EKF algorithm, the temperatures at each of the four thermocouple locations are also computed from the identified states at each time step, for validation against the thermocouple measurements.
The time update processes are then given by: wherex − k andx k are the a priori and a posteriori estimates of the state, and (P x k ) − and P x k−1 are the corresponding error covariances. Since the relationship between impedance and the cell state is non-linear, the measurement model must be linearised about the predicted observation at each measurement. The measurement update equations are: where K x k is the Kalman gain for the state, and H x k is the Jacobian matrix of partial derivatives of f with respect to x: The above algorithm can be simplified to a standard Kalman Filter (KF) by omitting the linearisation step (eq. 18) and replacing f (·) in eqs. 12 and 16 with an appropriate linear operator. Such a KF is applied later using the surface temperature, T 3 , as measurement input for comparison with the EKF using Z .

Setup
Experiments were carried out with a 4.4 Ah cylindrical cell (A123 Model AHR 32113 Ultra-B) with LiFePO 4 positive electrode and a graphite negative electrode. A large form factor cell was used in order to ensure measurable 2-D effects. The properties of the cell are given in Table 1. The cell was fitted with four thermocouples, three on the surface and another inserted into the core via a hole which was drilled in the positive electrode end. The thermocouple locations correspond to the model output locations described previously. Cell cycling and impedance measurements were carried out using a Biologic HCP-1005 potentiostat/booster. The impedance was measured using Galvanostatic Impedance Spectroscopy (GEIS) with a 200 mA peak-to-peak perturbation current. The environmental temperature was controlled with a Votsch VT4002 thermal chamber. The chamber includes a fan which runs continuously at a fixed speed during operation. Photos of the test equipment and a schematic of the experimental setup are shown in Figure 4.
Two different cooling configurations were tested (see Figure 4d). In Config. 1 the right end and entire curved surface of the cell are thermally insulated, and a heat sink fixed to the left end. An additional (auxiliary) fan is placed inside the chamber as shown in Figure 4b. In this case, both the built-in chamber fan and the auxiliary fan run continuously throughout the duration of the experiments. This setup is designed to achieve maximum cooling from the left end of the cell (heat sink) whilst minimising radial heat transfer. In Config. 2, only the right end of the cell is insulated, and the auxiliary fan is switched off. This setup allows for greater radial heat flux. 100 mm Jelly-roll outer radius (r out ) 16 mm Jelly-roll inner radius (r in ) 1 mm

Procedure: Full spectrum impedance
The procedure described here applies to the results presented in Section 6.1. In order to verify that the impedance at the selected frequency was independent of SoC, full spectrum impedance measurements were carried out at a uniform cell temperature of 10 • C, over a range of 10-90% SoC in intervals of 20% 3 . The EIS frequency range spanned 5 kHz to 0.1 Hz with 10 frequencies per decade and 3x averaging. The cell capacities at each temperature were first determined with a constant current constant voltage (CCCV) charge/discharge. The SoC was then adjusted to the required values by drawing a 1 C current. At each SoC, the cell was allowed to rest to ensure its temperature had equilibrated, which typically occurred in less than 2 hrs. This was verified when the temperatures registered by the internal and surface thermocouples were within 0.2 • C and repeated impedance measurements over a 20 minute period yielded identical results.

Procedure: validation experiments
The procedure described here applies to the results presented in Sections 6.2 -6.4. Dynamic experiments were conducted using three different current excitation profiles, denoted HEV-I, HEV-II, and HEV-III. Each cycle was generated by looping over different portions of an Artemis HEV drive cycle, and scaling the applied currents to the range ±50 A. HEV-I was used to parameterise the thermal model, and HEV-II and HEV-III were used to validate the identified parameters and to demonstrate the temperature estimation technique (see Section 6.2).
The procedure for the experiments was as follows: impedance measurements at 215 Hz were carried out every 24 s with 4 s pauses before each impedance measurement, and the four thermocouple temperatures were also monitored. Before each experiment, the SoC was adjusted to 50% by drawing a 1 C current. The temperature of the thermal chamber was set to 8 • C in order to maintain the cell within a typical operating range (cell temperatures of ∼20 to 30 • C are generally optimal [29]). The cell was allowed to rest until its temperature equilibrated before each test began.

Full spectrum impedance
Our earlier work demonstrated independence of the impedance at 215 Hz with respect to SoC for a 26650 LiFePO 4 cell (A123 ANR26650 m1-A) over a range of temperatures from -20 to 45 • C. Since the cell used in the present study is from the same manufacturer and has the same chemistry, it was deemed sufficient to verify SoC independence at a single temperature close to the middle of this range. Hence, the full spectrum impedance tests were carried out at a uniform cell temperature of 10 • C (over a range of 10-90% SoC in intervals of 20%). The results are provided in the supplementary material. They confirm that the impedance at 215 Hz is approximately independent of SoC, as required.

Parameterisation and validation
The full set of thermal parameters include: ρ, c p , k r , k z and the three convection coefficients (one on each end, and one on the cell curved surface), denoted by h l , h r and h t (see Fig. 3). However, three of these parameters were known a priori from the literature: Fleckenstein et al. [27] identified values for the density, ρ, the specific heat capacity, c p , and the radial thermal conductivity, k r , for an identical cell using thermal impedance spectroscopy, and so these values are used in the present case. Hence, parameterisation is only required to estimate the remaining four parameter values: the axial thermal conductivity, k z , and the three convection coefficients: h l , h r and h t .
The measurements from HEV-I (comprising cell current, voltage, and surface and core temperatures, and the chamber temperature) were used for the parameter estimation for both cooling configurations. HEV-I was deliberately chosen for the parameterization since it results in slightly higher cell temperatures than HEV-II or HEV-III, and hence ensures that the polynomial fit is applied using the largest temperature range possible. For Config. 1, all four parameters were identified. For Config. 2, the axial thermal conductivity, k z , was assumed known a-priori, using the identified value from Config. 1 (since k z is the same in each case), and hence it was only necessary to identify the three convection coefficient parameters.
The parameterisation was carried out using fmincon in Matlab to minimise the magnitude of the Euclidean distance between the measured and estimated temperatures for each of the four thermocouples. Concretely, the error between the measured (subscript 'exp') and model predictions (subscript 'm') at each time step, k, is given by,  (19) and the parameters are identified by, where N f is the number of time steps in the cycle. Table 2 presents the thermal parameters, including those known a-priori and those identified via parameterisation. The convection coefficient values are within the range expected of forced convection air cooling [30]. The left coefficient is greatest in both cases as expected due to the presence of the copper heat sink. The values of h r and h l are greater in Config. 1, due to the presence of the auxiliary fan, whereas the value of h t is greater in Config. 2 since the curved surface is uninsulated. The value of k z is ∼ 55 times greater than k r ; this is typical of cylindrical cells with wound jelly-roll constructions [27]. Uncertainty in the parameter estimation may be attributed to manufacturing variability, error in the heat generation calculation (due to the omission of entropic heating), heat generation in the contact resistances between the cell and connecting wires and/or measurement uncertainty in the temperature. The measured core and surface temperatures (subscript 'exp') and the corresponding model predictions (subscript 'm') for the parameterised model for Config. 1 are shown in Fig. 5a. The model with identified parameters was validated against the second current excitation profile, HEV-II (Fig. 5b). For Config. 2, HEV-I was used as the parameterisation cycle and HEV-III for validation; these results are shown in Figures 5c and 5d respectively. Note that, for clarity, T 4 is omitted from these and subsequent plots since it is very close in value to T 3 . The root-mean-square errors (RMSE) in each case are shown in table 3. The errors in the validation tests are only marginally greater than those in the parameterisation tests, indicating that the estimation is satisfactory.

Impedance calibration
The calibration of the impedance-temperature polynomial coefficients (a 1 , a 2 and a 3 ) is achieved as follows. The  current/voltage data from Config. 1, HEV-I is applied to the parameterised model (Fig. 6a. A KF is applied to the model using the surface temperature (T 3 ) as measurement input 4 . The predicted T m output from the KF at each measurement step is paired with the measured Z exp value from the HEV-I experimental data. A second order polynomial fit is then applied to the Z exp vs. T m data (Fig. 6b). Fig. 6 shows that a second order fit is capable of closely approximating the measured data. It should be noted that an Arrhenius fit could also be obtained but the polynomial fit is sufficiently accurate and facilitates faster online computation.

State estimation
We now compare the performance of the EKF with that of a KF based on the same thermal model but with T 3 as the measurement input rather than Z . We chose T 3 for comparison since this is the location of the measurement input used in earlier studies [2,3,19]. The initial state estimate provided to the battery is a uniform temperature distribution at 25 • C, whereas the true initial battery state is a uniform temperature distribution at 8 • C. The covariance matrices are calculated as R n = σ 2 n and R v = β 2 v diag(2, 2). The first tuning parameter for the EKF is chosen as σ n = 3 × 10 −5 Ω, based on the estimated standard deviation of the imaginary impedance measurement. The second tuning parameter was chosen as β v = 5 × 10 −3 , by trial and error. For the KF with T 3 as measurement input the tuning parameters were chosen as σ n = 5 × 10 −4 • C and σ n = 0.05, which are the same values as those used in our previous work [19].
Figs. 7 a)-c) and d)-f) show the results for Config. 1 and Config. 2, respectively. For each configuration, the core and surface temperatures quickly converge to the correct values and are accurately estimated throughout the rest of the excitation profile. Fig. 7b shows that, for the HEV-III cycle, the EKF and KF underestimate the true temperature, in particular during periods with very low or zero applied currents. This error may be a result of a limitation in the thermal model: the heat generation term (eq. (4)) only accounts for ohmic heating, whereas additional electrochemical heating may occur due to the relaxation of concentration gradients, which would result in continued internal heating after the removal of an applied load [29]. This error may be mitigated by the inclusion of an electrochemical heat generation model, although this is beyond the scope of the current work. Figs. 7c) and f) show 2-D contour plots of the cell at the end of the HEV-I cycle and HEV-II cycle respectively. These contour plots are generated by the evaluating the cell temperature on a fine grid throughout the entire (r-z) domain using the thermal model described in Part I. They show that in each case, a significant radial temperature distribution develops within the cell, but in Configuration 2, a larger axial temperature gradient arises, as expected. Fig. 8 shows histogram plots of the errors in the estimates of T 1 and T 3 in each drive cycle for the EKF with Z and the KF with T 3 . Plots a)-d) apply to Config. 1 and plots e)-h) apply to Config. 2. Although the errors may include contributions from various sources (such as measurement uncertainty, variations in ambient temperature or fan flow rate, etc.), we can distinguish between impedance related and model related errors by considering the difference in error arising from the method using impedance measurement input versus surface temperature measurement input. In each case the errors are mostly peaked around a zero mean with a standard deviation in the range 0.2 • C − 0.7 • C. The RMS error in each case is displayed in the top right corner of the plot. It can be seen that the KF with T 3 performs slightly better in most cases, with a lower offset from the zero mean and a narrower distribution. This is to be expected given the higher accuracy of the thermocouple measurement. In particular, the KF performs better in the estimate of T 3 , which is unsurprising since it uses T 3 as a measurement input. However, in general the performance of the EKF is satisfactory and comparable with that of the KF. Fig. 8e shows that the error distribution for the estimate of T 1 for the HEV-III cycle is multi-modal: there are peaks at ∼ −0.2 and ∼ 0.5 • C. This is a manifestation of the errors discussed previously arising from the limitation of the purely ohmic heat generation term.
We note that the RMSE of T 1 in Table 3 is larger than that of the other three sensors (T 2 , T 3 and T 4 ). This is perhaps explained as follows: (i) the parameterization scheme minimizes the error over all four sensors; but (ii) the error in each of the three surface temperatures may be strongly correlated since the temperature gradient in the radial direction is greater than that in the axial; and hence (iii) the parameterization may be biased towards minimizing the error in the three surface temperatures at the expense of the error in the core.
It should also be noted that since the uncertainty of the impedance measurement increases as impedance decreases, the temperature estimates become more uncertain at higher temperatures. Hence, the implementation of this technique could be more challenging at higher ambient temperature conditions than those studied here. This issue was addressed by Spinner et al. [21], who applied a secondary, empirical fit for the upper temperature range to improve accuracy. A similar approach could be applied  using the hybrid method presented here, for applications involving higher temperature.

Conclusions
ITD can be used in conjunction with a thermal model to enable estimation of the spatially-resolved temperature distribution of lithium ion battery cells. The extension of this method to 2-D conditions represents a significant improvement over earlier work, and opens up the possibility of applying this approach to a more general set of conditions than previously possible. For instance, the method could be used to monitor the internal temperature of cells in EV/HEV battery packs, which may have various different cooling configurations. Moreover, it could also be applied using alternative thermal models and/or battery geometries to provide robust estimation of the spatially resolved temperature field of batteries in a range of applications. The calibration of the impedance coefficients using data from a single drive cycle is also a significant improvement over earlier methods. The reduced time and effort required for the calibration step makes practical implementation a feasible goal.
The accuracy of ITD is still slightly inferior to that of conventional methods based on surface thermocouples; however, the reduction in accuracy may be justified by the concurrent reduction in instrumentation cost/complexity. The application of this method online in a vehicle is perhaps the most important area of future work. A recent study has already demonstrated ITD (for average internal temperature estimation) in a vehicle [31]. An interesting area for future work could involve applying the hybrid ITD approach presented here in an on-board setting.
An important consideration for the practical implementation of this approach is the independence of the impedance response to SoH. Whilst such independence has been verified for certain cell chemistries [18], verification for each specific use case would be necessary. Moreover, the accuracy and precision of the impedance measurement will be a crucial factor in on-board applications, and devices capable of achieving high accuracy using low-cost equipment will become increasingly important.