- Email: [email protected]

Contents lists available at ScienceDirect

Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour

Spatiotemporal modeling of internal states distribution for lithiumion battery Mingliang Wang a, b, Han-Xiong Li a, * a

Department of Systems Engineering and Engineering Management, City University of Hong Kong, Tat Chee Avenue, Kowloon, Hong Kong Special Administrative Region b State Key Laboratory of High Performance Complex Manufacturing, School of Mechanical and Electrical Engineering, Central South University, Changsha, China

h i g h l i g h t s Data based spatiotemporal models are developed for lithium-ion batteries. Low computation cost and high accuracy of the model enables the online applications. The internal states distribution can be estimated in real-time with the developed model. New potential indicators are proposed to prevent over-voltage and under-voltage.

a r t i c l e i n f o

a b s t r a c t

Article history: Received 31 January 2015 Received in revised form 8 September 2015 Accepted 28 September 2015

Electrochemical properties of the battery are described in partial differential equations that are impossible to compute online. These internal states are spatially distributed and thus difﬁcult to measure in the battery operation. A space-time separation method is applied to model the electrochemical properties of the battery with the help of the extended Kalman ﬁlter. The model is efﬁciently optimized by using LASSO adaptation method and can be updated through data-based learning. The analytical model derived is able to offer a fast estimation of internal states of the battery, and thus has potential to become a prediction model for battery management system. © 2015 Elsevier B.V. All rights reserved.

Keywords: Lithium-ion batteries Lower order modeling Spatiotemporal modeling KarhuneneLoeve decomposition Extended Kalman ﬁlter

1. Introduction Hybrid electrical vehicles (HEV) are prevailing worldwide nowadays due to their environmental friendliness. Lithium-ion batteries are especially suitable for powering the HEV given their advantages in high energy density, weak memory effect and low environmental pollution. Battery management system aims to carry out the real-time monitoring of the battery status to ensure safety and health of batteries. Since many of battery states are not measurable, modeling of these immeasurable states become critical to the online monitoring. For physical modeling of electrochemical property of lithiumion battery, the one dimensional model is ﬁrstly proposed by

* Corresponding author. E-mail address: [email protected] (H.-X. Li). http://dx.doi.org/10.1016/j.jpowsour.2015.09.107 0378-7753/© 2015 Elsevier B.V. All rights reserved.

Doyle, Fuller and Newman [1], and further developed by Ramadass in the name of “Pseudo 2D model” by adding the capacity fading mechanism [2,3]. Though their work can simulate the small-scale battery cell well, it cannot work for the large-scale battery because the potential drop will be big to inﬂuence the current distribution in larger space. For a better modeling of battery physics, a 2D scale-up model is ﬁrst developed for the large-scale battery [4], and then, the multi-scale and multi-dimension model (MSMD) is further proposed for better description of electrochemical property [5]. In summary, both the electrochemical property and the thermal property of the battery belong to distributed parameter systems (DPS) that are described in partial differential equations (PDE), and have temporal dynamics distributed over space. The above physical models are unsuitable for real-time operation due to heavy computation. For fast computation online,

262

M. Wang, H.-X. Li / Journal of Power Sources 301 (2016) 261e270

2. Problem description A two dimensional model of Lithium ion battery is popularly used in modeling the large-scale battery [4,17]. As illustrated in Fig. 1, the current ﬂux between the electrodes would be perpendicular to the electrodes and homogeneous along its ﬂowing direction under the assumption that the distance between the two current collectors is very close. The domains Up and Un denote positive electrode and negative electrode, respectively. In real operations, the known inputs would be the boundary conditions: the applied current at the tab (iapp) measured by a Galvanometer and the terminal voltage (vterminal) by a Voltmeter. The mathematical equations for this model are presented in Appendix A. More Detailed information can be found in Refs. [4,17,18]. Difﬁculties of battery modeling lie in the following aspects:

Fig. 1. Schematic diagram of current ﬂow in battery cell.

analytical models simpliﬁed in ordinary differential equations (ODE) are needed. Principle component analysis (PCA) or ve (K-L) decomposition method is popularly used to KarhuneneLoe develop such models for DoyleeFullereNewman model [6e9]. However, no work has been reported for two dimensional scale-up model and MSMD model. Since the internal states of the battery are not measurable, the state observer is needed for the analytical modeling. The extended Kalman ﬁlter (EKF) [10e12], Sigma-point Kalman Filter [13,14], and Unscented Kalman Filter [15,16] are presented in modeling the depth-of-discharge (DOD) or state-of-charge. However, all these work are applied to the equivalent circuit model that is described in ODE. Since the battery is a typical spatiotemporal system, and up to now, no state observer has been reported for such kind of system. In this paper, we developed a spatiotemporal modeling method for online estimation of internal states of batteries, including DOD distribution, potential distributions, etc. First, the spatiotemporal dynamics of the battery will be decoupled into spatial basis functions and temporal variables. Space basis functions can be learned using KL method. Since the internal states DOD are immeasurable, the extended Kalman ﬁlter will be used and implemented in the forward-backward structure. The parameters can be identiﬁed more efﬁciently by using the LASSO adaptation method. Finally, after synthesizing the space functions and the estimated temporal variables, the spatiotemporal dynamics of the battery can be effectively constructed.

