An Adaptive Reduced-Order-Modeling Approach for Simulating Real-Time Performances of Li-Ion Battery Systems

In this work, a reduced-order-modeling(ROM) approach is introduced to simulate a battery system with multiple cells. The ROM is described by only two state variable equations and several explicit expressions. While running a stand-alone pseudo-2D (P2D) model by the nonlinear state variable modeling (NSVM) algorithm at a speciﬁed current, several parameters can be estimated through the NSVM procedure; and with these parameters, the cell voltage under additional current perturbations can be computed approximately by the ROM. This modeling scheme was tested for a 24P battery module in a co-simulation process: the stand-alone P2D model is solved at the battery-average current to estimate certain parameters to share with its ROM, then the ROM is implemented to evaluate the voltage response of each cell at its respective cell current. The ROM shows excellent accuracy and a much higher speed than the full-order-model with P2D models used for each cell.

The battery management system (BMS) plays a crucial role in electrical vehicles (EV). During a standard vehicle development cycle, in-loop simulation is a useful technique to test and verify complex real-time embedded systems. [1][2][3] To conduct in-loop simulation, the embedded system to be tested is connected to a mathematical representation of all related dynamic systems, because this arrangement lowers cost and improves safety as compared with using the real system. The mathematical representation interacting with embedded system is referred as "plant simulation". For the BMS software development, the plant simulation is performed by system-level battery models while the corresponding embedded system includes certain battery control and management algorithms. Therefore an accurate plant model for the battery system is vital for a successful BMS design. Besides accuracy, the plant battery model is also expected to work under high frequency inputs with noise and generate real-time output signals at low computational costs.
The pseudo-2D (P2D) model developed by Fueler et al. 4 is a physical model that can be used to provide reliable predictions for the transport and electrochemical phenomena in Li-ion cells. However, it is still difficult to use the P2D model in plant simulation. The first challenge is the numerical issues in solving the nonlinear partialdifferential-equations (PDEs) contained in the P2D model at highfrequency inputs. During practical driving modes, the power requests from the powertrain to the EV battery change very quickly with time, and the currents and voltages are processed as discrete-time electrical signals. Most multi-physics PDE solvers (e.g., COMSOL) are designed to model continuous-time processes (e.g., constant current) and are slow for high-frequency simulation cases. Another challenge is the high computational costs to scale-up the P2D model from celllevel to pack-level. A single-cell P2D model normally includes more than 100 discretized state variables even using a very coarse mesh. To meet the high power and energy demands, typical EV battery packs are made of hundreds to thousands of cells with extra connecting auxiliaries. Most of these EV packs employ a "mPnS" connection pattern to achieve effective battery requirements and management. For example, Tesla's newest production version, Model 3, uses 2,976 Panasonic cylindrical cells (the 21700 type) to form its 31P96S battery pack; that is, 31 cells are connected into a parallel module, and 96 such modules are connected in series to make the pack. In addition to cells, there are other components such as busbars and circuit boards in the pack, and the impedances of these components can generate current and voltage imbalance between the cells. If each of these battery cells is repre-sented by a P2D model, the memory and processing-power requests would be extremely huge for running a pack-level simulation.
Toward the first challenge, several control-oriented reformulations for the P2D model have been proposed, including the state-variablemodeling (SVM) approach by Smith et al. 5 and the discrete-time realization algorithm (DRA) by Leek et al. 6 and Lee et al. 7 These methods require linearization of the electrochemical equations and thus lose accuracy. Also, the gain and pole values used by the SVM or DRA method need to be pre-determined through implicit approaches and there are parameterization issues to implement these algorithms. Due to the above-mentioned limitations, we recently developed nonlinearstate-variable-modeling (NSVM) method, 8,9 which maintains all the nonlinear features while solving the P2D model on discrete highfrequency time space. The accuracy and time-efficiency of NSVM have been confirmed at cell-level by the work in References 8 and 9. In addition, the NSVM provides a broader operation window and much more freedom for users to customize parameters.
Reduced-order-modeling (ROM) is a promising solution to the second challenge. In recent years, many different ROM models have been reported toward the simplification of physics-based Li-ion cell model. These ROM approaches can be generally classified into two types. The first type of ROM maintains the basic physical features of P2D model and uses simplification for some of the governing equations; for example, Subramanian et al. 10 applied polynomial approximation to the electrolyte diffusion, and Kim et al. 11 derived the G-H function approach to simplify charge conservation in the electrode domain. These models still contain more than 20 differential algebraic equations (DAEs) and may not achieve sufficient cost-efficiency for large-scale simulations. Another type of ROM completely discards the physical approach, and regards the cell as a simple RC equivalent circuit 12 or even an internal resistance component. 13,14 The parameters in these ROMs need to be validated from the results of a rigorous model or directly from test data. For example, the RC components in an equivalent circuit can be evaluated by fitting the high-powerpulse-cycle (HPPC) profiles through various state-of-charge (SOC) and temperature ranges. These ROM approaches show poor robustness when working outside the SOC and temperature windows where the parameters can be effectively validated. For example, the RC models always fail for continuous discharges at SOC<10%, because the nonlinear mass transfer limitations in these operation conditions cannot be captured during the HPPC validation procedures.
The purpose of this work is to introduce an adaptive ROM approach developed based on our NSVM algorithm. We found that while solving a P2D cell model using NSVM with certain instantaneous current input, some operators can be further utilized to approximate the voltage responses when different current perturbations are added to the A3603 system. For a mPnS battery configuration, the average current is calculated by dividing the total pack current by the parallel number m, and the deviation between the individual cell currents and the average current can be regarded as a current perturbation. Therefore, by feeding the average current to a stand-alone P2D model at each time instance, a ROM can be derived from the NSVM procedure and then can be applied to calculate the voltage responses of each individual cell. As this workflow is implemented with real-time input signals, the ROM parameters can also be optimized on a real-time basis to obtain improved accuracy. The details for this mathematical method are shown in the following sections.

