A Padé Approximate Model of Lithium Ion Batteries

We present a reduced order model for a lithium ion battery in which Pade approximants are used to simplify complex, transcendental, transfer functions associated with the linearized, pseudo 2-dimensional (P2D) electrochemical model of the battery. The resulting transfer functions take the form of simple, rational polynomial functions, which can be used to compute any variable at any location within a one-dimensional representation of the battery domain. Corrections for nonlinear behavior are also incorporated into the reduced model. The results obtained using the reduced model compare favorably to those from the full (nonlinear) P2D model and the computational time required to produce these results is significantly reduced. Importantly, the form of the reduced model means that variables can be evaluated at specific discrete locations within the cell domain, without the need to compute all values of the variable at all discrete locations, as is the case with the spatial discretization methods most commonly used to implement the P2D model. We show that this results in further significant time savings and enhances the suitability of the model for variety of applications.

Lithium ion batteries are now commonly used for large renewable energy storage systems, as either as part of the power grid or in stand-alone installations. 1 A Battery Management System (BMS) is essential in order to utilize battery storage optimally, achieve safe charging and discharging of the battery system and prolong the battery lifetime. 2 The optimal control algorithms within a BMS require a battery model that can describe the voltage response and internal electrochemical parameters of the battery system during the charge/discharge process. An accurate battery model that can represent the complex and nonlinear internal characteristics of a battery is desirable in order to achieve optimal utilization, maintain safety and prolong battery life. Thus, a highly accurate battery model is an advantage in advanced BMSs.
Depending on the application, the energy capacity of a battery storage system can vary from kWh to MWh. 1 Given that a single cell has typically only 7.5Wh energy capacity, 3 large capacity systems are usually built using many battery packs, which in turn consist of large numbers of lithium ion cells connected in series and parallel to achieve the required voltage and capacity level. 4 As a result, not only is a highly accurate battery model desirable but the model itself must also be able to be evaluated quickly, using a highly efficient computational algorithm.
Current BMSs rely significantly on equivalent circuit models due to their simplicity and the inherently low computational burden required to implement them. 5 However, equivalent circuit models require extensive battery testing in order to validate them in a particular operational range in which the battery circuit parameters have been identified. 6 In addition, the models give only limited information on the electrochemical characteristics and degradation response of a battery. 6 As an alternative to equivalent circuit models, the pseudo 2-dimensional electrochemical (P2D) model, first developed by Doyle,Fuller,and Newman,7,8 does provide detailed insights on these characteristics. 6 However, the P2D model consists of partial differential equations (PDEs). To solve these, spatial discretization methods are applied to yield a system of differential algebraic equations (DAEs) in the time domain. 9 These methods are usually implemented on a desktop computer, often using commercial software, and are therefore not suitable for embedded BMSs, which require much lower computational workloads. For example, as reported by Lee, Chemistruck and Plett, 5 the P2D model implemented in COM-SOL Multiphysics 10 requires about 13 min on a desktop computer, to produce 25 min of simulated voltage response to the Urban Dy-namometer Driving Scheduler (UDDS) 11 current profile. Such computational requirements render a direct implementation of the P2D model infeasible for use in a BMS.
Padé approximations provide a way of greatly simplifying complicated transcendental transfer functions in order to produce simpler rational polynomial transfer functions that can be used in embedded applications. 12 Previously, Padé approximations have been applied to the transcendental transfer functions that model lithium ion diffusion in a single, solid, spherical particle in the so called Single Particle Model (SPM) 13,14 for lithium ion batteries. Marcicki et al., 15 Zhang et al. 16 and Yuan et al. 17,18 used Padé approximations for an extended version of the SPM, which incorporates concentration change within the electrolyte solution while maintaining the uniform current density assumption of the SPM. These authors report that their Padé approximate model is computationally more than 50 times faster than the P2D model. 18 The batteries used in storage systems for grid applications are designed to have a high energy density; as opposed to cells that are specifically designed to deliver high power (like those used in electric vehicle applications 9 ). Such high energy density cells are generally made with thicker electrodes, 9 which results in non-uniform, lithium ion concentration and current density profiles across the electrode. 9 However, in both the SPM and extended SPM the current density is assumed to be constant, spatially. Thus, these models can yield inaccurate predictions when applied to the thick electrodes of high energy density cells. 9 In this work, we propose an improvement to the previous SPM and extended SPM Padé approximation models discussed above. To do this we consider the transcendental transfer functions of the linearized P2D model and the corresponding corrections for nonlinear behavior given by Lee, Chemistruck and Plett. 5 These transfer functions describe the electrochemical dynamics of the lithium ion battery and the corrections for nonlinear behavior are used to improve the accuracy of the linearized P2D model. 5 We apply Padé approximations to these complicated transcendental transfer functions in order to develop simpler, rational polynomial transfer functions, which form a reduced model of cell operation that is amenable to rapid computation and therefore embedded BMS applications. We believe that this is the first attempt at using the Padé approximation method to simplify all transcendental transfer functions of the linearized P2D model. Our reduced model output includes the spatial and temporal variation of all of the variables of the P2D model. In addition, it can be used to determine any variable at any specific spatial location without the need to compute variables at all of the discretized spatial locations; as is the case when using spatial discretization methods to solve these models. Importantly, the physical meaning of the approximation model can be preserved since the coefficients of these rational polynomial transfer functions directly link to the physical parameters of the cell such as, SEI layer resistance, particle size of the active material, lithium ion diffusivity and solid and solution phase volume fractions. If the operating conditions of the cell, and hence the values of the physical parameters, change (for example, due to cell aging) these coefficients are readily updated as they themselves have computationally simple algebraic forms. Our results show that the reduced model delivers high accuracy for low computational overhead when compared to the P2D model.