1) Nonlinear space/time coupled dynamics: the governing equations of the battery are nonlinear partial differential equations as given in Appendix A. This spatiotemporal dynamics is difﬁcult to solve analytically and consumes extremely high computation for discrete solutions. 2) Uncertainty: arises due to the environmental disturbance, measurements inaccuracy, and battery aging. All these timevarying disturbances are difﬁcult to model. 3) Immeasurable states: internal states of the battery system, such as, depth-of-discharge and potentials, etc. cannot be measured online. This makes the modeling more difﬁcult. In rest of this paper, we will propose a spatiotemporal model based state prediction for the lithium-ion battery. First, the spacetime separation scheme will be applied to decouple spatiotemporal dynamics [19]. Then the extended Kalman ﬁlter will be designed to estimate the internal states of the battery in a forwardbackward structure [20,21]. The parameters in the ﬁlter will be efﬁciently optimized by using the LASSO adaptation method with uncertainty compensated from data-based learning. 3. Space-time separation based online modeling There are two critical internal states, DOD, and potential of positive/negative electrodes, which need to model. Since these two variables are spatially distributed, the spatiotemporal models are needed for a realistic estimation. The measurable boundary condition, (vterminal) and (iapp) will be the critical input to the modeling. Since the DOD has no direction connection to these boundary conditions, the extended Kalman ﬁlter is needed. 3.1. Modeling of depth-of-discharge (DOD) Since the DOD is an immeasurable state variable in the battery,

Fig. 2. Space-time separation based online modeling.

M. Wang, H.-X. Li / Journal of Power Sources 301 (2016) 261e270

the modeling will consist of space/time separation method and the extended Kalman ﬁlter as shown in Fig. 2. The modeling process includes three major steps. A space-time separation is ﬁrst applied to have the space variables fi(z) identiﬁed by the KL decomposition method (see Appendix B). The time variables qi(t) will then be estimated by the extended Kalman Filter [20,21] and optimized with LASSO adaptation method (see Appendix D), all these will be implemented in a forward-backward structure.

Z

Zt Z Dðz; tÞ,fj ðzÞdz ¼

U

263

Since D(z,t) is a function of J(z,t) as shown in Eq. (A.6), substituting Eq. (3) into Eq. (A.6), we will have

Dðz; tÞ ¼

Z

h a4 þ a5 Dðz; tÞ þ a6 Dðz; tÞ2 hðz; tÞ a0 þ a1 Dðz; tÞ i þ a2 Dðz; tÞ2 þ a3 Dðz; tÞ3 dt=QT (4)

By multiplying fj(z) to each side of the Eq. (4) and integrating them spatially, we will have,

h i fðzÞ a4 þ a5 Dðz; tÞ þ a6 Dðz; tÞ2 hðz; tÞ a0 þ a1 Dðz; tÞ þ a2 Dðz; tÞ2 þ a3 Dðz; tÞ3 dtdz=QT

(5)

0 U

Synthesis of the identiﬁed space variable and the estimated time variable will provide a good prediction of the battery states.

Z fj ðzÞ U

X

!

Zt Z

fi ðzÞqi ðtÞ dz ¼

i

2 fj ðzÞ4a4 þ a5

! fi ðzÞqi ðtÞ

i

0 U

þ a2

X

X

fi ðzÞqi ðtÞ

X

0 !2 3 2 X X 5 4 fi ðzÞqi ðtÞ fi ðzÞzi ðtÞ @a0 þ a1 fi ðzÞqi ðtÞ

i

!2

i

þ a6

The space/time decoupled form of D(z,t) and h(z,t) are substituted into Eq. (5) to have Eq. (6),

þ a3

i

!3 13 X fi ðzÞqi ðtÞ A5dtdz=QT

i

i

(6)

3.1.1. Forward modeling The depth-of-discharge D(z,t) can be expressed in space/time decoupled form as Eq. (1),

Dðz; tÞ ¼

X∞

f ðzÞqi ðtÞ i¼1 i

(1)

Where fi(z) are spatial basis functions and qi(t) are the corresponding time variables. The potential difference between electrodes is assumed to have same spatial basis functions as D(z,t), and can be denoted as follows,

hðz; tÞ ¼ Vp ðz; tÞ Vn ðz; tÞ ¼

X∞

f ðzÞzi ðtÞ i¼1 i

(2)