Mathematical Model
Newman's P2D model equations are given in Appendix A, including the illustrative schemes ( Figure A1), governing equations (Table A1), and expressions for some physical variables (Equations A1 through A3).
General form of reduced-order model (ROM).-In the P2D model, the cell voltage is defined as equal to the solid phase reference potential for cathode: where V P2D is the cell voltage estimated by P2D model at a specific current I , and φ ref S is the solid phase potential defined at the boundary of cathode current collector (x = l n +l s +l p ). In this work, the variable φ ref S is split into three terms: the equilibrium voltage U eq , the electrical over voltage V e , and the diffusional over voltage V d : [2] and therefore V P2D is given by: [3] Let I stand for a small perturbation in the current, and the cell voltage at current I = I + I will be estimated using a reduced order model (ROM) as follow: [4] where V ROM is the voltage estimated by the ROM, U eq , V e , and V d are respectively the equilibrium voltage, electrical over voltage, and diffusional over voltage estimated at the current I = I + I . I and I take positive values for charge and negative values for discharge. The derivation of these variables will be given later.
Equilibrium voltage.-According to Reference 9, the dimensionless concentration of Li ions in active material particle, θ avg , is governed by the following equation: [5] where c s,max is the maximum Li ion concentration in the active material, F is the Faraday's constant, j s is the local electrochemical current density, R s is the radius of active material particle, and t is time. The electrode average concentration of Li ions through anode and cathode domains, θ n and θ p , are defined by: where l n , l s , and l p are respectively the thickness of anode, separator, and cathode. The overall charge conservation in the anode and cathode domains are respectively given by: where a s,n and a s,p are respectively the specific surface area of anode and cathode, and A c is the coating area of the cathode. By substituting Equation 5 into 6 and applying Equation 7, the governing equations for θ n and θ p can be derived as follows:  [8] where the subscripts "n" and "p" respectively denote anode and cathode. When a cell is fully equilibrated, there are no concentration gradients in the solid and electrolyte phases, and the dimensionless Li ion concentration in anode and cathode particles are uniform and respectively equal to θ n and θ p . Therefore, the equilibrium cell voltage is calculated by: where U n and U p are half-cell open circuit potentials for anode and cathode whose profiles are presented in Appendix A. The following constraint for θ n and θ p can be derived from Equation 8: Equation 10 can also be expressed as follows by integrating with time: a s,n R s,n c s,n,max θ n + a s,p R s,p c s,p,max θ p = C 0 [11] where C 0 is a constant which can be evaluated at initial time t = 0 : therefore θ p is a linear function of θ n : θ p = C 0 − a s,n R s,n c s,n,max a s,p R s,p c s,p,max θ n [13] According to Equations 9 and 13, the equilibrium cell voltage can be expressed as a single variable function of θ n : For the ROM, the following equations are used to evaluate U eq at current I = I + I ,: [15] U eq = F θ n [16] where θ n is the anode average dimensionless Li ion concentration for current I = I + I . Next, we need to express Equations 15 and 16 in state-space form; from discretized time instance variable θ n can be evaluated by the forward time-integral method: t A c a s,n R s,n Fc s,n,max [17] where "[k]" and "[k − 1]" in variables stand for values at time instances t[k] and t[k − 1] . Therefore U eq is computed at each discretized time instance: Journal of The Electrochemical Society, 164 (14) A3602-A3613 (2017) Electrical over voltage.-According to Equation 58 of Reference 9, through a linearized Butler-Volmer equation, the reference solid phase potential can be approximated bỹ whereφ ref S is the approximated reference solid phase potential, and i B = − I Ac is the inward current density at the cathode/current collector boundary(x = l n + l s + l p ). The expressions for u B and z B are given in Reference 9. The Butler-Volmer equation turns from hyperbolic to linear as the current density goes close to zero, which means for currents close to 0. As discussed above, u B stands for φ ref S evaluated at zero current. According to Reference 9, the value of u B only depends on the concentrations of solid and electrolyte phases; therefore at time instance k, the value for u B can be expressed as function of discretized solid/electrolyte concentrations: [20] where x[k] is the vector for the discretized concentration variables. In our NSVM, the concentration states are updated by the forward time integral approach: where , , and B are matrix operators, and d where R e [k] is the instantaneous kinetic resistance of the cell given by: According to Equations 22 and 23, R e [k] can also be expressed as follows: [26] To calculate R e [k], we start from Equations 79 and 81 in Reference 9 which can be written as follows when the Newton loop converges at time instance t[k] : where Jac is the Jacobian matrix, j s is the vector for electrochemical current densities and φ ref L is the reference potential of the electrolyte. In Equation 27, variables dj s , dφ ref L , dφ ref S , and di B are all close to zero. The Jacobian matrix Jac can be expressed blockwise as follow: where J 2 is the right-most column vector of matrix Jac and J 1 is the remaining part of Jac . Equation 27 can be rewritten as follows by substituting Jac with [ J 1 J 2 ] and i B with − I Ac : Left-multiply Equation 29 with an operator [ 0 0 · · · 1 ](J τ and according to Equations 26 and 30, the value for R e can be calculated by: In state-space form, Equation 24 is written as: Diffusional over voltage.-According to the definition in Equation 59 of Reference 9, u B is not necessarily equal to U eq as concentration gradients might exist unless the cell is fully stabilized. Therefore the diffusional over voltage V d is defined as the difference between u B and U eq : According to the above discussion, the diffusional over voltage V d is determined by the mass transfer processes in solid and electrolyte phases, therefore we assume that V d can be approximately described by the following dynamic model: where τ d and k d are, respectively, the time constant and gain for the model. Rewrite Equation 34 into the following state space format: where the parameters b 1 [k] and b 2 [k] are respectively expressed by: can be regarded as parameters in a linear model: [37] where the dependent variable vector y[k] is defined as follows: The independent variable matrix X[k] is defined as follows: and the parameter vector b[k] is defined as follows: The values for b[k] can be optimized by the linear least square method: In this work, we use the recursive least square (RLS) procedure 15  can be estimated at current I = I + I by the following state space equation: [43] and V d is calculated by: Battery module implementation.-In this work, the ROM algorithm was implemented on a battery module as illustrated in Figure 1. This battery module includes 24 cells connected in parallel (24P module) by two busbars, and the positive and negative electrodes of each cell are welded to busbars through tabs. The cells in this battery module are labeled as cell #01 through cell#24, I i and V i (i = 1, 2, · · · , 24) respectively stand for current and voltage of each cell, p,i and n,i respectively stand for the positive and negative busbar potentials of each cell, T,p and T,n are respectively positive and negative busbar potentials at the external battery terminals, R Busbar is the busbar resistance between two adjacent cells, R Tab is the tab resistance of each cell, and I T is the total current passing through module. All the cells are assumed to have identical design and properties, and thermal effects are neglected. According to Figures 1b and 1c, the following equations can be derived for the busbar potentials: [45] [48] The overall current balance for the module is described as follows: and the terminal of negative busbar is chosen for potential reference For each cell, the correlation between busbar potentials and cell voltages are given by: p,i − n,i = V i + I i R Tab (i = 1, 2, · · · , 24) [51] To estimate cell voltage V i , a straightforward method is to apply full-order P2D model to each cell. The flowchart for a P2D model block is illustrated by Figure 2a  which equals to V P2D,i . An alternative approach is to replace the P2D models by the ROM, but before implementing ROM for each cell, certain parameters need to be evaluated first. As shown in Figure 2b,

