Order Reduction of Lithium-Ion Battery Model Based on Solid State Diffusion Dynamics via Large Scale Systems Theory

A simple and accurate model is essential for the model-based control design. The 1D electrode-averaged battery model maintains the prediction accuracy and depth of physical interpretations. The most critical and complicated dynamics in the system is solid phase concentration diffusion, which has a direct impact on battery status. The ﬁnite difference method based on uneven discretization beneﬁts to reduce the complexity. Then, the simpliﬁcation theory of large scale systems is applied to retain the dominant variable and make the model more compact. The simpliﬁcation process meets the requirements of fast computation and retention of electrochemical diffusion dynamics. This approach was experimentally veriﬁed under 0.2C to 2Cdischarge tests. The simulation output converges to the real terminal voltage within an error of ± 0.1V. The simpliﬁed model enjoys advantages including high accuracy, low computational cost, and therefore is suitable for deployment and use in real-time embedded battery management systems (BMSs). © The Author(s) 2016. by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, 10.1149/2.1251607jes]

Recently, Lithium-ion batteries have been widely used in the electric vehicles (EVs), hybrid electric vehicles (HEVs) and plug-in electric vehicles (PHEVs), 1 thanks to its high energy density, no memory effect, and low self-discharge when not in use. However, Lithium-ion batteries are exposed to the potential risk of damage or explosion under extreme operating conditions, such as high temperatures and over-charging/discharging. 2 Hence, the design of a sophisticated BMS is necessary to ensure safe battery operations under wide variety of driving conditions. 3 Most BMS designs are based on simple equivalent circuit models (ECM) because they have fast computation and relatively good battery state of charge (SOC) prediction capability. 1 However, these models lack the predictiveness of lithium-ion electrochemical diffusion dynamics, so the predicted responses cannot provide sufficient accuracy over various driving conditions. Unlike equivalent circuit models, electrochemical models describe the chemical processes that take place in battery in great detail, which makes them the most accurate of battery models. However, the highly detailed description results in considerable computational complexity, and it may take hours to simulate a charge-discharge cycle of a detailed battery model if no model reduction approach is used. 4 Literatures on electrochemical modeling of batteries is quite extensive, including both full order and simplified models. Newman and Tiedemann presented a first electrochemical approach to porous electrodes modeling for battery applications in Ref. 5. Then, a fullorder model was developed by Doyle et al. in Ref. 6, they built a lithium anode/solid polymer separator/insertion cathode cell model using concentrated solution theory. The model is general enough to include a wide range of polymeric separator materials, lithium salts, and composite insertion cathodes. The same approach was later expanded to the dual lithium ion insertion (rocking-chair) cell in Ref. 7. Based on this model, a micro-macroscopic coupled model was introduced by Wang et al. in Ref. 9. It incorporates solid-state physics and interface chemistry and can be adapted to a wide range of active materials and electrolyte solutions. In fact, it was presented for Ni-MH batteries at the beginning. 10 Then, it was expanded to lithium-ion batteries in Ref. 11 where the thermal behavior was also described. In addition, Ramadass et al. in Ref. 12 took into account the decay in capacity; Sikha et al. in Ref. 13 included the change in the porosity of electrode material using a pseudo-two-dimensional approach, and further improved in the follow-up studies. [15][16][17] Full order models can predict concentration profiles across the electrodes and electrolyte, but their long computational time and the z E-mail: xuxing@mail.ujs.edu.cn extensive complexity prevent their application in control design and real-time on-board estimator. In Refs. 8,14 several approximations between the input and output of dynamical system are introduced to reduce the system order and complexity. To reduce the model order for fast simulation, a residue grouping approach was proposed by Smith et al. 18,19 This approach truncated and grouped the eigenvalues of the diffusion dynamics system. Another approach is the model reductions using electrode-averaged lithium-ion battery model proposed by Di Domenico et al. [20][21][22] This model assumes that the concentration distribution along the length of electrodes and separator are constant under mild charging/discharging cases. The model was further modified by applying uneven discretization of the ordinary differential equations (ODEs) in Refs. 1,2, to enable fast computation for model-based control design.
In this work, the simplification theory of large scale systems is introduced to make the electrode-averaged model more compact for the requirements of fast computation. The paper is organized as follows. First, the electrode-averaged battery model is extracted from the recent electrochemical model researches for lithium-ion batteries. Then, the finite difference method is used to solve the solid concentration diffusion equation, where two different discretization approaches are compared. The uneven discretization turns out to be benefit to reduce model complexity while maintaining the prediction accuracy. Next, the simplification theory of large scale systems is applied, based on uneven discretization, to make the model more compact. Experiments are followed to validate the model accuracy. Finally, summary and conclusions are presented.