Where Vp(z,t) and Vn(z,t) are potentials of positive and negative electrode respectively, their formulations are summarized in Table A1 of Appendix A. After the space basis functions are identiﬁed by the KL method, the dimension reduction is required to reduce the inﬁnite form of ffi ðzÞg∞ i¼1 in Eq. (1) and Eq. (2) intoffi ðzÞgki¼1 , where k is determined when x > 0.99 in Eq. (B.15). After substitution Eqs. (A.4) and (A.5) into Eq. (A.3) the intraelectrodes potential can be derived as Eq. (3),

Jðz; tÞ ¼ a4 þ a5 Dðz; tÞ þ a6 Dðz; tÞ2 h i hðz; tÞ a0 þ a1 Dðz; tÞ þ a2 Dðz; tÞ2 þ a3 Dðz; tÞ3 (3)

The basic idea of the KL approach is to ﬁnd modes that represent the dominant character of the system. The maximum number of nonzero eigenvalues is n N and the KL expansion yields empirical Eigen functions (EEFs) fi(z) in a descent order manner according to the magnitude of the corresponding eigenvalues. The ratio of the sums of k largest eigenvalues to the total sum gives the fraction of the variance retained in the k-dimensional space deﬁned by the Eigen-function associated with those eigenvalues. This ratio can be mathematically described as Eq. (B.17) in the paper. Typically, a value of more than 0.9 is often taken for h. In this paper, the P Pn ratioh ¼ k¼1 i¼1 li = i¼1 li 0:99, this means the ﬁrst order of the EEF captures more than 99% dominant behavior of the distributed parameter system. This characteristic is determined by the specialty of physical property of the battery. Thus the ﬁrst order of fi(z) will be used in modeling of D(z,t) and shown in Fig. 3. For convenience, the termffi ðzÞgki¼1 can be simply written as f(x), fqi ðtÞgki¼1 as q(t) andfzi ðtÞgki¼1 as z(t) respectively. Substituting one-order form spatial basis functions into Eq. (6),

Z QT , U

vqðtÞ ¼ fðzÞ2 dz, vt

Z

i h fðzÞ a4 þ a5 ðfðzÞqðtÞÞ þ a6 ðfðzÞqðtÞÞ2

U

h fðzÞxðtÞ a0 a1 ðfðzÞqðtÞÞ

i a2 ðfðzÞqðtÞÞ2 a3 ðfðzÞqðtÞÞ3 dz (7)

264

M. Wang, H.-X. Li / Journal of Power Sources 301 (2016) 261e270

By reformulating the left hand side of the Eq. (7), Eq. (R8) as follows,

Z

The q(t) and z(t) in Eq. (7) are unknown variables. From Fig. 1, we can see that the terminal voltage vterminal is actually the boundary point measurement of the potential difference at electrodes. Then

vqðtÞ fðzÞ2 dz, vt

QT , U

Z h i h i a4 fðzÞ þ a5 fðzÞ2 qðtÞ þ a6 fðzÞ3 qðtÞ2 fðzÞxðtÞ a0 a1 ðfðzÞqðtÞÞ a2 ðfðzÞqðtÞÞ2 a3 ðfðzÞqðtÞÞ3 dz ¼ U

0 1 a4 fðzÞ2 xðtÞ a0 a4 fðzÞ a1 a4 fðzÞ2 qðtÞ a2 a4 fðzÞ3 qðtÞ2 a3 a4 fðzÞ4 qðtÞ3 Z B C B C ¼ B þa5 fðzÞ3 ,xðtÞqðtÞ a0 a5 fðzÞ2 qðtÞ a1 a5 fðzÞ3 qðtÞ2 a2 a5 fðzÞ4 qðtÞ3 a3 a5 fðzÞ5 qðtÞ4 Cdz @ A U þa6 fðzÞ4 ,xðtÞqðtÞ2 a0 a6 fðzÞ3 qðtÞ2 a1 a6 fðzÞ4 qðtÞ2 a2 a6 fðzÞ5 qðtÞ4 a3 a6 fðzÞ6 qðtÞ5 Z Z Z Z Z ¼ a4 fðzÞ2 dz ,xðtÞ þ a5 fðzÞ3 dz ,xðtÞqðtÞ þ a6 fðzÞ5 dz ,xðtÞqðtÞ2 a0 a4 fðzÞdz fðzÞ2 ða1 a4 þ a1 a4 Þdz ,qðtÞ U

U

U

U

(8)

U

Z Z Z fðzÞ3 ða2 a4 þ a1 a5 þ a0 a6 Þdz ,qðtÞ2 fðzÞ4 ða3 a4 þ a2 a5 þ a1 a6 Þdz ,qðtÞ3 fðzÞ5 ða3 a5 þ a2 a6 Þdz ,qðtÞ4 U

U

Z fðzÞ6 a3 a6 dz ,qðtÞ5

U

U

its time variable z(t) can be described as a polynomial function of vterminal as below.

The simpliﬁed form will be obtained as follows.

zðtÞ ¼ k0 þ k1 vterminal þ k2 v2terminal þ k3 v3terminal þ :::