Lithium Ion Battery Electrochemical Model
The P2D model is applied to 4 domain elements of a lithium ion battery, namely, the negative electrode, the separator, the electrolyte and the positive electrode. 9 During discharge, lithium ions in the negative electrode diffuse from the interior to the surface of the spherical solid particles that constitute the porous electrode, where they undergo an electrochemical reaction and are transferred into the electrolyte phase. Thereafter, the ions travel through the electrolyte solution to the positive electrode where they again react at the surfaces of the solid particles (constituting that porous electrode) and are intercalated into these particles. These processes are shown schematically in Figure 1 along with the definition of our modeling domain in which the negative electrode (neg) is defined from x = 0 to x = L neg ; the separator (sep) is defined from x = L neg to x = L neg + L sep and the positive electrode (pos) is defined from x = L neg + L sep to x = L tot . The reverse of the depicted discharge process occurs in the charge process, whereby lithium ions flow from the solid phase of the positive electrode to the solid phase of the negative electrode. As the separator forms an electrically insulated barrier between the electrodes, the flow of electrons associated with the charge and discharge processes occurs via an external circuit connecting the two compartments, doing useful work as they do so. 19 The P2D model that describes the above processes consists of five state variables including the electric potential φ s (x, t) in the solid electrode, the electric potential φ e (x, t) in the electrolyte, the lithium concentration of the active material c s (x, r, t) of the positive and negative electrodes, the lithium concentration c e (x, t) in the electrolyte, and the molar fluxes j(x, t) of the charge that flows between the active material in each electrode and the electrolyte. The governing equations of the model are given as follows 9 with the definitions of the symbols shown in the equations being given in the List of Symbols.
The potential in the solid phase of each (negative and positive) electrode is derived from the principle of conservation of charge and is given by (for i ∈ {neg, pos}), with the boundary conditions The electrolyte potential is defined along the whole cell domain and is given by (for i ∈ {neg, sep, pos} and j sep (x, t) = 0), where, with the boundary conditions Lithium ions are transported in the solid particles of each electrode by diffusion and their concentration is given by (for i ∈ {neg, pos}), [4] where the boundary conditions and initial condition are given by [5] The boundary conditions and initial condition are respectively. The transfer of charge at the solid/electrolyte interfaces in each electrode is given by (for i ∈ {neg, pos}), where, c i s,e is the surface concentration of lithium in a spherical electrode particle.
The overpotential of the intercalation reactions in each electrode is given by (for i ∈ {neg, pos}), The potential difference between the positive and the negative current collectors yields the cell voltage, namely, [8] Equations 1-8 constitute a system of coupled PDEs and algebraic equations that represent the P2D model for a lithium ion battery. Under galvanostatic conditions, the input of the model is the applied current I app (t) and the output of the model is the cell voltage V cell (t). In the separator (L neg ≤ x ≤ L neg + L sep ), the model involves two coupled PDEs which are the electrolyte concentration given by Equation 5 and the electrolyte potential given by Equation 2. In the electrodes (0 ≤ x ≤ L neg ,L neg + L sep ≤ x ≤ L tot ), the model involves the coupled PDEs for electrolyte concentration, given by Equation 5, electrolyte potential, given by Equation 2 and solid phase potential, given by Equation 1 on the macro (cell) scale (x) and Equation 4 on the micro (particle) scale (r) (which exists at every x in the electrode domains). 20 The nonlinear nature of the P2D model means that it must be solved numerically. As exemplified by the work of Farrell and coworkers, 21,22 this is generally achieved using a desktop computer running finite volume or finite element software. As noted earlier, this complexity prevents the implementation of the P2D model in embedded real-time applications such as BMSs.