Electrode-Averaged Model
The electrode-averaged battery model proposed by Di Domenico et al. [20][21][22] is derived from the micro-macroscopic coupled model introduced by Wang et al. in Ref. 9. The micro-macroscopic coupled lithium ion cell model considers dynamics along only one axis(the X-axis), and both macroscopic and microscopic physics have been considered. This approximation depicted in Figure 1 consists of three domains-the negative composite electrode (with Li x C 6 active material), separator, and positive composite electrode (consist of Li y Mn 2 O 4 , Li y CoO 2 , Li y NiO 2 , or some combination of metal oxides).
Almost all of the electrochemical Li-ion battery models are based on a set of differential algebraic equations that describe the timedependent electrochemical phenomena between electrodes. 4 Even after assuming fixed electrolyte concentration, these models are of high order and complexity. 20 The complete set of equations are summarized in Table I, which describe the battery system with four quantities, i.e., Table I. Summary of model equations for a Li-ion battery.

Model equations Boundary conditions
Electrolyte phase potential Solid phase potential Species in solid phase ∂r ∂cs ∂r | r =0 = 0 ∂cs ∂r | r =Rs = − j Li as F Ds [5] Butler-Volmer current density ) − R f I [7] solid and electrolyte concentration (c e , c s ) and solid and electrolyte potential (ϕ e , ϕ s ). The battery parameters used for constructing the electrodeaveraged model are listed in Table II. The parameters are cited from Reference 20 and revised by providing fits to experimental data.
In the electrode-averaged model, the solid concentration distribution along electrode is neglected, and the concentration in electrolyte is assumed to be constant. This approximation is similar but higher order than the single-particle electrode model in Ref. 24. Although this simplified model may reduce the prediction accuracy under highcurrent discharging and charging conditions, it can be useful in control and estimation applications, and the important diffusion dynamics in solid particles are still captured.
In accordance with the average solid concentration, the spatial dependence of Butler-Volmer current is assumed as a constant valuē j Li n , which is evaluated from the spatial integral.  The battery terminal voltage Eq. 7, using average values at the positive and negative electrode instead of the boundary values, can be rewritten as V (t) =η p −η n + φ e, p −φ e,n + U p c se, p − U n c se,n − R f I [9] Finally, the battery voltage can be written as a function of battery current and of the average solid concentration, where ξ n =¯j In order to obtain the battery terminal voltage, it is essential to figure out the Lithium-ion concentration profiles in electrode particles, which can be evaluated from Eq. 5. The partial differential equation is based upon the Fick's second law under spherical coordinates. It cannot be applied to the control system if no simplification approach is used. Then, our work focuses on solving this issue via the finite difference method and the simplification theory of large scale systems.

The Finite Difference Method
The diffusion equation can be expressed as a set of ordinary differential equations by using the finite difference method for spatial variable r. Then, we can get a state-space formulation of the solid concentration diffusion equations, which is primary for the battery control oriented model. The sphere radius is discretized into r = (r 1 ,r 2 ,. . . ..,r Mr−1 ). There are two different discretization approaches: even discretization and uneven discretization, depicted in Figure 2.