vqðtÞ ¼ b0 þ b1 qðtÞ þ b2 qðtÞ2 þ b3 qðtÞ3 þ b4 qðtÞ4 þ b5 qðtÞ5 vt

Since the terminal voltage vterminal is measurable, by substituting Eq. (10) to Eq. (9), only one time variable needs to identify in Eq. (11),

2

þ g1 zðtÞ þ g2 qðtÞzðtÞ þ g3 qðtÞ zðtÞ ¼

5 X

bi ,qðtÞi þ

3 X

i¼0

gi ,½qðtÞði1Þ ,zðtÞ

(9)

i¼1

where R R R a0 a4 fðzÞdz fðzÞ2 ða1 a4 þa1 a4 Þdz fðzÞ3 ða2 a4 þa1 a5 þa0 a6 Þdz R R ; b1 ¼ U ; b2 ¼ U , b0 ¼ U R QT fðzÞ2 dz QT fðzÞ2 dz QT fðzÞ2 dz U U U R R R fðzÞ4 ða3 a4 þa2 a5 þa1 a6 Þdz fðzÞ5 ða3 a5 þa2 a6 Þdz fðzÞ6 a3 a6 dz R R ;b4 ¼ U ;b5 ¼ U R , b3 ¼ U QT fðzÞ2 dz QT fðzÞ2 dz QT fðzÞ2 dz U U R R R U 5 2 3 a4 fðzÞ dz a5 fðzÞ dz a6 fðzÞ dz ;g2 ¼ UR ;g3 ¼ UR . g1 ¼ UR 2 2 2 QT

U

fðzÞ dz

QT

U

fðzÞ dz

QT

U

fðzÞ dz

(10)

b q k ¼ f ðqk1 ; vterminal Þ

(11)

where the b q k denotes the variable estimated at time k. After synthesizing the identiﬁed b q k and the spatial function fi(z), the spatiotemporal variable D(z,t) can be constructed.

3.1.2. Backward modeling Since the depth of discharge D(z,t) is an unmeasurable variable, so the modeling is needed to identify qi(t), which is the temporal

Fig. 3. (a) Spatial basis function f(z), (b) Spatial basis function 4(z), (c) time variable a(t)

M. Wang, H.-X. Li / Journal of Power Sources 301 (2016) 261e270 20

5

(a)

15

raw value noise added value

265

(b)

raw value noise added value

4.8

10 4.6

5

Value

Value

0 -5 -10 -15

4.4 4.2 4

-20 3.8

-25 -30

0

20

40

60

80

3.6

100

0

20

40

60

Time

80

100

Time

Fig. 4. (a) Applied Current iapp proﬁles and measurements, (b) Terminal voltage vterminal measurements.

3

2.98 Real value Estimation

(a)

2.96

Real value Estimation

(b)

2.98 2.96

2.94

2.94 Value

Value

2.92 2.9

2.92 2.9

2.88

2.88

2.86 2.84

2.86

0

10

20

30

40

50

60

70

80

2.84

90

0

10

20

30

40

50

60

70

80

90

Time

Time

Fig. 5. Estimation of b qðtÞ for D(z,t) under (a) correct initial condition, (b) incorrect initial condition.

variable of D(z,t). By combining qi(t) with the spatial basis function f(z), the depth of discharge D(z,t) can be constructed. As shown in Fig. 2, identiﬁcation of qi(t) requires two real-time measurement iapp and Vterminal with their relationship clearly shown in Fig. 1. Without considering the spatial effects, the depth of discharge (DOD) can be calculated by coulomb counting method as [22].

Zt DðtÞ ¼ Dðt ¼ 0Þ þ

iapp dt

(12)

0

It is obvious that the temporal effects of DOD are the cumulative Rt effect of applied currentIðtÞ ¼ 0 iapp dt. Since I(t) is known due to measurement of iapp, so q(t) can be identiﬁed with these two known measurement iapp and vterminal under the forward-backward structure in Fig. 2. the relationship between I(t) and q(t) is identiﬁed in data based approach as Eq. (13).

bIðtÞ ¼ hðqðtÞÞ ¼ p þ p qðtÞ þ p qðtÞ2 þ p qðtÞ3 þ ::: 0 1 2 3

(13)

The LASSO method can be used to reduce complexity of the identiﬁcation matrix and make the parameter estimation more efﬁcient. 3.1.3 Online parameter estimation It is always difﬁcult to ﬁgure out model structure, like number of parameters, in the regression method. LASSO, which is short for least absolute shrinkage and selection operator [23], can automatically select the important parameters by shrinking the trivial parameters to zero. So, the LASSO is used to enhance the extended Kalman ﬁlter (see Appendix C). Model uncertainties such as environmental disturbance and noises will be simulated by adding disturbance terms wk and vk, which are assumed to follow Gaussian distributions N(0,Pk) and N(0,Rk) respectively. The parameter identiﬁcation in the Kalman ﬁlter can be expressed in a discrete form as follows.

qkþ1 ¼ f ðqk ; vterminal Þ þ wk

(14)

Ik ¼ hðqk Þ þ vk

(15)

Table 1 Comparison of RMSE ofb qðtÞ under different conditions.

RMSE

Correct initial condition without noise

Correct initial condition with noise

Incorrect initial value with noise

0.0050

0.0060

0.0173

266

M. Wang, H.-X. Li / Journal of Power Sources 301 (2016) 261e270

0.126

0.13

0.124

0.128 0.126

0.122 0.1

0.1

0.12

0.124 0.122

0.118

0.12

0.116

0.05

0.05

0.118

0.114

0.116

0.112 0

0

0

0.02 0.04 0.06 0.08

0.114 0

0.02 0.04 0.06 0.08

(b)

(a) 0.15

0.15

4.592

4.131

4.591

4.13

4.59

0.1

4.129

0.1

4.128

y

y

4.589 4.588 0.05

4.587

4.127

0.05

4.126

4.586 0

0 0.02 0.04 0.06 0.08 (c)x

4.585

4.125 0

0 0.02 0.04 0.06 0.08 x (d) -3

0.15

x 10 0

0.15 -0.005

-2

-0.01 -0.015

0.1

-0.02

-6

-0.025

-8

-0.03

0.05

-4

0.1

0.05

-10

-0.035 -12

-0.04 0

-0.045 0

0.02 0.04 0.06 0.08

0

-14 0

0.02 0.04 0.06 0.08

(f)

(e)

Fig. 6. Predictions of critical variable at 40s and 80s: D(z,t) at (a) 40s and (b) 80s, Vp(z,t) at (c) 40s and (d) 80s, Vn(z,t) at (e) 40s and (f) 80s.

With the online estimated{b q k } and ofﬂine identiﬁed space basis function {f(z)}, the spatiotemporal model of D(z,t) can be constructed and updated in the proposed forward-backward learning structure. 3.2. Modeling of potentials

Since the terminal voltage vterminal is the potential difference that can be measured at the boundary, so the extended Kalman ﬁlter is not needed. The potential of positive electrode can be expressed in the following form of time/space separation.

Vp ðz; tÞ ¼

In this section, space-time separation based online models are developed for the potentials of positive and negative electrode.

k X

4i ðzÞai ðtÞ

(16)

i¼1

Table 2 Root Mean square error of critical variables with incorrect initial value.

RMSE

Depth-of-discharge D(z,t)

Potential of positive electrode Vp(z,t)

Potential of negative electrode Vn(z,t)

0.0080

0.0129

0.0095

M. Wang, H.-X. Li / Journal of Power Sources 301 (2016) 261e270

wheref4i ðzÞgki¼1 are spatial temporal basis functions derived by KL method. Similarly in the previous section, the ﬁrst order off4i ðzÞgki¼1 is sufﬁcient to cover dominant dynamics of Vp(z,t) measured in Eq. (B.8). Thus, spatial functionsf4i ðzÞgki¼1 can be simpliﬁed into 4(z) and the corresponding temporal variablesfai ðtÞgki¼1 into a(t). The terminal voltage vterminal can be expressed as the potential at the measurement point in Eq. (17),

Vp ðz; tÞz¼z

pmea

¼ 4 zpmea aðtÞ

(17)

Where zpmea is the measurement point at the tab where the terminal voltage can be measured by Voltmeter and the potential of negative electrode at tab is usually set as zero. Since vterminal and 4(zpmea) are known, thus a(t) can be derived analytically as

aðtÞ ¼ vterminal =4 zpmea

(18)

The modeling of a(t) is shown in (c) of Fig. 3 that matches well with the real measurement. From Eq. (2), the potential of negative electrode can be expressed as.

Vn ðz; tÞ ¼ Vp ðz; tÞ hðz; tÞ

(19)

In traditional approaches, vterminal is used as the index to prevent over-voltage (during discharging) and under-voltage (during charging). However, vterminal is the temporal variable measured at the ﬁxed location and the voltage over the battery space is the spatiotemporal variable that cannot be measured online but can be simulated using the physical model. With the spatial distribution of potential difference obtained, we can propose a new indicator to prevent the over-voltage and under-voltage over entire working space as given below.

upper limit indicator ¼ maxðhðx; tÞÞ

(20)

lower limit indicator ¼ minðhðx; tÞÞ

(21)

4. Model validations The proposed spatiotemporal models will be evaluated in the real-time battery test experiment. The current proﬁles are produced by the UDDS test of hybrid electrical vehicles. A Gaussian noise with 10% signal-to-noise ratio per sample in dB is added to

267