Results and Discussion
According to the discussion above, the 24P battery model can be solved by either running 24 full-order P2D models (FOM method) or using 24 ROMs with a stand-alone P2D model (ROM method). Comprising only two state variable equations and a few explicit algebraic expressions, the ROM significantly reduces the complexity of a P2D model which involves several nonlinear PDEs. To verify the strength of our approach, both FOM and ROM methods were first implemented to simulate an aggressive drive cycle. The parameter values for battery cells are given in Table I and the input current profile for the drive cycle is presented in Figure 3. The peak discharge current for the drive cycle is equivalent to 2.0C (12A in average for all the cells) and the sampling frequency is 1 Hz. This simulation was performed by MATLAB/Simulink software, the computation time for FOM is 104.8 sec and for ROM is 1.95 sec. Figure 4 shows the simulated voltage results by FOM and ROM, where the voltage profiles of battery terminals and three cells (#01, #12, #18) are compared. Zoom-up plots are presented in each case to catch the details of voltage comparisons. The root-square-mean-deviations (RSMD) between FOM and ROM are calculated for each plot in Figure 4. The simulated FOM and ROM current profiles for cell#01, #12, and #18 are shown in Figure 5 with zoom-up plots and calculated RSMD. The simulation results in Figure  4 and Figure 5 suggest that the ROM achieves about 97% accuracy. In addition to the FOM vs ROM comparisons, Figure 5 also shows a significant current distribution through this battery module: the maximum discharge current is −8.05 A for cell #01 while −20.0 A for cell #12; these values deviate from the average current by 70%, so the ROM has good confidence to predict the electrical imbalance in a battery. Besides the high-frequency drive cycle, FOM and ROM are also used to simulate constant current (CC) discharges with 0.5C, 1C, and 2C. The profiles for terminal voltage vs depth-of-discharge (DOD) for different CC cases are displayed in Figure 6, and the computation time and RSMD values for each case are listed in Table II. The above simulation results show that ROM is 50∼70 times faster than FOM. However, note that in our simulation, the FOM has already been sped up by the NSVM method (50∼100 times faster than COMSOL solver according to References 8 and 9), therefore the time-efficiency of ROM is significant. Although its accuracy decreases with current magnitude, the ROM method still exhibits outstanding performance considering the low computational cost.
In Figures 7a through 7e Figure 7c, V d [k] drops significantly when DOD exceeds 90%, which is caused by the depletion of intercalation sites at the surface of  cathode. The two diffusional parameters b 1 [k] and b 2 [k] take smooth profiles during all four cases except for a short unstable initial period. The reason is that b 1 [k] and b 2 [k] are estimated recursively from initial guesses, and it may take several iterations to settle these parameter values. As shown in Figure 7, the profiles for a specific parameter vary dramatically in different operation cases. These results suggest that our ROM is an adaptive algorithm, which means the parameter values can be optimized according to the real-time current inputs; if the input signal changes, the algorithm is able to re-evaluate parameters to achieve higher accuracy. A reduced-order-modeling (ROM) method was developed for the P2D Li-ion cell model based on the specific solution procedure of the NSVM algorithm. By sharing the real-time parameters estimated by a stand-alone P2D model, the ROM, which shows great accuracy and adaptiveness for various operation inputs, can be distributed into different cells in a battery pack. Compared with the full-order-model (FOM), which employs distributed P2D models and is solved by the NSVM, the ROM improves the computation efficiency by 50∼70 times. The ROM only contains two state equations (one for θ n [k] and the other one forṼ d [k]), and its numerical complexity is just comparable to a 1-RC equivalent circuit model. As a stand-alone P2D model is a pre-request for running the ROM, only when used for simulating large battery systems can the advantages of ROM be realized.
In the future, this ROM approach can potentially be used to couple with heat transfer equations in the Multi-Scale-Multi-Dimensional (MSMD) framework, 16 and model the distributed electrical and thermal behavior inside the cell volume. The feasibility for this proposed transformation is based on the model structure in Reference 17, where the planar electrode pair is treated as a 2D circuit network, and the ROM would replace the electrochemical model to compute the current flow at each grid.

Appendix A
Illustrations of the porous electrode domains for the P2D model are shown in Figure  A1 and the governing equations are summarized in Table A1. The open circuit potential expressions are given below:  where U n and U p are respectively the open circuit potentials of the anode and cathode and θ is the dimensionless Li concentration at the surface of solid particles. The expression for the electrical conductivity of the electrolyte, κ L,bulk , is as follows: where T is temperature and c L is the electrolyte concentration.

Appendix B
For time instance k, Equation 41 can be converted into the following form; and a similar equation can be derived for time instance k − 1 :

Physics Equations
Electrolyte diffusion Defined only for 0 ≤ x ≤ l n and l n + l s ≤ x ≤ l n + l s + l p ) Solid phase diffusion c S,max = 0 (Defined only for 0 ≤ x ≤ l n and l n + l s ≤ x ≤ l n + l s + l p ) Butler-Volmer equation RT η) − exp(− 0.5F RT η)] (Defined only for 0 ≤ x ≤ l n and l n + l s ≤ x ≤ l n + l s + l p ) The following equations can be derived according to Equation B3: [B4] where χ[ j] is the j th row vector of the matrix X . Substitute Equation B4 into B2 to yield and the following equation can be derived by subtracting Equation B5 from Equation B1: Left multiply Equation B6 by (X τ [k]X[k]) −1 and rearrange terms to get the following recursive expression for b[k]: Instead of recording all history, the value for X τ [k]X[k] can be updated at each time instance using: In our specific case, X τ [k]X[k] is a 2 × 2 symmetric matrix which can be expressed as follows: and the elements f [k], g [k], and h[k] are estimated by: where χ 1 [k] and χ 2 [k] are respectively the 1 st and 2 nd element of vector χ[k] . Therefore, the covariance matrix is estimated as follows: [B11] Substitute Equations B11 into B7, and the individual parameters b 1 [k] and b 2 [k] can be calculated using the following equations: where the expression for the gain value K RLS is given by: