Sensorless Battery Internal Temperature Estimation using a Kalman Filter with Impedance Measurement

4.9 MB·English

by  Robert R. Richardson

1
Sensorless Battery Internal Temperature Estimation
using a Kalman Filter with Impedance Measurement
Robert R. Richardson and David A. Howey, Member, IEEE
Abstract—This study presents a method of estimating battery parameters. However, using additional online measurements -
cellcoreandsurfacetemperatureusingathermalmodelcoupled typicallyofthecellsurfacetemperatureandofthetemperature
withelectricalimpedancemeasurement,ratherthanusingdirect
of the cooling fluid - coupled with state estimation techniques
surface temperature measurements. This is advantageous over
such as Kalman filtering, the cell internal temperature may be
5 previous methods of estimating temperature from impedance,
1 which only estimate the average internal temperature. The estimated with high accuracy [4], [5], [6], [7]. However, large
0 performance of the method is demonstrated experimentally on battery packs may contain several thousand cells [8], and so
2 a 2.3 Ah lithium-ion iron phosphate cell fitted with surface the requirement for surface temperature sensors on every cell
and core thermocouples for validation. An extended Kalman
n represents substantial instrumentation cost.
filter, consisting of a reduced order thermal model coupled
a An alternative approach to temperature estimation uses
with current, voltage and impedance measurements, is shown to
J
accurately predict core and surface temperatures for a current electrochemical impedance spectroscopy (EIS) measurements
5
excitationprofilebasedonavehicledrivecycle.Adualextended at one or several frequencies to directly infer the internal cell
2
Kalman filter (DEKF) based on the same thermal model and temperature, without using a thermal model [9], [10], [11],