the applied current iapp and terminal voltage vterminal as shown in Fig. 4. The spatial basis function is identiﬁed ofﬂine as shown in Fig. 3. The estimated temporal variable q(t) in Eq. (1) will be compared between the analytical modeling and the physics modeling in COMSOL. Since the initial value of DOD is usually unknown, so initial value selected for b q k may not be correct. The proposed modeling method can have a good prediction performance for both correct and incorrect initial value of b q k as shown in Fig. 5. The root mean square errors (RMSE) under different conditions are summarized in Table 1. Since both DOD and potential are spatiotemporal variables, the predicted spatial distributions are presented in Fig. 6, with corresponding errors summarized in Table 2. The total computation time of the proposed model is within 1 s, while the physics-based model in COMSOL cost about 20 min. Thus the spatiotemporal model developed is efﬁcient and reliable in estimating the internal states of lithium-ion batteries. The voltage measurement vterminal cannot represent the potentials of the battery over entire space. The proposed spatiotemporal modeling can provide internal information over entire working space. To better compare the value of vterminal and h(z,t), we plot yterminal over the working space. The comparison of vterminal and potential difference h(z,t) is plotted in Fig. 7. The Fig. 7 shows that the largest potential is bigger than vterminal and lowest potential difference is smaller than vterminal. Thus the extreme values of “potential difference distribution” can be used as the indicators to improve the battery protection. The predicted h(z,t) in Fig. 7, demonstrates that Eqs. (20) and (21) can be better indicators to prevent over-voltage and under-voltage. 5. Conclusion A spatiotemporal modeling method has been developed to online estimate the internal states of lithium-ion battery. The spatial functions are identiﬁed using KL method. Then the temporal variables are identiﬁed with help of the extended Kalman ﬁlter, and optimized by the LASSO adaptation method. The space-time separation based analytical models developed can predict the dynamic performance over entire working space in a very small computation time compared with the physics-based model in COMSOL. Thus, the developed model can work for online prediction of the battery state. The real experiment demonstrates the feasibility and effectiveness of the proposed modeling method.

Fig. 7. Potential differences of the cell and terminal voltage (vterminal) at time 65s.

268

M. Wang, H.-X. Li / Journal of Power Sources 301 (2016) 261e270

Acknowledgment The work presented in the paper is partially supported by a GRF project from RGC of Hong Kong SAR (CityU: 11205615).

numerical solutions in space-time dimensions. The dominant spatial feature of these snapshots can be described as a series of spatial basis functionsffi ðxÞgni¼1 . They are mathematically derived by optimizing the following objective function [24,25],

D

Appendix A. Two-dimensional model for battery

max f

jðy; fÞj2

E (B.1)

kfk2

The two-dimension model for Lithium-ion batteries only considers the Ohm's law and the charge conservation on the electrodes [6]. As shown in Fig. 1, the applied current ﬂows into the battery though one tab and ﬂows out though another tab. Due to the ohm resistance of the electrodes, the current is not distributed evenly over space. Based on the continuity of current in the electrodes, the following equations can be derived,

Where (,,,) denotes the inner product, 〈,〉 denote the ensemble average and k,kis the induced norm kfk¼(f,f)1/2. Then the problem turns to maximize 〈j(y,f)j2〉 subject to kfk2 ¼ 1. The constraint kfk2 ¼ 1 is imposed to ensure f is unique. Then the corresponding Lagrangian function can be established as Eq. (B.2).

! V, i p ðz; tÞ Jðz; tÞ ¼ 0 z2Up ;

E D J ¼ ðfi ðxÞ; yðx; tÞÞ2 lððf; fÞ 1Þ

¼0

! V, i n ðz; tÞ Jðz; tÞ

z2Un

(A.1)

! ! where i p ðz; tÞ and i n ðz; tÞ are the linear current density vectors in the electrodes and J(z,t) is the inter-electrode current density. According to Ohm's law, the following equations can be derived,

! ! i p ðz; tÞ ¼ 1 r VVp ðz; tÞ z2Up ; i n ðz; tÞ n

A necessary condition for extrema is that the functional derivative vanish for all variation f þ dj,

d Jðf þ djÞ ¼0 dd d¼0

(B.3)

From Eq. (B.2), we have

p

¼ 1=r VVn ðz; tÞ z2Un

(B.2)

(A.2)

d d Jðf þ djÞ ¼ ½〈ðy; f þ djÞðf þ dj; yÞ〉 lðf þ dj; f þ djÞjd¼0 dd dd d¼0

where rp and rn are resistance at positive and negative electrodes, Vp(z,t) and Vn(z,t) are potentials of positive and negative electrodes respectively, as provided in the Table A1. The inter-electrode current density J(z,t) can be written as a polarization expression of potential difference Vp(z,t)Vn(z,t) [6],

Interchanging the order of the averaging operation and the inner product, the quantity in the brackets can be written as Eq. (B.5)

Jðz; tÞ ¼ Yðz; tÞ Vp ðz; tÞ Vn ðz; tÞ Uðz; tÞ

ð〈ðf; yÞy〉; jÞ lðf; jÞ

(A.3)

¼ 2Re½〈ðy; jÞðf; yÞ〉 lðf; yÞ ¼ 0 (B.4)

(B.5)

Where Y(z,t) and U(z,t) are ﬁtting variables that are functions of depth of discharge D(z,t) as shown in Eqs. (A.4) and (A.5) respectively, with parameters a0 a6 to be determined from experiments.

To further simpliﬁed the quantity in the brackets, the linear operator R is deﬁned by Rf¼〈(f,y)y〉, then Eq. (B.5) can be simpliﬁed as Eq. (B.6),

Uðz; tÞ ¼ a0 þ a1 Dðz; tÞ þ a2 Dðz; tÞ2 þ a3 Dðz; tÞ3

(A.4)

ð

Yðz; tÞ ¼ a4 þ a5 Dðz; tÞ þ a6 Dðz; tÞ2

(A.5)

Finally, since j is an arbitrary variation, the condition reduced to the eigenvalue problem

According to the reference [6], the depth of discharge is deﬁned as follows, with QT as the total capacity per unit area of the electrodes.

Dðz; tÞ ¼

Jðz; tÞdt=QT

(A.6)

(B.7)

For the deﬁnition of ensemble average, (R7) becomes

Z

Zt

(B.6)

〈yðx; kÞyðx0 ; kÞ〉fj ðx0 Þdx0

(B.8)

U

0

By deﬁning the two-point correlation function,

Table A1 Potential variables of the 2D model Domain

Potential variable

Governing equation

Positive electrode Up

Vp(z,t)

V2Vp(z,t) ¼ rpJ(z,t)

Boundary condition at tabs vV ðz;tÞ iapp ¼ rLp p! ; vterminal ¼ Vp ðz; tÞz¼z Vn ðz; tÞjz¼znmea pmea vn

Negative electrodeUn

Vn(z,t)

V2Vn(z,t) ¼rnJ(z,t)

Vn ðz; tÞjz¼zn ¼ 0

z¼zp

M 1 X yðx; kÞyðx0 ; kÞ M

ve decomposition method Appendix B. KarhuneneLoe

Kðx; x0 Þ ¼ 〈yðx; kÞyðx0 ; kÞ〉 ¼

The interested spatiotemporal observationsfyðxi ; tÞgN;M are i¼1;t¼1 generally called snapshots. They are uniformly sampled from the

Therefore, the optimal condition Eq. (B.8) reduces to eigenvector/eigenvalue problem of the integral equation:

(B.9)

k¼1

M. Wang, H.-X. Li / Journal of Power Sources 301 (2016) 261e270

Z

269

with covariance

Kðx; x0 Þfj ðx0 Þdx0 ¼ lj fj ðxÞ

(B.10)

U

T Pkjk1 ¼ Fk1 Pk1 Fk1 þ Qk1

By the methods of SchmidteHilbert technique or the method of snapshots, fi(x) can be expressed as a linear combination of the snapshots,

where the Jacobian Fk1 is the derivative of function f with respect to x at the value ofb x k1jk1 . In the innovation phase, the innovated state is estimated by

M X

fi ðxÞ ¼

git yðx; tÞ

(B.11)

b x kjk1 þ Kk zkþ1 h b x kjk ¼ b x kjk1

(C.4)

(C.5)

t¼1

with covariance

Substituting Eq. (B.13) and Eq. (B.13) to Eq. (B.10) yields,

Z

M M M X X 1X yðx; tÞyðc; tÞ git yðc; tÞdc ¼ li git yðx; tÞ U M t¼1 t¼1 t¼1

Pkjk ¼ Pkjk1 Kk Sk KkT (B.12)

The eigenvalue problem of Eq. (B.5) can be equivalently written as

Cgi ¼ li gi

(B.13)

where C is the two-point temporal correlation function is deﬁned in following equation

Ctk ¼

1 L

Z yðx; tÞyðx; tÞdx

Xk

The innovation covariance Sk and the near-optimal Kalman gain Kk are given by

Sk ¼ Hk Pkjk1 HkΤ þ Rk

(C.7)

Kk ¼ Pkjk1 HkΤ S1 k

(C.8)

where the Jacobian Hk is the derivative of function h with respect to x at valueb x kjk1 .

(B.14) Appendix D. Introduction of LASSO

U

The above matrix Ctk is symmetric and positive deﬁnite. Eigenvectors gi¼[gi1,...,giL] corresponding ith eigenvalues can be derived by the standard method of eigenvalue decomposition. The above eigenvector of matrix eigenvalue in Eq. (B.15) is then substitute into Eq. (B.13) to generate the spatial basis functions fi(x). Arranging the eigenvalues l1>l2>l3/>ln in decent order, and the ratio between cumulative sums of ﬁrst k eigenvalues to the total can be determined as Eq. (B.17).

x¼

(C.6)

.X n

l i¼1 i

l i¼1 i

(B.15)

For x, a value of more than 0.99 is often taken which means the 99% dominance feature has been captured by the k-dimension spatial basis functions.