A1432
Journal of The Electrochemical Society, 163 (7) A1429-A1441 (2016) Even discretization.-By dividing the solid material sphere radius in M r same intervals, with M r large enough, of size r = R s / M r , the first and the second derivative of the diffusion equation can be approximated with finite difference. So the Eq. 5 becomes [11] For boundary conditions the same procedure was applied. Hence, a state-space formulation of the solid concentration diffusion equations can be obtain. The model input I governs boundary conditions through the Butler-Volmer current. where However, the number of discretization steps M r must be large enough to satisfy accuracy requirements. This way, the state-space formulation would be a large scale system, which cannot be applied to control systems. Besides, the state matrix A 1 turns out to be close to singular or badly scaled in the MATLAB operation, where RCOND is close to zero and rank( Since each interval has the same length, which makes the elements in matrix A 1 having some certain regularity, as a result, A 1 is quasisingular.
Generally, using even discretization, we can obtain the state-space equation of lithium-ions distribution in electrode particles. But the state matrix A 1 is a large quasi-singular matrix, so we cannot apply the simplification methods for large system.
Uneven discretization.-Then, the uneven discretization method is used to change the relations between those elements in matrix A 1 . Each interval can be selected randomly, the dominant area should be discretized specifically, while the non-dominant area can be discretized roughly. This way, the elements in matrix A 2 is hard to have any regularity. In the literature, 1 the optimal uneven discretization results indicate that the radius close to surface must be finely discretized, while the discretization step width increases toward the center of particle. Hence, the difference between two adjacent r k and r k+1 should be smaller toward surface, which means the change rule of r k is similar to a parabolic. Then, we propose the function y = −(x − b) 2 + b 2 to create a series of r k , which meets the requirement, Figure 3 shows the uneven discretization results. However, in order to describe the concentration distribution more precisely, it is necessary to use some optimization methods to obtain the r k . In this paper, we pay more attention to the whole process of model reduction rather than the model accuracy, and the results of this uneven discretization is quite acceptable.
Through the same procedure described in the even discretization approach, a similar state-space formulation of solid concentration diffusion equations can be obtained. In this system, the matrix A 2 is no more close to singular and can be calculated in MATLAB without any warning, rank(A 2 ) = M r − 1. Then, we can use the simplification methods in the following part to build a more compact model suitable for model-based control design.
) unless CC License in place (see abstract Other model reduction approaches could show the similar eigenvalue distributions, which describe the dominant battery dynamics. However, the uneven discretization uniquely maintains the physical interpretation related to Li-ion diffusion dynamics. Simulation result.-The tridiagonal matrix algorithm (TDMA), also known as Thomas algorithm, is a simplified form of Gaussian elimination that can be used to solve tridiagonal systems of equations. It is an adaptation of LU decomposition idea to solution of a linear system in which the coefficient matrix A is populated with zeroes except for the diagonal, sub-diagonal and super-diagonal elements. Given the extreme sparsity of coefficient matrix, a considerable reduction in memory requirements and an increase in computational speed can be gained by storing the nonzero elements in three 1D vectors instead of saving the entire n × n coefficient matrix. The results of Thomas method can be seen as the most accurate numerical solution. In this part, three different methods of solving the partial differential equation has been compared for positive electrode, namely, Thomas method and transform the differential equation into a state space equation with even or uneven discretization. The numerical solution is implemented in MATLAB/SIMULINK.
The battery used in this research is Boston-Power's Swing 5300, which can offer high usable energy density combined with unmatched safety features and longer cycle life at broad operating temperatures. It is ideal for a wide range of applications. According to the specification, a charge/discharge test is implemented. The test procedures, shown in the top plot of Figure 4, consist of a 1C battery discharge till 2.7 V, an open circuit relaxation for 15 min, and a 1 C charge till 4.2 V followed by open circuit relaxation, the initial SOC is 100%, and the environment temperature is 25 ± 2 • C. In the simulation, we repeat this test to analyze the performance of each method.
The cathode Li-ion concentration profiles of three different methods are demonstrated in Figure 4. The solution of uneven discretization is more accurate than even discretization with less difference dimension, which benefits to simplify the large system.
Hence, using uneven discretization, we can get a high-precision, low-dimension state-space equation of the concentration distribution in solid particles.
For model-based control design, the state-space equation should be as simple as possible. However, as we can see in Figure 5, with the discretization number changes, the bias of concentration profiles is noticeable, with the number of discretization decreasing to 20, the concentration converges, and further decreases in the solution accuracy diminish. The uneven discretization number should be neither too small nor too large. Hence, uneven discretization with 50 steps is selected for the further research. The bias between uneven discretization with 50 steps and Thomas method is acceptable. Figure 6 illustrates the terminal voltage sensitivity to the number of uneven discretization steps. Most of the data are highly consistent with the result of Thomas method, except for the inflexions, and the discretization with 50 steps still turns out to be the best choice.

Simplification of Large Systems
The model simplification of large system is essentially a mathematical method. The basic idea of the simplification is to retain the dominant variable and correct steady-state behavior, and ignore the non-dominant variables that decay very fast. It does not change the physical meaning of original model variables, and does not introduce new parameters into the system. Through the optimization theory, the dominant modes selection principle can be established, which makes the least deviation between the simplified model and original model.

Model reduction process.-Let
Eq. 1 be the original model, the state variables are c si (i = 1,2, . . . ,n). Then, we can decompose the model into a mth-order (m < n) sub-model through the following simplification process. Firstly, select c s j ( j = 1, 2, · · · , m) as the basic state variables for simplified model, which should be the indispensable or vital variables among c si .ẋ = Ax + Bu [14] The state matrix A needs to be diagonalizable, which means there exists an invertible n × n matrix V, such that = V −1 AV . Then, let x = V z, so Eq. 14 can be rewritten as, The matrix is a diagonal matrix having eigenvaluesλ i (i = 1, 2, · · · , n) of matrix A down the main diagonal. The column vectros v i (i = 1, 2, · · · , n) of matrix V are the eigenvectors of matrix A corresponding to λ i . In order to obtain an accurate sub-model, a set of λ j ( j = 1, 2, · · · , m) need to be selected and retained. The reserved eigenvalues are related to the dominant modes of original model. See   "Selection of dominant eigenvalues " for more details. As the selected eigenvalues arranged side by side in a new matrix , the matrix G and eigenvectors in matrix V also need be adjusted in proper order. This way, the matrices can be expressed as several block matrices for convenience.
Furthermore, the Eq. 16 is decomposed into dominant mode and non-dominant mode from block multiplication. [18] Similarly, the original state variables x can be expressed as, The matrix x 1 consists of the basic state variable c s j ( j = 1, 2, · · · , m) selected at the beginning, while the remaining variables grouped in matrix x 2 can be neglected at following reduction. Then, the error matrix E is introduced to eliminate the non-dominant variables. See "Derivation of error matrix E" for more details. [20] Finally, the simplified sub-model can be expressed as, Selection of dominant eigenvalues.-The dominant eigenvalues (that is the dominant pole) are generally considered to be the points close to imaginary axis. Because the point which is far away from imaginary axis in left half plane has a smaller time constant, so its dynamic process disappeared quickly, and the dynamic response can be neglected in simplified model. The response amplitude also should be considered when determine the dominant mode of retention.
Taking the Laplace transform of Eq. 15 yields, ) unless CC License in place (see abstract With Eq. 19, we can get the transfer function of the basic state variable x 1 and analyze the dynamic response of system. When the input signal is step one u = 1/s, the output can be expressed as x 1 (s). And x 1 j (s) represents the corresponding response of state variable x j . Finally, the time-domain response x 1 j (t) is obtained with the Laplace transform.
Let c j,k denote the amplitude of response caused by λ k . Then, the dynamic response effect caused by λ k can be inflected through c j,k . This way, not only the position of the poles in the complex plane is considered, but the observability and controllability of system can be analyzed.
Every each of the basic state variables c s j ( j = 1, 2, · · · , m) has a corresponding relationship with the input, and there are c j,k (k = 1, 2, · · · , n) for each ' j − path'. Thus, the main criteria M k and S k are defined to evaluate the significance of each eigenvalues λ k .
Under the guidance of M k and S k , the dominant eigenvalues λ j ( j = 1, 2, · · · , m) can be picked out.

Derivation of error matrix E.-Taking
When the input signal is step one u = 1/s, the time-domain responses of dominant and non-dominant mode can be expressed as z 1 (t),z 2 (t).
[29] Obviously, the deviation ε i is related to the row vector e T i of error matrix E. Based on this, we propose a cost function as follow, The const in Eq. 29 is constrained to be zero, otherwise the integration does not converge.
[32] At this point, the derivation of error matrix E is transformed into a problem of restricted optimization for minimum. Then, Lagrange multiplier method is applied to strive for the Constrained extremum. The Lagrange multiplier isη i (i = 1, 2, · · · , n − m)and the objective function is expressed as, ) unless CC License in place (see abstract Finally, the error matrix E can be calculated by Eq. 33.
Application of the model reduction method.-After the uneven discretized with 50 steps, we can obtain a 49 × 49 non-singular matrix A 2 , which can predict the concentration distribution accurately. But the matrix A 2 is too large for control system. So, we try to use the simplification method of large system to retain the dominant variable and ignore the non-dominant variables.
Specific reduction steps are as follow: Step 1: Determine 5 dominant state variables and transform the matrix based on the order of dominant state variables:[ c s12 c s28 c s33 c s37 c s49 ].
Step 4: Transform the matrices V, G based on the new order of λ i and v i . Then, divide them into 4 and 2 parts respectively.
Step 5: Calculate the matrices S pq , T pq and E.
Step 6: Using error matrix E to construct the simplified model. The matrix A 2 is a tridiagonal matrix, its eigenvalues with the corresponding main criteria can be obtained by numerical calculation in step 3, the results are shown in Table III The simplification method is based upon the uneven discretization no = 50. Hence, we can also use the uneven discretization no = 5 with the same state variables:[ c s12 c s28 c s33 c s37 c s49 ]directly. Then, these two different methods are compared in the simulation. Figure 7 shows the comparison of uneven discretization and large system simplification method, both of them are 5 dimensions. The results of the simplification method are of higher accuracy than that of uneven discretization no = 5 and the uneven discretization with the same state variables. The close following of the Thomas method result implies that the electrochemical diffusion dynamics is also well captured. Besides, the SIMULINK running time of different methods are displayed in Table IV. Finally, the calculation time has been shorten by 97% compared with the uneven discretization no = 50.

Experiment and Discussion
After all the simplification, the complete mathematic model of Lithium-ion battery can be presented by the terminal voltage equation involving two simplified diffusion equations. In addition, it is difficult to measure the solid concentration inside batteries. So, it will not work to validate the simplified diffusion system directly. However, the battery terminal voltage is available. Then, we implement the same test profile shown in the top plot of Figure 4 with the battery test system shown in Figure 8. The experiment results are demonstrated in Figure  9. The overall variation tendency of the terminal voltage can be well captured. During the process of charging, the simulation output has obvious deviation with the experimental data. This is mainly because the model parameters in Table II are fitted through 1 C discharge data.
Then, we consider to analyze the accuracy of the simplified model under different discharge rate. A series of tests has been implemented, which consist of 0.2 C, 0.5 C,1 C, 1.5 C and 2 C discharge, both of them are started with 100% SOC and followed by a 15 min open circuit relaxation. The results shown in Figure 10 implies that the simplified model still has good consistency with the original electrodeaveraged model (thomas method), while the accuracy of the original model need to be improved. Generally, the simulation results and the experimental data have good agreement during the middle range of SOC. Fortunately, the batteries used in vehicles are also mainly operated in this range. Since the electrode-averaged model cannot cover all battery information as the full order model, the simulation data cannot track the actual values during 0.2 C long-term discharging. The influence of the battery temperature also has not been considered, which should be improved in the following researches. However, the simplification process propose in this paper is feasible to reduce the computational time significantly without losing too much information of the original model.

Conclusions
The finite difference method and simplification theory of large scale systems have been applied to construct the 1D electrodeaveraged Li-ion battery models for fast computation and accurate prediction of electrochemical diffusion dynamics. The work focuses on simplifying the solid Li-ion concentration diffusion to make it suitable for model-based control design. The finite difference method transforms the diffusion equation into a set of ordinary differential equations. Comparisons between two different discretization approaches indicated that uneven discretization benefits to enables significant reduction of model complexity while maintaining prediction accuracy. The simplification method reduces the state matrix to 5 × 5, which retains the dominant variable and correct steady-state behavior. The method is validated by a full 1C charge/discharge test and different C-rate discharge tests. The average deviation of simulation voltage in the middle range SOC is less than ±0.1 V. The model operation time is shortened by 97%. Hence, the model simplification method proposed in this paper is able to reduce the computational time significantly without any degradation of prediction accuracy and maintains the depth of physical interpretations.