impedance measurement input is capable of estimating the [12], [13]. This exploits the fact that impedance is related
]
convection coefficient at the cell surface when the latter is
Y to a type of volume averaged cell temperature, which we
unknown. The performance of the DEKF using impedance as
S the measurement input is comparable to an equivalent dual define later in this article. For brevity, we refer to the use
. Kalman filter using a conventional surface temperature sensor of impedance to infer such a volume-averaged temperature as
s
c as measurement input. ‘Impedance-Temperature Detection’ (ITD). This has promise
[ Index Terms—Lithium-ion battery, impedance, temperature, for practical application, since methods capable of measuring

1 thermal model, Kalman filter, state estimation. EIS spectra using existing power electronics in a vehicle or
v other application have been developed [14], [15], [16]. How-
0 ever, just like conventional surface temperature sensors, ITD
6 I. INTRODUCTION alone does not provide a unique solution for the temperature
1
6 THE sustainable development of transportation relies on distribution within the cell. Our previous work showed that
0 the widespread adoption of electric vehicle (EV) and by combining ITD with surface temperature measurements
. hybrid electric vehicle (HEV) technology. Lithium-ion bat- the internal temperature distribution could be estimated [17].
1
0 teries are suitable for these applications due to their high However, this approach still requires each cell to be fitted
5 specific energy and power density. However, their widespread with a surface temperature sensor. Moreover, whilst the ITD
1 deployment requires reliable on-board battery management technique was validated under constant coolant temperature
v: systems to ensure safe and optimal performance. In partic- conditions, the accuracy of the technique may be reduced if
i ular, accurate on-board estimation of battery temperature is the temperature of the cooling medium is varied rapidly, as
X
of critical importance. Under typical operating conditions, discussed in Section IV. Thus, if the cell electrical/thermal
ar such as a standard vehicle drive cycle, cells may experience properties are known or can be identified, this information
temperature differences between surface and core of 20 ◦C can be exploited to improve the estimate of the thermal state
or more [1]. High battery temperatures could trigger thermal of the cell or reduce the number of sensors required.
runaway resulting in fires, venting and electrolyte leakage. In this study we demonstrate that ITD can be used as the
Whilesuchincidentsarerare[2],consequencesincludecostly measurement input to a thermal model in order to estimate
recalls and potential endangerment of human life. the cell temperature distribution. First, an extended Kalman
The conventional approach to temperature estimation is to filter (EKF) is used to estimate the cell temperature distribu-
use numerical electrical-thermal models [3], [4], [5], [6], [7]. tion when all the relevant thermal parameters, including the
Such models rely on knowledge of the cell thermal proper- convectioncoefficient,areknown.Thethermalmodelconsists
ties, heat generation rates and thermal boundary conditions. ofapolynomialapproximation(PA)tothe1Dcylindricalheat
Models without online sensor feedback are unlikely to work equation. The measurement input consists of the cell current
in practice since their temperature predictions may drift from andvoltage,alongwithperiodicmeasurementsoftherealpart
thetruevaluesduetosmalluncertaintiesinmeasurementsand of the impedance at a single frequency. Second, a dual ex-
tendedKalmanfilter(DEKF)isusedtoidentifytheconvection
The authors are with the Department of Engineering Science, Uni- coefficient online when the latter is not known. The predicted
versity of Oxford, Oxford, UK, e-mail: {robert.richardson, david.howey}
core and surface temperatures in each case are validated
@eng.ox.ac.uk.
Copyright(cid:13)c 2015IEEE against core and surface thermocouple measurements, with
2
agreement to within 0.47 ◦C (2.4% of the core temperature 700
Data
increase). The performance of the combined thermal model 600 Fit: a T2 + a T + a
1 2 3
plus ITD estimator is comparable to the performance of
the same thermal model coupled with conventional surface 500
tienmthpeeractounrteexmteoafsuerxeimstienngtst.emTapbelreatIursehoewstsimtahteiopnretesecnhtnisqtuuedsy, −1ΩY’ / 400
300
which highlights that this is the first study to use impedance
as the measurement input to a thermal model. 200
100
Study Model Measurement 0 10 20 30 40
T (°C)
Tsurf ITD
Forgezetal.[3] (cid:88) (cid:88)
Figure1. Polynomialfittoexperimentaldataofadmittanceatf =215Hz
Linetal.[5],[7] (cid:88) (cid:88) vs.uniformcelltemperature.
Kimetal.[6],[4] (cid:88) (cid:88)
Srinivasanetal.[12],[9] (cid:88)
Raijmakersetal.[11] (cid:88)
Richardsonetal.[17] (cid:88) (cid:88) by offline impedance measurements on the cell at multiple
Presentstudy (cid:88) (cid:88) uniform temperatures [17].
TableI ITD can be viewed as identifying an EIS-based volume
COMPARISONOFONLINETEMPERATUREESTIMATIONTECHNIQUES. average temperature T , which is defined as the uniform
EIS
cell temperature that would give rise to the measured EIS1.
Thus the impedance input is similar to a conventional temper-
ature measurement since it is a scalar function of the internal
II. MEASUREMENTPRINCIPLE temperature distribution but does not uniquely identify the
Theelectrochemicalimpedance,Z(ω)=Z(cid:48)(ω)+jZ(cid:48)(cid:48)(ω), temperature distribution. Either measurement can therefore be
of lithium ion cells is a function of temperature, state of used in conjunction with a thermal model, as shown in Fig.
charge(SOC)andstateofhealth(SOH).Withinanappropriate 2, to estimate core temperature.
frequency range, however, the dependence on SOC and SOH
a b
is negligible and the impedance can thus be used to infer
information about the cell temperature [12]. Previous ITD
studies have used as a temperature-dependent parameter the
real part of the impedance at a specific frequency [10], the
phase shift at a specific frequency [12], [9], and the intercept
frequency [11]. To demonstrate our technique we use the real
part of the impedance at f = 215 Hz. Our previous work
showedthattherealpartofthecelladmittance(theinverseof
thecellimpedance)at215Hzcanberelatedtothetemperature
distributionusingasecondorderpolynomialfit.Foranannular Figure 2. Schematic of (a) the conventional approach to temperature
estimationand(b)theproposedapproachbasedonITD.
cell with inner radius r and outer radius r the real part of
i o
the admittance is given by [17]:
ˆ
Y(cid:48) = 2 ror(cid:0)a +a T(r)+a T2(r)(cid:1)dr (1)
r2 1 2 3 III. THERMAL-IMPEDANCEMODEL
o ri
where a , a and a are the 1st, 2nd and 3rd coefficients of A. Thermal Model
1 2 3
thepolynomialrelatingimpedancetouniformcelltemperature
Thecellthermalmodelconsistsoftheheatequationfor1D
(Y(cid:48) = a + a T + a T2 ), provided that the
1 2 uniform 3 uniform unsteadyheatconductioninacylinder,givenbythefollowing
admittivityvariesintheradialdirectiononly.Thisassumption
Boundary Value Problem (BVP) [2]:
isvalidiftheheattransferfromthetopandbottomendsofthe
cell is negligible, which is approximately true for cylindrical ∂T(r,t) ∂2T(r,t) k ∂T(r,t) Q(t)
cells connected in series with identical cells on either end ρcp ∂t =kt ∂r2 + rt ∂r + V (2a)
b
[18], a configuration which may apply to the majority of
cells in a large battery pack. However, the application of where ρ, cp and kt are the density, specific heat capacity and
this approach to cooling configurations involving substantial thermalconductivityrespectively,Vb isthecellvolume,andQ
end cooling would require a more involved expression for the is the heat generation rate. The boundary conditions are given
admittance than eq. 1, as well as an appropriate modification
to the 1D thermal model described in the following section. 1Notethat,sincetheimpedancetemperaturerelationshipisnon-linear(as
Moreover,althoughthemethodisappliedtoacylindricalcell, demonstrated in [19]), the EIS-based volume average temperature, TEIS,
is not necessarily equal to the volume average temperature, T. Although, it
the proposed approach could be applied to other geometries
shouldalsobenotedthatthetwoaretypicallycloseinvalue,particularlyif
in a similar fashion. The polynomial fit (Fig. 1) was obtained thetemperaturevariationwithinthecellissmall.
3
by: By obtaining the volume-average of eq. 2a and of its
(cid:12) partial derivative with respect to r, a two-state thermal model
∂T(r,t)(cid:12) h
∂r (cid:12)(cid:12) = −k (T(ro,t)−T∞(t)) (2b) consisting of two ODEs is obtained:
r=ro t
(cid:12) x˙ =Ax+Bu
∂T(r,t)(cid:12)(cid:12) = 0 (2c) (8)
∂r (cid:12) y=Cx+Du
r=0
whereT∞ isthetemperatureoftheheattransferfluid,andhis where x=(cid:2)T γ(cid:3)T, u=[QT∞]T and y =[Tcore Tsurf]T
the convection coefficient. A commonly employed expression arestate,inputsandoutputsrespectively.Thesystemmatrices
for the heat source in a lithium ion battery is A, B, C, and D are defined as:
Q=I(V −U )+IT∂UOCV (3) (cid:34) −48αh −15αh (cid:35)
OCV ∂T A= ro(24kt+roh) 24kt+roh
−320αh −120α(4kt+roh)
which is a simplified version of the expression first proposed ro2(24kt+roh) ro2(24kt+roh)
(cid:34) (cid:35)
α 48αh
by Bernardi et al [20]. The first term is the heat generation
B= ktVb ro(24kt+roh)
due to ohmic losses in the cell, charge transfer overpotential 0 320αh
athnidsmexapsrsestrsaionnsfearrelimmietaastiuornesd.oTnhleinceu.rTrehnetoIpaenndcviroclutaitgevoVltafgoer C=(cid:34) 2244kktt−+3rroroho2h(24k−t+1r28o0(h2r)o4kktt++1ro5hro2)h (cid:35) (9)
U is a function of SOC but is approximated here as a 24kt 15rokt
OCV 24kt+roh 48kt+2roh
constant value measured at 50 % SOC, since the HEV drive (cid:34) (cid:35)
0 4roh
cycles employed in this study operate the cell within a small D= 24kt+roh
0 roh
range of SOC (47−63%) and therefore OCV variation. If 24kt+roh
necessary, an estimator of U could also be constructed
OCV where α=k /ρc is the cell thermal diffusivity.
t p
(for example using a dynamic electrical model [21]), but for
clarity and brevity we neglect this here. The second term, the
C. Impedance Measurement
entropic heat, is neglected in this study because (i) the term
∂U /∂T is small (0 < ∂U /∂T < 0.1 mVK-1) within Eq. 1 applies to an annulus with inner radius r and outer
avg avg i
the operated range of SOC [3], and (ii) the net reversible heat radiusr .Iftheinnerradiusissufficientlysmall,thecellmay
o
would be close to zero when the cell is operating in HEV be treated as a solid cylinder, and eq. 1 becomes
mode. ˆ
Y(cid:48) = 2 r0r(cid:0)a +a T(r)+a T2(r)(cid:1)dr (10)
r2 1 2 3
B. Polynomial Approximation o 0
Substituting eq. 6 in eq. 10, the real admittance can be
A polynomial approximation (PA) is used to approximate
expressed as a function of of T , T, and γ
the solution of eq. 2a. The approximation was first introduced surf
in [6] and is described in detail in that article, although the Y(cid:48) = a +a T +3a T2+2a T2 −4a TT
1 2 3 3 surf 3 surf
essential elements are repeated here for completeness.
15a r2γ2 15a r Tγ 15a r T γ (11)
The model assumes a temperature distribution of the form + 3 o + 3 o − 3 o surf
32 8 8
(cid:18) r (cid:19)2 (cid:18) r (cid:19)4
T(r,t)=a(t)+b(t) +d(t) (4) Noting from eq. 7 that T is itself a function of T, γ and
surf
r r
o o T andthecellparameters,admittanceisultimatelyafunction
∞
The two states of the model are the volume averaged temper- of T and γ, along with the known thermal parameters and
ature T and temperature gradient γ: environmental temperature. In other words, for known values
ˆ ˆ
2 ro 2 ro (cid:18)∂T(cid:19) of ro, kt, cp, ρ, and h, the impedance is a function of the cell
T = rTdr, γ = r dr (5) state and T , thus:
r2 r2 ∂r ∞
o 0 o 0
Z(cid:48) =f(T,γ,T ) (12)
The temperature distribution is expressed as a function of T, ∞
γ, and the cell surface temperature, T :
surf
IV. FREQUENCYDOMAINANALYSIS
15r
T(r,t)=4T −3T − oγ In this section we analyze the error associated with the
surf 8
polynomial approximation, by comparing the frequency re-
(cid:20) 15r (cid:21)(cid:18) r (cid:19)2
+ −18T +18T + oγ sponse of the PA model to the frequency response of a
surf 2 ro full analytical solution of eq. (2a). We also examine the
(cid:20) 45r (cid:21)(cid:18) r (cid:19)4 approximation employed in [17], which used a combination
+ 15Tsurf −15T − 8oγ r (6) of impedance and surface temperature measurements but no
o
thermalmodel.Toachieveauniquesolutioninthatcase,itwas
Using 2b, the surface temperature can be expressed as
necessary to impose the assumption of a quadratic solution to
24k 15k r r h the temperature distribution, and so we refer to this here as
T = t T+ t o γ+ o T (7)
surf 24k +r h 48k +2r h 24k +r h ∞ the ‘quadratic assumption’ (QA) solution.
t o t o t o
4
Asin[6],thefrequencyresponsefunctionofthePAthermal mm) with LiFePO positive electrode and graphite negative
4
system, H(s), is calculated by electrode. The cell was fitted with two thermocouples, one
on the surface and another inserted into the core via a hole
H(s)=D+C(sI−A)−1B (13)
which was drilled in the positive electrode end (Fig. 4). Cell
where s = jω and I is the identity matrix. The frequency cyclingandimpedancemeasurementswerecarriedoutusinga
response of the analytical model is derived in [22], and that Biologic HCP-1005 potentiostat/booster. The impedance was
of the QA solution is derived in Appendix A. measuredusingGalvanostaticImpedanceSpectroscopywitha
200mApeak-to-peakperturbationcurrent.Theenvironmental
temperature was controlled with a Votsch VT4002 thermal
30 0 chamber. The chamber includes a fan which operates contin-
dB 0 dB −5 uously at a fixed speed during operation.
e / e / In order to calibrate the impedance against temperature,
nitud−30 nitud−10 EIS measurements were first conducted on the cell in thermal
g g
Ma−602 H11(s) = Tcore(s) / Q(s) Ma−1105 H12 = Tcore(s) / T∞(s) ethqeuisluibbrsieuqmueanttaidreanntgieficoaftitoenmopfertahteurpeosl.yTnhoemsiealexcpoeerfifimcieennttssaan1d,
Error / dB 0 Error / dB 0 a2,Daynndama3icaerxepdereismcreinbtesdwienre[1t7h]e.nconductedusingtwo3500s
−120−4 10−2 100 −1100−4 10−3 10−2 currentexcitationprofiles-thefirsttoparameterisethethermal
20 0 model, and the second to validate the identified parameters
dB 0 dB and to demonstrate the temperature estimation technique. The
e / e / −5 Analytical profiles were generated by looping over different portions of
ud−20 ud Polynomial approx.
gnit gnit−10 Qaudratic assumption anArtemisHEVdrivecycle.Theappliedcurrentswereinthe
Ma−40 Ma range −23 A to +30 A. For the duration of the experiments,
H21(s) = Tsurf(s) / Q(s) H22(s) = Tsurf(s) / T∞(s)
0.5 0.5 single frequency (215 Hz) impedance measurements were
B B
d d carriedoutevery24sandthesurfaceandinternaltemperatures
Error / 0 Error / 0 were also monitored. In order to minimise heat loss through
−0.5 −0.5 the cell ends, these were insulated using Styrofoam (Fig. 4).
10−4 10−2 100 10−4 10−3 10−2
Before each experiment, the SOC was adjusted to 50% by
Frequency (Hz) Frequency (Hz)
drawing a 0.9 C current. The temperature of the thermal
Figure 3. Comparison of frequency responses of (i) analytical solution to chamber was set to 8 ◦C and the cell was allowed to rest
the 1D cylindrical heat transfer problem, (ii) the polynomial approximation
until its temperature equilibrated before experiments began.
usedinthecurrentstudyand(iii)thequadraticassumptionusedin[17].
Besides being a function of temperature and SOC, the
Fig. 3 shows the impact of changes in heat generation on impedance is also a function of DC current, mainly due to
T andT (H andH respectively),andtheimpactof the charge transfer polarization decreasing with increasing
core surf 11 21
changes in cooling fluid temperature T on T and T current [23]. Previous results confirmed that when the EIS
∞ core surf
(H and H respectively), for each model along with the perturbationcurrentissuperposedoveranappliedDCcurrent,
12 22
error relative to the analytical solution. The results in these theimpedancemeasurementisaltered[17].Toovercomethis,
plotswereobtainedusingthethermalparametersofthe26650 the cell was allowed to rest briefly for 4 s before each EIS
cell used in the present study (Table II). The response of measurement was taken, i.e. 20 s periods of the excitation
H and H for both the PA and QA models are in good current followed by 4 s rests. The duration of this rest period
11 21
agreement with the analytical solution (note that the error was kept as short as possible to ensure that the thermal
in the response of H and H for the QA model is 0 response of the cell to the applied cycle was not significantly
21 22
since it takes the measured surface temperature as one of altered,anditwasfoundthatthecorecelltemperaturedropped
its inputs). However, the responses for both cases to changes by at most 0.25◦C during these rest periods. The issue of
in T are less satisfactory. In particular, the response H impedance measurement under DC current warrants further
∞ 12
for the QA model shows a rapid increase in error relative to investigation.
the analytical solution at frequencies above ∼ 10−3 Hz. The
PA performs well up to higher frequencies, although its error
VI. MODELPARAMETERISATION&VALIDATION
in H and H become unsatisfactory above ∼ 10−2 Hz.
21 22
However, the frequency range at which the PA agrees with Parameterisation is performed to estimate the values of
the analytical solution is broader than that of the QA model, kt, cp, and h. The density ρ was identified in advance by
andisconsideredsatisfactorygiventheslowrateatwhichthe measuring the cell mass and dividing by its volume. The
coolingfluidtemperaturefluctuatesinatypicalbatterythermal valuesfromthefirstexcitationprofile(comprisingcellcurrent,
management system. voltage plus surface, core and chamber temperatures) were
used for the estimation.
V. EXPERIMENTAL The parameterisation was carried out offline using fmin-
Experiments were carried out with a 2.3 Ah cylindrical cell search in Matlab to minimise the magnitude of the Euclidean
(A123 Model ANR26650 m1-A, length 65 mm, diameter 26 distancebetweenthemeasuredandestimatedcoreandsurface
5
a
Temperaturedsensors all of these studies), heat generation in the contact resistances
(fordvalidation): T ,T between the cell and connecting wires and/or measurement
core surf
uncertainty in the temperature. The convection coefficient is
Thermaldchamber
within the range expected of forced convection air cooling
[24].
Parameter Units Reference Initial Identified
ρ kgm-3 2047-2118 - 2107
[4],[25],[26]
cp Jkg-1 K-1 1004.9-1109.2 1050 1171.6
[3],[7],[4]
e- kt Wm-1 K-1 [202.4],88[2-05.]6,9[04] 0.55 0.404
h Wm-2 - 20 39.3
TableII
40 30 COMPARISONOFREFERENCE&ESTIMATEDPARAMETERS
20
Id(A)-200 T20
-40 10
The measured core and surface temperatures (subscript
0 500 1000 1500 2000 0 500 1000 1500 2000
td(s) td(s) ‘exp’) and the corresponding model predictions (subscript
Powerdsupplydvdimpedance Datadacquisitiondmodule ‘m’) for the parameterised model are shown in Fig. 5. The
root-mean-square errors (RMSE) in the surface and core
b c
temperatures are 0.19 ◦C and 0.18 ◦C respectively.
30
25
20
C
°T / 15
T T
core,exp core,m
10
T T
surf,exp surf,m
5
0 500 1000 1500 2000 2500 3000 3500
time / s
Figure 5. Model parameterisation: Comparison between measured and
d
predictedcoreandsurfacetemperaturesintheparameterisedmodel.
The model with identified parameters was validated against
the second current excitation profile (Fig. 6). The RMSEs
in the core and surface temperatures were 0.21 ◦C and 0.16
◦C respectively in this case. These errors are only marginally
greater than those in the parameterisation test, indicating that
the estimation is satisfactory.
30
25
20
C
°T / 15
T T
Figure4. Experimentalsetupforparameteridentificationandvalidation.(a) core,exp core,m
schematic diagram with cutaway view showing cell core and jelly roll, (b) 10 T T
surf,exp surf,m
cell drilling procedure, (c) prepared cell inside thermal chamber, (d) power 5
supply&thermalchamber. 0 500 1000 1500 2000 2500 3000 3500
time / s
Figure 6. Model validation. Comparison between measured and predicted
temperatures, as in [4]. Table II compares the identified core and surface temperatures in the parameterised model applied to the
secondcurrentexcitationprofile.
parameters with the initial guesses and parameters for the
same cell from the literature. The estimated parameters are
close to those reported in the literature. The deviations may VII. STATEESTIMATION
beattributedtomanufacturingvariability,errorintheheatgen- Kimetal.[6]showedthattheeffectofchangestothevalue
eration calculation (due to the omission of entropic heating in of h on the predicted surface and core temperatures is greater
6
than the effect of changes to the other thermal parameters. where Kx is the Kalman gain for the state, and Hx is the
k k
Moreover, h depends strongly on the thermal management Jacobian matrix of partial derivatives of f with respect to x:
system settings and its calculation often relies on empirical (cid:12)
∂f(x ,h )(cid:12)
correlations between coolant flow rates and heat transfer. Hx = k k (cid:12) (26)
k ∂x (cid:12)
Tonhluins,etdhuerreingisoapenreaetidont.oTidheisntsifeyctitohne coountlvinecetsiotnhecouesfeficoifenat k (cid:12)xk=xˆ−k
Themeasurementupdateprocessesfortheparameterfilterare:
dual extended Kalman filter (DEKF) [27] for estimating the
core and surface temperatures and the convection coefficient. Kh =(Ph)−(Hh)T (cid:0)Hh(Ph)−(Hh)T +Rn(cid:1)−1 (27)
k k k k k k
The DEKF reduces to an EKF if the convection coefficient is (cid:16) (cid:17)
hˆ =hˆ−+Kh z −f(xˆ , hˆ−) (28)
assumed known and provided to the model in advance. k k k k k k
We firstly modify eq. 9 by rewriting it as a discrete Ph =(I−KhHh)(Ph)− (29)
k k k k
time model, setting the impedance as the model output, and
explicitly including the dependency on the parameter, hk: where Hkh is the Jacobian matrix of partial derivatives of f
with respect to h:
xk+1 =A¯(hk)xk+B¯(hk)uk+vk (14) (cid:12)
∂f(x ,h )(cid:12)
yk =f(xk, hk)+nk (15) Hkh = ∂hk k (cid:12)(cid:12) (30)
hk+1 =hk+ek (16) k (cid:12)hk=hˆ−k
where y = Z(cid:48) and f(x , h ) is the non-linear function TheabovealgorithmcanbesimplifiedtoastandardEKFby
k k k
omitting the parameter update processes (eqs. 16, 18-19 and
relating the state vector to the measurement (i.e. eq. 12),
27-30)andassumingtheconvectioncoefficientisfixed.Inthe
and v , n and e are the noise inputs of the covariance
k k k
matrices Rv, Rn and Re. The states, inputs and measured following section we investigate the performance of both the
outputs are thus x = (cid:2)T γ(cid:3)T, u = [QT ]T and y = Z(cid:48). baseline EKF and the full DEKF algorithm.
∞
Note that, although the impedance is now the model output,
VIII. RESULTS
the core and surface temperatures are also computed from the
identified states and parameter at each time step using eq. 8, We first investigate the performance of an EKF estimator
for validation against the thermocouple measurements. A¯ and whereby the convection coefficient is provided and assumed
B¯ are system matrices in the discrete-time domain, given by fixed. We then compare the performance of the DEKF algo-
rithm with that of the baseline EKF when an incorrect initial
A¯ =e(A∆t), B¯ =A−1(A¯ −I)B (17)
estimate of the convection coefficient is provided. Lastly, we
compare the performance of the DEKF with that of a dual
where ∆t is the sampling time of 1 s. The update processes
Kalman filter (DKF) based on the same thermal model but
are then given as follows. The time update processes for the
with T as the measurement input rather than Z(cid:48).
parameter filter are: surf
hˆ− =hˆ (18)
k k−1 A. Convection Coefficient Known
(Ph)− =Ph +Re (19)
k k−1 The initial state estimate provided to the battery is xˆ0 =
wherehˆ− andhˆ aretheaprioriandaposterioriestimatesof [250], i.e. the battery has a uniform temperature distribution
k k at25◦C.Thetrueinitialbatterystateisauniformtemperature
the parameter h, and (Ph)− and Ph are the corresponding
k k−1 distributionat8◦C.Thecovariancematricesarecalculatedas
error covariances.
Rn =σ2 and Rv =β2diag(2,2). The first tuning parameter
The time update processes for the state filter are: n v
is chosen as σ = 1×10−4 Ω, which is a rough estimate
n
xˆ−k =A¯k−1xˆk−1+B¯k−1uk−1 (20) of the standard deviation of the impedance measurement. The
(Pxk)− =A¯k−1Pxk−1A¯Tk−1+Rv (21) second tuning parameter was chosen as βv =0.1, by trial and
error. Fig. 7 shows that, using the EKF, the core and surface
where xˆ−k and xˆk are the a priori and a posteriori estimates temperatures quickly converge to the correct values and are
ofthestate,and(Pxk)− andPxk−1 arethecorrespondingerror accurately estimated throughout the entire excitation profile.
covariances. The matrices A¯k−1 and B¯k−1 are calculated by: The RMSEs of core and surface temperature are 1.35 ◦C and
A¯k−1 =A¯(h)(cid:12)(cid:12)h=hˆ−k, B¯k−1 =B¯(h)(cid:12)(cid:12)h=hˆ−k (22) l1o.3o4p ◦mCodreelsp(escutbivscerlyip.tIn‘mc’o)nwtriatsht,ntohemReMasuSrEesmefonrt tfheeedobpacekn
Since the relationship between impedance and the cell state are 6.66 ◦C and 4.42 ◦C respectively. It should be noted that
is non-linear, the measurement model must be linearised since the uncertainty of the impedance measurement typically
about the predicted observation at each measurement. The increases as impedance decreases, the temperature estimates
measurement update equations for the state filter are: could be more uncertain at higher temperatures. Hence, the
implementation of this technique could be more challenging
Kx =(Px)−(Hx)T (cid:0)Hx(Px)−(Hx)T +Rn(cid:1)−1 (23)
k k k k k k at higher ambient temperature conditions than those studied
(cid:16) (cid:17)
xˆ =xˆ−+Kx z −f(xˆ−, hˆ−) (24) here.
k k k k k k
Itshouldbenotedthatwealsoachievedsimilarperformance
Pxk =(I−KxkHxk)(Pxk)− (25) using a simpler EKF based on Z(cid:48) with the assumption that
7
45
17 27
40 16 C 4 (Tcore,EKF − Tcore,exp) (Tcore,DEKF − Tcore,exp)
35 15 25 °or / 2
err
30 0

C 0 500 1000 1500 2000 2500 3000 3500
°T / 25 45
20 40 22
18
15 T TT T 35 16 20
10 core,exp ccoorree,,mm core,EKF 30
T TT T
surf,exp ssuurrff,,mm surf,EKF C 25
50 500 1000 1500 2000 2500 3000 3500 °T / 20
Time / s
15
Figure7. TemperatureresultsforEKFusingZ(cid:48) asmeasurementinput 10 Tcore,exp Tcore,EKF Tcore,DEKF
T T T
5 surf,exp surf,EKF surf,DEKF

0 500 1000 1500 2000 2500 3000 3500
the impedance is related directly to T rather than to TEIS. 100
However,thisassumptionmayleadtounsatisfactoryresultsfor −1k 80
cellswithalargerradiusorwhenlargertemperaturegradients −2Wm 60 htrue hEKF hDEKF
exist within the cell. Moreover, since this approach assumes h / 40
that the impedance is a function of the state only (and not 20
0 500 1000 1500 2000 2500 3000 3500
the parameter h), it is not suitable for the application of the t / s
DEKF discussed in the following section.
Figure8. TemperatureresultsforDEKFusingZ(cid:48) asmeasurementinput
B. Convection Coefficient Unknown
EKF in this case overpredicts the core temperature by a
Next we investigate the performance of the DEKF. The
greatermarginthantheEKFbasedonZ(cid:48),althoughthesurface
same incorrect initial state estimate is provided to the battery,
temperature estimate is much more accurate. This is because
xˆ = [250]. Moreover, an incorrect initial estimate for the
0 the surface thermocouple ensures an accurate surface temper-
convection coefficient is provided, hˆ = 2 × h . This
0 true ature estimate and to reconcile this with the overestimated
value of h is also provided to the EKF. The error covariance
convection coefficient, the core temperature estimate is forced
matrix for the parameter estimator is Re = β2 where the
e to be much greater. The DKF correctly identifies the correct
tuning parameter is chosen as β = 2.5. Fig. 8 compares
e convection coefficient, in the same way as the DEKF. Since
the results of both of these cases against the thermocouple
the thermocouple measurement exhibits less noise than the
measurements. The EKF is shown to overestimate the core
impedance measurement, the model converges to the correct
temperature and underestimate the surface temperature for
estimate for h more quickly than in the DEFK case, as shown
the duration of the experiments. This is expected, since the
by the RMSE values in Table III.
impedance measurement ensures the accuracy of the volume
averaged cell temperature but the model assumes that the
Method 0<t<3500s 1200<t<3500s
convection coefficient is higher than in reality, and therefore Tcore Tsurf Tcore Tsurf
the temperature difference across the cell is overestimated. In EKF+Z(cid:48) 2.04 2.06 1.79 1.98
contrast, the DEKF corrects the convection coefficient, and KF+Tsurf 2.49 1.44 2.90 1.58
thusimprovestheaccuracyofthesubsequenttemperaturepre- DEKF+Z(cid:48) 1.43 1.24 0.47 0.42
dictions.Thisisevidentfromtheerrorsinthecoretemperature DKF+Tsurf 0.36 0.33 0.16 0.14
estimate(topplotofFig.8),whichinitiallyaresimilarinboth TableIII
COMPARISONOFRMSESFORCOREANDSURFACETEMPERATURES(◦C)
cases but drop to much smaller values for the DEKF once the
WITHUNKNOWNCONVECTIONCOEFFICIENT.
correctconvectioncoefficientisidentified.TheRMSEsofcore
and surface temperature in each case are shown in Table III,
along with the values for the time period, 1200 < t < 3500 s In conclusion, the temperature and convection coefficient
(i.e. after the convection coefficient value has converged). estimators using Z(cid:48) as measurement input are capable of ac-
Finally, we investigate the case where T is used as the curately estimating the core and surface temperatures and the
surf
measurement (y = T ) to the estimator rather than Z(cid:48). convection coefficient. The performance is comparable to that
surf
This results in a linear KF and DKF exactly equivalent to of an estimator using the same thermal model coupled with
that studied in [4]. The tuning parameters for the covariance surfacetemperaturemeasurements.Moreover,theperformance
matrices are also chosen to be the same as those employed of the present method is superior to that of previous methods
in [4]. The same initial state and parameter estimates are based on impedance measurements alone, which only provide
provided as for the DEKF. Fig. 9 shows that the standard an estimate of the average internal temperature of the cell.
8
ACKNOWLEDGEMENTS
8 This work was funded by a NUI Travelling Scholarship,
C 6 (Tcore,EKF − Tcore,exp) (Tcore,DEKF − Tcore,exp) a UK EPSRC Doctoral Training Award and the Foley-Bejar
°or / 24 scholarship from Balliol College, University of Oxford. This
err 0 publicationalsobenefitedfromequipmentfundedbytheJohn
−2
0 500 1000 1500 2000 2500 3000 3500 Fell Oxford University Press (OUP) Research Fund. Finally,
45 the authors would like to thank Peter Ireland and Adrien

40 Bizeray for valuable comments, and Robin Vincent for the
18 22
35 facilities to prepare the instrumented battery.
16 20
30
°T / C2205 A. Frequency Domain AnAalPyPsEisNoDfIXQuadratic Assumption
15
The analysis leading to the frequency response plots of the
10 Tcore,exp Tcore,EKF Tcore,DEKF
QA model in Fig. 3 is outlined in this section. The QA model
5 Tsurf,exp Tsurf,EKF Tsurf,DEKF
wasusedin[17]toobtainauniquesolutionforthetemperature
0 500 1000 1500 2000 2500 3000 3500
distributionwhentheimpedanceandsurfacetemperaturewere
100
−1k 80 measured but the cell thermal properties and heat generation
−2m 60 h h h rates were assumed unknown. To achieve a unique solution
W true EKF DEKF in that case, it was necessary to impose the assumption of a
h / 2400 quadratictemperaturedistributionbasedonthesolutionofthe
0 500 1000 1500 2000 2500 3000 3500
1D heat equation at steady state.
t / s
Muratori et al. [22] showed that the time domain PDE of
Figure9. DKFusingTsurf asmeasurementinput. eq.2a,canbetransformedintoanequivalentODEproblemin
thefrequencydomainwhichcanbesolvedanalytically.Using
thisapproachthesolutionforthetemperaturesatthecoreand
surface of the cell are shown to be:
IX. CONCLUSIONS
T (s)= 1 Q(s)+ kht(T∞(s)− k1a2Q(s)) (31)
core k a2 hJ (ar )−aJ (ar )
Impedancetemperaturedetectionenablesbothcoreandsur- t kt 0 o 1 o
facetemperatureestimationwithoutusingtemperaturesensors. T (s)= 1 Q(s)+ kht(T∞(s)− kt1a2Q(s))J (ar )
In this study, the use of ITD as measurement input to a cell surf k a2 hJ (ar )−aJ (ar ) 0 o
thermal model is demonstrated for the first time. t kt 0 o 1 o (32)
Previously, we estimated cell temperature distribution by
where a2 = sα−1, Q(s) is the transform of Q(t), and J
i
combining ITD with a surface temperature measurement and
is the ith-order Bessel function of the first kind [28]. Eqs.
imposing a quadratic assumption on the radial temperature
31 and 32 can be interpreted as the outputs of a continuous
profile.FrequencydomainanalysisshowsthattheQAsolution
time dynamic system [22], where u(t)=[Q(t), T (t)]T and
∞
may be inaccurate if the temperature of the cooling fluid has
y(t) = [T (t), T (t)]T, such that the solution of the
fluctuationsontheorderof10−3Hzorhigher.ThePAthermal core surf
BVP is equivalent to the impulse response of the system:
model used in the present study is robust to higher frequency
(cid:20) (cid:21) (cid:20) (cid:21)(cid:20) (cid:21)
fluctuations (∼10−2 Hz). Tcore(s) = H11(s) H12(s) Q(s) (33)
T (s) H (s) H (s) T (s)
AnEKFusingaparameterisedPAthermalmodelwithITD surf 21 22 ∞
measurement input is shown to accurately predict core and where the H matrix is formed by the transfer functions:
surface temperatures for a current excitation profile based on
an Artemis HEV drive cycle. A DEKF based on the same H (s)= 1 khtJ0(aro)−aJ1(aro)− kht (34)
thermalmodelandmeasurementinputiscapableofaccurately 11 kta2 khtJ0(aro)−aJ1(aro)
identifying the convection coefficient when the latter is not h
provided to the model in advance. The performance of the H (s)= kt (35)
12 hJ (ar )−aJ (ar )
DEKF using impedance as measurement input is comparable kt 0 o 1 o
to an equivalent DKF estimator using surface temperature as H (s)= 1 −aJ1(aro) (36)
measurementinput,althoughthelatterisslightlysuperiordue 21 k a2 hJ (ar )−aJ (ar )
t kt 0 o 1 o
to the higher accuracy of the thermocouple.
hJ (ar )
Future work will investigate the application of ITD to H22(s)= hJ (akrt )0−aoJ (ar ) (37)
multiple cells in a battery pack, as well as investigating kt 0 o 1 o
methodsofcombiningimpedancewithconventionalsensorsto This system of transfer functions gives frequency responses
enablemorerobusttemperaturemonitoringandfaultdetection, for the analytical solution results plotted in Fig. 3. Using a
and self-calibration. similar approach, we can develop an analytical solution to the
9
QA model used in [17], where we denote the new transfer where the H matrix is formed by the functions:
function HQA(s). In this case, both the volume-averaged cell k a2r J (ar )−2hr aJ (ar )+4hJ (ar )
temperature (identified via the impedance2) and the surface HQA,11(s)=− t ok1r a3o(hJ (aro)−0k aoJ (ar ))1 o
temperature were measured directly, and it is assumed that t o 0 o t 1 o
(44)
no other information was available. Since these two inputs
alone are not sufficient to achieve a unique solution for the H (s)= −hkta3J0aro+4hkta2J1aro (45)
temperature distribution, it was also necessary to impose the QA,21 k r a3(hJ ar −k aJ (ar ))
t o 0 o t 1 o
following constraint on the temperature profile: H (s)=H (s) (46)
QA,12 21
(cid:18) r2(cid:19) HQA,22(s)=H22(s) (47)
T (r)=T +(T −T ) 1− (38)
QA surf core,QA surf ro2 REFERENCES
[1] U. S. Kim, C. B. Shin, and C.-S. Kim, “Modeling for the scale-up of
This is the 1D steady-state solution of the heat equation for
a lithium-ion polymer battery,” J. Power Sources, vol. 189, no. 1, pp.
a cylinder with uniform heat generation [24], with T(r) and 841–846,Apr.2009.
T replaced by T (r) and T . [2] Q. Wang, P. Ping, X. Zhao, G. Chu, J. Sun, and C. Chen, “Thermal
core QA core,QA runaway caused fire and explosion of lithium ion battery,” J. Power
The volume average of the QA temperature distribution is Sources,vol.208,pp.210–224,Jun.2012.
set equal to the true volume average temperature, giving: [3] C.Forgez, D.Vinh Do,G.Friedrich, M.Morcrette,and C.Delacourt,
“Thermal modeling of a cylindrical LiFePO4/graphite lithium-ion bat-
ˆ
2 ro tery,”J.PowerSources,vol.195,no.9,pp.2961–2968,May2010.
T = rT (r)dr =T (39) [4] Y.Kim,S.Mohan,S.Member,J.B.Siegel,A.G.Stefanopoulou,and
QA r2 QA Y. Ding, “The Estimation of Temperature Distribution in Cylindrical
o 0
BatteryCellsUnderUnknownCoolingConditions,”IEEETrans.Con-
trolSyst.Technol.,pp.1–10,2014.
Substituting eq. 38 in eq. 39 and integrating, the QA core
[5] X.Lin,H.E.Perez,J.B.Siegel,A.G.Stefanopoulou,R.D.Anderson,
temperature becomes and M. P. Castanier, “Online Parameterization of Lumped Thermal
Dynamics in Cylindrical Lithium Ion Batteries for Core Temperature
EstimationandHealthMonitoring,”IEEETrans.ControlSyst.Technol.,
T =2T −T (40)
core,QA surf vol.21,no.5,pp.1745–1755,Sep.2013.
[6] Y. Kim, J. B. Siegel, and A. G. Stefanopoulou, “A computationally
Thus, an estimate for the core temperature is obtained di- efficient thermal model of cylindrical battery cells for the estimation
of radially distributed temperatures,” in American Control Conference
rectly from the surface and volume averaged temperature
(ACC),Washington,DC,2013,pp.698–703.
measurements. Since the surface temperature is measured, [7] X. Lin, H. E. Perez, S. Mohan, J. B. Siegel, A. G. Stefanopoulou,
the values of H (s) and H (s) are identical to Y. Ding, and M. P. Castanier, “A lumped-parameter electro-thermal
QA,21 QA,22
model for cylindrical batteries,” J. Power Sources, vol. 257, pp. 1–11,
the corresponding values of the unsteady thermal model.
Jul.2014.
The values of HQA,11(s) and HQA,12(s) can be obtained [8] A.Pesaran,G.H.Kim,andM.Keyser,“IntegrationIssuesofCellsinto
as follows: Substituting T from eq. 5 into eq. 40, the QA Battery Packs for Plug-In and Hybrid Electric Vehicles,” in EVS-Â24
InternationalBattery,HybridandFuelCellElectricVehicleSymposium,
approximation of the core temperature can be expressed as a
no.May,Stavanger,Norway,2009.
function of the temperature distribution: [9] R. Srinivasan, “Monitoring dynamic thermal behavior of the carbon
ˆ anode in a lithium-ion cell using a four-probe technique,” J. Power
4 ro Sources,vol.198,pp.351–358,Jan.2012.
T = (rT(r)dr−T )dr (41) [10] J.P.Schmidt,S.Arnold,A.Loges,D.Werner,T.Wetzel,andE.Ivers-
core,QA ro2 0 surf Tiffée, “Measurement of the internal cell temperature via impedance:
evaluation and application of a new method,” J. Power Sources, vol.
Substituting eqs. 31 and 32, f´or T(r) and Tsurf, respectively, 243,pp.110–117,Jun.2013.
and integrating (noting that rorJ (ar)dr = r J (ar )/a), [11] L. Raijmakers, D. Danilov, J. van Lammeren, M. Lammers, and
0 0 o 1 o P. Notten, “Sensorless battery temperature measurements based on
we obtain: electrochemical impedance spectroscopy,” J. Power Sources, vol. 247,
pp.539–544,Feb.2014.
[12] R. Srinivasan, B. G. Carkhuff, M. H. Butler, and A. C. Baisden,
Q(s)
T = + “Instantaneous measurement of the internal temperature in lithium-ion
core,QA k a2 rechargeablecells,”ElectrochimicaActa,vol.56,no.17,pp.6198–6204,
t
kht(T∞(s)− kt1a2Q(s))(cid:20)4J1(aro) −J (ar )(cid:21) (42) [13] JJ.ul.Z2h0u1,1Z.. Sun, X. Wei, and H. Dai, “A new lithium-ion battery
hJ (ar )−aJ (ar ) ar 0 o internaltemperatureon-lineestimatemethodbasedonelectrochemical
kt 0 o 1 o o impedancespectroscopymeasurement,”JournalofPowerSources,vol.
274,pp.990–1004,Jan.2015.
Thus, the QA system model is given by: [14] D. A. Howey, P. D. Mitcheson, S. Member, V. Yufit, G. J. Offer, and
N.P.Brandon,“On-linemeasurementofbatteryimpedanceusingmotor
(cid:20) (cid:21) (cid:20) (cid:21)(cid:20) (cid:21) controller excitation,” IEEE Trans. Veh. Technol., vol. 63, no. 6, pp.
T (s) H (s) H (s) Q(s)
core,QA = QA,11 QA,12 2557–2566,2014.
Tsurf,QA(s) HQA,21(s) HQA,22(s) T∞(s) [15] N. Brandon, P. Mitcheson, D. A. Howey, V. Yufit, and G. J. Offer,
(43) “Battery Monitoring in Electric Vehicles, Hybrid Vehicles and other
Applications,”2012.
[16] W.HuangandJ.AbuQahouq,“AnOnlineBatteryImpedanceMeasure-
2AsdiscussedinSectionI,theimpedanceactuallyidentifiesTEIS,which mentMethodUsingDC-DCPowerConverterControl,”IEEETrans.Ind.
is not necessarily equal to T. However, the assumption that TEIS = T Electron.,vol.0046,no.c,pp.1–1,2014.
issatisfactoryforthepurposeofidentifyingthefrequenciesatwhicherrors [17] R. R. Richardson, P. T. Ireland, and D. a. Howey, “Battery internal
relativetotheanalyticalsolutionbecomesignificant,andisonlyusedforthis temperatureestimationbycombinedimpedanceandsurfacetemperature
purposeinthisarticle. measurement,”J.PowerSources,vol.265,pp.254–261,May2014.
10
[18] M. Fleckenstein, O. Bohlen, M. a. Roscher, and B. Bäker, “Current
densityandstateofchargeinhomogeneitiesinLi-Ionbatterycellswith
LiFePO4 as cathode material due to temperature gradients,” J. Power
Sources,vol.196,no.10,pp.4769–4778,Jan.2011.
[19] Y. Troxler, B. Wu, M. Marinescu, V. Yufit, Y. Patel, A. J. Marquis,
N.P.Brandon,andG.J.Offer,“Theeffectofthermalgradientsonthe
performance of lithium-ion batteries,” Journal of Power Sources, no.
247,pp.1018–1025,Jun.2013.
[20] D. Bernardi, E. Pawlikowski, and J. Newman, “A General Energy
Balance for Battery Systems,” J. Electrochem. Soc., vol. 132, no. 1,
pp.5–12,1985.
[21] C.R.BirklandD.A.Howey,“Modelidentificationandparameteresti-
mationforLiFePO4,”inIETHybridandElectricVehiclesConference
2013(HEVC2013),2013,pp.1–6.
[22] M.Muratori,N.Ma,M.Canova,andY.Guezennec,“AModelOrderRe-
ductionMethodfortheTemperatureEstimationinaCylindricalLi-Ion
Battery Cell,” ASME 2010 Dynamic Systems and Control Conference,
Volume1,no.7,pp.633–640,2010.
[23] B. Ratnakumar, M. Smart, L. Whitcanack, and R. Ewell, “The
impedancecharacteristicsofMarsExplorationRoverLi-ionbatteries,”
J.PowerSources,vol.159,no.2,pp.1428–1439,Sep.2006.
[24] F. P. Incropera and D. P. De Witt, Fundamentals of Heat and Mass
Transfer,6thed. Wiley,2007.
[25] H. Khasawneh, J. Neal, M. Canova, and Y. Guezennec, “Analysis of
heat-spreadingthermalmanagementsolutionsforlithium-ionbatteries,”
in ASME 2011 International Mechanical Engineering Congress and
Exposition,2011,pp.421–428.
[26] X.Li,F.He,andL.Ma,“Thermalmanagementofcylindricalbatteries
investigatedusingwindtunneltestingandcomputationalfluiddynamics
simulation,”J.PowerSources,vol.238,pp.395–402,Sep.2013.
[27] E. A. Wan and A. T. Nelson, “Dual extended Kalman filter methods.”
inKalmanfilteringandneuralnetworks. NewYork,NY:Wiley,2001,
pp.123–173.
[28] E.Kreyszig,Advancedengineeringmathematics. JohnWiley&Sons,
2010.
Robert R. Richardson receivedtheB.Eng.degree
in mechanical engineering from the National Uni-
versityofIreland,Galway,in2012.Heiscurrently
workingtowardtheD.Phil(PhD)degreeintheEn-
ergyandPowerGroup,DepartmentofEngineering
Science,UniversityofOxford,Oxford,U.K,where
hisresearchisfocusedonimpedancebasedbattery
thermalmanagement.
David A. Howey (M’10) received the B.A. and
M.Eng. degrees from Cambridge University, Cam-
bridge, U.K., in 2002 and the Ph.D. degree from
Imperial College London, London, U.K., in 2010.
HeiscurrentlyaLecturerwiththeEnergyandPower
Group,DepartmentofEngineeringScience,Univer-
sity of Oxford, Oxford, U.K. He leads projects on
fast electrochemical modeling, model-based battery
managementsystems,batterythermalmanagement,
andmotordegradation.Hisresearchinterestsinclude
condition monitoring and management of electric-
vehicleelectric-vehiclecomponents.