LASSO is short for least absolute shrinkage and selection operator [23]. The main idea of LASSO is that it adds the L2-norm of the parameter vector to the objective function and maintains it within a preset value. The objective function for parameter (q) estimation of LASSO is

2 1 b q ¼ argmin y FT q þ lkqk1 2

P where y is the output, F is the input features and kqk1 ¼ D i¼1 jqi j is the L1-norm, l is a hyper parameter that controls its inﬂuence. To optimize the objective function, q can be written as a difference between two vectors with positive entries

q ¼ qþ q

Extended Kalman ﬁlter is the nonlinear version of Kalman ﬁlter which deals with nonlinear state estimation problem under the noisy measurements [20,21]. The problem can be formulated as follows.

xk ¼ f ðxk1 ; uk1 Þ þ wk1

(C.1)

zk ¼ hðxk Þ þ vk

(C.2)

Thus the objective function can be rewritten as

2 X 1 b qþ þ q q ¼ argmin y FT qþ q þ l i i 2 qþ ;q i þ

Where wk1 and vk are assumes to be known zero-mean Gaussian process noise with wk~N(0,Q) and vk~N(0,R), uk1 is the control vector. The function f predicts system states xk from previous estimate xk1 and function h provides the expected measurement zk. The Jacobian matrix is evaluated at each time step. The implementation of the EKF includes prediction phase and innovation phase. In the prediction phase, the system state is predicted by

b x kjk1 ¼ f b x k1jk1 ; uk1

(D.2)

s:t:qþ 0; q 0

Appendix C. Extended Kalman Filter

(C.3)

(D.1)

s:t:q 0; q 0

(D.3)

qþ , the objective function can be written in a q standard form of quadratic program.

By deﬁning x ¼

1 min xT Hx þ f T x 2

(D.4)

FFT f ¼ l1 Fy , and 1 is the vector of FFT T T Fy FF FF ones. Thus, the parameter q is derived by solving the above optimization problem. Therefore L1-norm forces parameters of trivial feature to zero so that the complexity of the model can be reduced. whereH ¼

270

M. Wang, H.-X. Li / Journal of Power Sources 301 (2016) 261e270

References [1] M. Doyle, T.F. Fuller, J. Newman, J. Electrochem. Soc. 140 (1993) 1526e1533. [2] P. Ramadass, B. Haran, P.M. Gomadam, R. White, B.N. Popov, J. Electrochem. Soc. 151 (2004). A196. [3] S. Santhanagopalan, Q. Guo, P. Ramadass, R.E. White, J. Power Sources 156 (2006) 620e628. [4] K.H. Kwon, C.B. Shin, T.H. Kang, C.-S. Kim, J. Power Sources 163 (2006) 151e157. [5] G.-H. Kim, K. Smith, K.-J. Lee, S. Santhanagopalan, A. Pesaran, J. Electrochem. Soc. 158 (2011). A955. [6] L. Cai, R.E. White, J. Electrochem. Soc. 156 (2009). A154. [7] K.A. Smith, C.D. Rahn, C.-Y. Wang, J. Dyn. Syst. Meas. Control 130 (2008) 011012. [8] V.R. Subramanian, V.D. Diwakar, D. Tapriyal, J. Electrochem. Soc. 152 (2005). A2002. [9] T.-S. Dao, C.P. Vyasarayani, J. McPhee, J. Power Sources 198 (2012) 329e337.

[10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25]

G.L. Plett, J. Power Sources 134 (2004) 277e292. J. Lee, O. Nam, B.H. Cho, J. Power Sources 174 (2007) 9e15. R. Xiong, X. Gong, C.C. Mi, F. Sun, J. Power Sources 243 (2013) 805e816. G.L. Plett, J. Power Sources 161 (2006) 1369e1384. G.L. Plett, J. Power Sources 161 (2006) 1356e1368. S. Santhanagopalan, R.E. White, Int. J. Energy Res. 34 (2010) 152e163. F. Sun, X. Hu, Y. Zou, S. Li, Energy 36 (2011) 3531e3540. U.S. Kim, C.B. Shin, C.-S. Kim, J. Power Sources 189 (2009) 841e846. H. Gu, J. Electrochem. Soc. 130 (1983) 1459e1464. H.-X. Li, C. Qi, J. Process Control 20 (2010) 891e901. S.S. Haykin, S.S. Haykin, S.S. Haykin, Kalman Filtering and Neural Networks, Wiley Online Library, 2001. M. Hoshiya, E. Saito, J. Eng. Mech. 110 (1984) 1757e1770. K.S. Ng, C.-S. Moo, Y.-P. Chen, Y.-C. Hsieh, Appl. Energy 86 (2009) 1506e1511. R. Tibshirani, J. R. Stat. Soc. Ser. B Methodol. (1996) 267e288. H.-X. Li, C. Qi, Y. Yu, J. Process Control 19 (2009) 1126e1142. C. Qi, H.-X. Li, Chem. Eng. Sci. 64 (2009) 4164e4170.