Analytic Laplace Domain Transfer Functions of Linearized P2D Model
Truncated Taylor series expansions about a set point can be applied to each of the variables of the P2D model to obtain a linear form of the model. The details of this are given by Lee, Chemistruck and Plett. 5 The set-point, p * , used by these authors is given by the initial conditions of the cell, defined as p * = { φ s,e = U OC P (c s,0 ), c s,e = c s,0 , c e = c e,0 , j = 0}, where φ s,e = φ s − φ e . These authors also define a set of transfer functions for the cell, shown in their respective domains in Figure 1, where the dimensionless distance z is related to the spatial variable x via the transfor- . Subsequently, the "∼" notation denotes the variables φ s,e = φ s,e − U OC P (c s,0 ), andc s,e = c s,e − c s,0 . The transcendental transfer functions are formed by first linearizing the governing equations of the P2D model around the set point p * and then taking Laplace transforms of these linearized equations. The details of these transfer functions are as follows.
For the phase potential difference φ neg s,e in the negative electrode: φ neg s,e (z, s) [9] where: [11] and s is the frequency parameter in the Laplace domain.
For the molar flux in the negative electrode: . [12] For the over potential in the negative electrode: η neg (z, s) For the solid phase Li + concentration on the surface of the particles in the negative electrode: , [14] For the potential in solid phase of the negative electrode: . [15] The transfer functions for the electrolyte potential in the positive electrodeφ For the Li + concentration in the electrolyte: [16] where,C * e,k (s) , [17] Here λ k is the k th eigenvalue, and ψ(x; λ k ) the corresponding eigenfunction, forC e (x, s)/I app (s) and K is the order of the expansion. 5 The transfer functions, J neg * k (s)/I app (s)and J pos * k (s)/I app (s) as well as those for the electrolyte potential in each domain are given in the Appendix.
The nonlinear corrections for the internal variables are computed as follows: 5 [ ϕ e (x, . [23] Finally, the cell voltage is calculated from: [24] The advantage of these transfer functions is that they model the P2D variables at any spatial location in the cell domain based only on the applied currentI app (s). The transfer functions given in Equations 9-17 and Equations A1-A5 (in the Appendix) are complicated and computationally expensive to calculate which, like the P2D model from which they are derived, makes them not amenable to real time application. The Discrete Time Realization algorithm (DRA) proposed by Lee, Chemistruck and Plett 23 can reduce the transcendental transfer functions into an optimized, reduced order, discrete-time, state-space model. In the DRA, Lee et al. first compute the discrete-time pulse response of the transcendental, continuous-time, transfer functions (Equations 9-17 and Equations A1-A5). They then apply the Ho-Kalman algorithm to produce the state-space realization from this discrete-time impulse response. These authors show that the resulting simple, state-space model produced from the DRA can reduce the computation time requirements by more than 5000 times, compared to the P2D model. However, the DRA requires about 11 min on a desktop computer to produce the reduced state-space model from the transcendental transfer functions for a given set of parameters (e.g. R film , r eff , D s , D e , ε e and ε s ). There have been two major improvements to the classical DRA in terms of reducing the computation overhead, the first by Gopalakrishnan et al. 24 who implemented an improved singular value decomposition scheme on a virtual Hankel matrix to speed up the DRA. This resulted in a reduced order model that was approximately 100 times faster than that using the classical DRA. The second by Rodríguez et al. 25 who applied variation of parameters to speed up the computation of the electrolyte concentration transfer function. In this case, the modified transfer function can be computed 3800 times faster than the previous electrolyte concentration transfer function given by Lee, Chemistruck and Plett. 5 The speed up of an overall reduced order model using this modified transfer function is not considered by Rodríguez et al., 25 however. Regardless of these improvements, when the values of the parameters change, for example, due to cell aging, the DRA needs to be re-run to regenerate another state-space model. This results in an increased run time for DRA based reduced order models that may prevent their practical use in real-time applications where rapid parameter updating is required.

Padé Approximation Method
In order to simplify the above transfer functions, we apply Padé approximations. The aim of this approach is to reduce the transcendental nature of the above transfer functions to ones that contain only rational functions of simple polynomials in the Laplace domain. 12 We note that an (M, N) order Padé approximation of a given transfer function,G i (s), centered at s 0 has the form, 12 where the superscript i associates the quantity with the corresponding transfer function, namely, phase potential, molar flux, overpotential, lithium concentration and solid phase potential as given by where the coefficients γ i k can be determined from the power series expansion of the transfer function G i (s), namely, Equation 26 generates a polynomial of order N(N+M+2) in s. 12 Its right hand side equals zero for all s, therefore, the coefficients of s equal zero. 12 Consequently, we can create a set of M+N+1 linear equations given in Equation 28, which includes N equations dependent on b i n and M+1 equations dependent on a i m and b i n , namely, . [28] where the zeroth order term, b i 0 , is assumed to be equal to 1 to normalize the solutions. 12 Solving the equations in 28, we can obtain the M+1 coefficients, a i m , (m = 0, ..., M), and N coefficients, b i n , (n = 1, ..., N ). 13,26 There are several readily available routines that automate the above methodology. In this work we used the PadeApproximant [G i (s), s 0 , {M, N }] function in Wolfram Mathematica 27 to determine the Padé approximations to our transfer functions.
As a specific example we note that the transfer function of the molar flux in the negative electrode,J neg (z, s)/I app (s), given by Equation     . Current input and corresponding frequency content of a battery storage system that is being used to smooth the power generated from a wind farm connected to a grid. Current data sourced from Ref. 29. can be approximated to first order by the Padé approximation, where, Similarly, the remaining transcendental transfer functions, introduced above, can be approximated by their rational polynomial, Padé approximants, given in Figure 2. In this way we can develop a system of transfer functions that represent the variables of the P2D model for a lithium ion battery and which have the advantage of being easy to implement and compute in a real-time application. In addition, the physical meaning of the approximation model can be preserved since the coefficients of these rational polynomial transfer functions directly link to the physical parameters of the cell such as R film , r eff , D s , D e , ε e and ε s . If the operating conditions of the cell and hence the values of the physical parameters change, these coefficients are readily updated as they themselves have computationally simple algebraic forms as shown in Equations 30-33. This is in distinct contrast to the DRA approach discussed above. We note that the outputs of the Padé approximants are used within their corresponding nonlinear correction term, given by Equations 18-23, in order to improve the accuracy of the approximations to the purely linear P2D model. 5 Further improvement could be achieved by incorporating this approximation model with an observer, 28 however, we leave this as future work. Figure 3 compares the frequency response of the Padé approximant with their transcendental transfer function counterparts for the four internal electrochemical variables within the negative electrode. We note that the frequency response of each set of transfer functions  (where a set consists of a number of curves corresponding to the range of z values between z = 0 and z = 1) depicted in Figure 3 is for a particular set-point p * . We consider a range of p * values for each set, resulting in a family of sets, with each member of this family corresponding to a different p * . From this family of sets, the actual p * value that was used to produce a given set shown in Figure 3 was the value that maximized the root mean squared error (when compared to the linearized P2D model). This ensures that the frequency responses of the Padé approximations shown in Figure 3 are the worst case scenarios (i.e. the ones with maximum error when compared to the linearized P2D model).
We see that the error between the Padé approximants and the transcendental transfer functions increases as the input frequency is increased. However, we note that in grid applications low frequencies dominate the battery input current profile as shown in Figure 4. Figure  4a depicts the current profile, reported by Li et al., 29 of a battery storage system that is being used to smooth the power generated from a wind farm connected to a grid. Figure 4b shows the fast Fourier transform of the current signal in Figure 4a. Figure 4c plots the Power Spectrum Density (PSD) and identifies 99% occupied bandwidth, 30 which is calculated by determining where the integrated power crosses 0.5% and 99.5% of the total power in the spectrum. 30 As shown in the figure, 99% of the occupied bandwidth of the power signal have frequencies less than 8.225 mHz. Comparing the Padé approximants with the linearized P2D model outputs in the mHz frequency in Figure 3 we observe that the two responses show an acceptable match, suggesting that it is reasonable to apply the low order Padé approximants for grid scale applications.

Results and Discussion
Comparisons of Padé approximants of various cell variables with those obtained from the full P2D model under various applied current conditions are given in Figures 5, 6 and 7. The P2D model has been implemented using the open source code supplied by Torchio et al. 31 The cell parameters used for these simulations were obtained from Ref. 5 and are listed in Table I. Figure 5 compares the electrolyte concentration c e (L tot , t), c e (0, t) and the surface concentration in the solid phase c s,e (L tot , t), c s,e (0, t) and cell voltage for a pulse current profile. Figure 6 shows a similar comparison when the Urban Dynamometer Driving Schedule (UDDS) current profile is applied and Figure 7 shows the comparison when the input current is that of the wind farm, grid application from Ref. 29, as depicted in Figure 4a.
In all of these cases we observe that the Padé approximations compare very well to the P2D model outputs. The computation times (in seconds) for each of these applied current profiles is shown in Table II. We note that all simulations were implemented in MATLAB R2016b and run on the same desktop computer with an Intel Core i7-6700 CPU running at 3.4 GHz and 16GB RAM. In addition, the Padé approximants are determined at the same discrete spatial locations as given by the FVM grid used to solve the P2D model. We observe that the Padé approximations require approximately 9.6 s to simulate the pulse current profile in comparison to the P2D model, which takes approximately 15 s. The time savings of the Padé approximations are more evident in the cases of the UDDS and grid current profiles, with the Padé approximation model solving in approximately 8 and 42.5 s, respectively and the P2D model taking approximately 539 and 22,494 s, respectively. It is this feature of the reduced model coupled with the fact that the form of the Padé approximations, once obtained, only require relatively simple rational functions to be evaluated, which leads to the significant computational savings observed in Table II. Like the P2D model, the Padé approximation model can be used to predict how important cell variables will evolve both spatially and temporally. However, another significant advantage of the Padé approximation approach is that variables can be evaluated at discrete locations within the domain, without the need to compute all values of the variable at all discrete locations, as is the case with the FVM implementation of the P2D model. This can result in further  j(x, t), c s,e (x, t), ϕ e (x, t), c e (x, t) and η(x, t) at only two locations, x = L tot and x = 0, are required to compute the cell voltage. If we only wish to calculate the cell voltage then, as shown in Table III, the computation time of the Padé approximation model reduces to approximately 4.7 s for the pulse current profile, and approximately 4.3 and 29.1 s for the UDDS and grid-scale current profiles, respectively.
We now consider the application of our reduced model to battery charging. The traditional constant current-constant voltage (CC-CV) charging method is widely used to charge lithium batteries. 32 As the name suggests, this method applies constant current (CC) to charge the battery to a preset voltage at which time the charging mode switches to a constant voltage (CV) regime until the current reduces to a preset value. The CV regime is often lengthy resulting in a long charging time and in addition, as the battery ages, the preset voltage limit for the CC regime often results in the cell being overcharged. 9 When the cell is being overcharged, it can be damaged and the risk of explosion is potentially increased. 9 A faster and safer method of charging has been reported by Chaturvedi et al., 9 and is based on observing the value of the potential difference φ s − φ e at the negative electrode/separator boundary (x = L neg in Figure 1). We note that at this point, φ s − φ e , will be a minimum within the negative electrode during charging. 33 Thus, when φ s − φ e at x = L neg is negligible (point B shown on Figure  8a), we should stop charging to avoid overcharge and this is the reason for monitoring this point in the approach presented by Chaturvedi et al. 9 Figure 8a depicts the cell voltage and Figure 8b depicts the over potential, φ s − φ e , as given by the Padé and P2D models for a single charge step at a rate of 1C. In the voltage limit method, we monitor the cell voltage and switch from constant current charging to a constant voltage mode when the cell voltage reaches 4.2V. From Figure 8a we can see that at the point of switching (point A, 4.2V) this corresponds   to a state of charge (SOC) of 93.2%. This SOC will continue to increase under the constant voltage regime, but will do so very slowly. From Figure 8b, however, we can see that at the point of switching the potential difference φ s − φ e at x = L neg is greater than 0V. Therefore, under the method proposed by Chaturvedi et al. we would continue to charge at a constant current until φ s − φ e = 0, yielding a 6.8% extra SOC. The additional time required for gaining 6.8% extra SOC in this constant current regime is 243s, which is much faster than the time required to gain the same amount of charge if we had switched to a constant voltage regime (which would take on the order of hours). The pivotal point here is that the values of φ s − φ e at x = L neg predicted by our Padé approximation model compare very favorably to those predicted by the P2D model but have the advantage of being easily computed. An accurate and fast prediction of the φ s − φ e values at x = L neg , which cannot be easily measured empirically, means that our Padé approximation method could be useful for control purposes in fast charging. Finally, we note that the accuracy of our proposed Padé approximation model can be further improved by increasing the order of the rational functions used in the approximations. However, higher order models result in higher computational requirements. In order to investigate this tradeoff, we increase the approximation order as shown in Figure 9. Figure 10 depicts the voltage error (when compared with the P2D model predictions) for the low order approximations given in Figure 2, and the higher order approximations given in Figure 9 for our three previous applied current profiles (i.e. the pulse current profile, UDDS, grid-scale application). Table IV shows the root mean square value of the voltage error over the length of each current profile. From these figures we observe that the accuracy of the higher order approximation model increases by some 5.4% and 7.5% for the pulse current profile and the UDDS profile, respectively. In the pulse current profiles and the grid-scale profile, there are negligibly small differences between the low and higher order approximations models. However, as can also be seen in Table III, the computational time requirement of the higher order approximation model increases by approximately 20% for the pulse and UDDS current profiles and 7.7% for the grid current profile. We see that the higher order approximation model does not result in significant improvements for any of these applied current profiles and therefore, increasing the approximation order is unnecessary.

Conclusions
We presented a reduced order model for a lithium ion battery in which Padé approximants were used to simplify complex, transcendental, transfer functions associated with the linearized, pseudo 2-dimensional (P2D) electrochemical model of the battery. The resulting transfer functions take the form of simple, rational polynomial functions, which can be used to compute any variable at any location  within a one-dimensional representation of the battery domain. Corrections for nonlinear behavior are also incorporated into the reduced model. The proposed model substantially reduces the complexity of the P2D model whilst maintaining the functionality of predicting all of the variables of that model. The form of the Padé approximations, once resolved, only require relatively simple rational functions to be evaluated, which leads to significant computational savings when solving the reduced model. Our results demonstrate that the reduced model predictions match very closely with those obtained from the full (nonlinear) P2D model but in a fraction of the computational time. Importantly, variables can be evaluated at specific discrete locations within the domain, without the need to compute all values of the variable at all discrete locations, as is the case with the FVM implementation of the P2D model. This can result in further significant time savings and enhance the flexibility of the model for a variety of applications.    [A5] ORCID Ngoc Tham Tran https://orcid.org/0000-0001-9994-3620 Troy Farrell https://orcid.org/0000-0002-6629-4174