**Model-Based Adaptive Joint Estimation of the State of Charge and Capacity for Lithium–Ion Batteries in Their Entire Lifespan**

#### **Zheng Chen 1,2, Jiapeng Xiao 1, Xing Shu 1, Shiquan Shen 1, Jiangwei Shen 1,\* and Yonggang Liu 3,\***


Received: 31 January 2020; Accepted: 13 March 2020; Published: 18 March 2020

**Abstract:** In this paper, a co-estimation scheme of the state of charge (SOC) and available capacity is proposed for lithium–ion batteries based on the adaptive model-based algorithm. A three-dimensional response surface (TDRS) in terms of the open circuit voltage, the SOC and the available capacity in the scope of whole lifespan, is constructed to describe the capacity attenuation, and the battery available capacity is identified by a genetic algorithm (GA), together with the parameters related to SOC. The square root cubature Kalman filter (SRCKF) is employed to estimate the SOC with the consideration of capacity degradation. The experimental results demonstrate the effectiveness and feasibility of the co-estimation scheme.

**Keywords:** state of charge; available capacity; adaptive model-based algorithm; square root cubature Kalman filter; joint estimation

#### **1. Introduction**

Nowadays, energy crises and environmental damage have become the main concerns of society, and require being tackled with high attention [1]. Transportation electrification provides a possible manner to reduce emissions and dependence upon fossil fuels. Electric vehicles (EVs) and hybrid EVs (HEVs) are promising solutions, which however, require electrical energy storage systems to completely or partially replace propelling power supplied by traditional internal combustion engines [2]. In this context, applications of lithium–ion batteries have been intensively spurred due to their numerous advantages, such as their wide environmental temperature operation capability, high energy density, long lifespan and their large charge/discharge current [3]. For lithium–ion batteries, the state of charge (SOC) and available capacity, usually provided by battery management systems (BMSs), are crucial parameters for evaluation of the electrical performance of the battery, as well as for the control of the vehicle.

Typically, estimation methods of the SOC can be divided into four categories, including the coulomb counting method, and characterization parameter-based methods such as the open circuit voltage (OCV) method, model-based methods and data-driven methods. Amongst them, the coulomb counting method [4] and OCV method [5] have been widely applied in BMSs of EVs, because of their simplicity and ease of implementation, whereas the former is prone to the production of large accumulated errors, due to interferences or uncertainties of current sensors/transducers and inaccurate initial values, and the latter is not suitable for online estimation, as it usually costs long shelving time to acquire the OCV value. With the development of computation technologies and machine learning, a variety of artificial intelligence-based, data-driven methods, such as neural networks [6] and support vector machines [7], are proposed for SOC estimation by establishing black-box models. Data-driven methods feature a strong nonlinear mapping capability with high accuracy; however, these approaches show high complexity, and require a considerable amount of training data. Alternatively, model-based methods have been widely investigated and applied for SOC estimation, thanks to the capabilities of online application, high precision and the independence of initial values. Conventional modeling manners mainly include electrochemical models and the equivalent circuit models (ECMs). Compared with complicated electrochemical models, ECM is commonly used to describe the electrical behavior of batteries, and subsequently to estimate the SOC due to its simplification and preferable precision. Yanwen Li et al. proposed a multi-model probability fusion algorithm to describe the battery's electrical characteristics, and subsequently estimate the SOC [8]. In model-based methods, the combination of the battery model and the intelligent filtering algorithm is a hotspot in SOC estimation research. The frequently used filtering algorithms include Kalman filtering (KF) [9], the H-infinity filter (HIF) [10], particle filter (PF) [11], and their various extensions. In particular, the extended KF (EKF) is widely employed to execute SOC estimation using a first-order Taylor expansion on the basis of the battery's nonlinear model [12]. Nonetheless, the second and higher order expansion is usually neglected, thus leading to slow a convergence rate, and even divergence. The unscented KF (UKF) is exploited to estimate battery SOC, based on the recursive unscented transformation to approximate the nonlinear observation without Taylor polynomial expansions [13]. The UKF shows better estimation precision and robustness than the EKF in strong, nonlinear systems [14]. On the basis of the radial–spherical cubature criterion, the cubature Kalman filter (CKF) leverages a set of volume points to approximate the mean and covariance of states with additional Gaussian noise [15]. Although CKF outperforms EKF and UKF in terms of filtering divergence and estimation error, it is susceptible to inaccurate, initial difference and disturbances, and is difficult to guarantee a symmetric and nonnegative definition of the covariance matrix all the time. HIF is applied in state estimation and model parameters identification of batteries, due to its good, anti-interference performance in high nonlinear systems [16]. PF exhibits attractive advantages in solving nonlinear, non-Gaussian distribution problems, and highlights more application potential than EKF. Thus, it has been widely developed and applied in multifarious fields, such as batteries, robotics and navigation systems [17]; however, it is limited by strong dependence upon noise and time-varying parameters of the system.

In addition, battery aging is an irreversible process with operation, where the most intuitive appearance is a decline of capacity and the increase of internal impedance [18]. In general, the attenuation process is nonlinear, complex, and even difficult to predict. To attain it, a body of algorithms have been successfully proposed and applied to achieve capacity estimation, mainly including experimental analysis methods and model-based methods. The most direct and easiest manner of evaluating the capacity is to conduct the calibration test [19]. However, it is obviously time-consuming, and only supports offline estimation. Additionally, the battery's impedance variation also highlights the capacity degradation trend [20]. However, the online electrochemical impedance measurement is not suitable for practical applications, due to its exceptional complexity of experiment. Motivated by these difficulties and constraints, incremental capacity analysis (ICA) is introduced to conduct capacity estimation by evaluating the increment of capacity in a certain charging interval [21]. Similar algorithms also include differential voltage analysis (DVA), with the help of analyzing the variation characteristics of voltage curves in predetermined charge/discharge operations [22]. Yes, they can reflect the aging mechanism of batteries, and highlight preferable accuracy; nonetheless, they are intractable to apply in practice, as chances of encountering the interval with predetermined current are seldom. Since it is time-consuming to measure and determine the battery capacity directly; model-based estimation approaches may supply an indirect manner to evaluate it. The model-based methods can leverage adaptive algorithms (such as joint estimation approaches [23] and fuzzy logic algorithms [24]) to identify the battery capacity.

These algorithms are easy to be implemented, and meanwhile demonstrate preferable accuracy; whereas, in the model-based approaches, the battery capacity is regarded as a key parameter or state variable, and is obtained by parameter identification, or estimated together with the SOC, with an established circuit model or electrochemical model. From this point of view, the model accuracy can directly affect the estimation accuracy of our battery capacity.

The above-mentioned state estimation methods are mostly developed for either SOC or capacity estimation individually, rather than for both simultaneously. SOC refers to the residual capacity rate over nominal values, while SOH represents the nominal capacity value with operations. To a certain extent, battery capacity shows the same significance as SOC, and essentially, they are tightly coupled with each other [25]. Apparently, SOC estimation based on the known and unchanged capacity exhibits certain limitations in practice. A common knowledge is that the internal parameters of batteries change with degradation. The internal resistance will increase, and the capacity decreases gradually, thus resulting in the challenges and difficulties of estimating the SOC reliably and robustly. Consequently, it is critical to update the model parameters, particularly the capacity, in a timely manner. Motivated by this, a joint estimation scheme is proposed in this study to improve the estimation accuracy of the SOC and capacity in the entire lifespan of the battery. Firstly, a second-order resistance–capacitance (RC)-based ECM is established, and the co-estimation scheme of SOC and the battery capacity is presented. In it, the square root cubature Kalman filter (SRCKF) algorithm is employed to estimate the SOC; meanwhile, the battery capacity, as one of the key model parameters, is identified by the genetic algorithm (GA), based on the constructed three-dimensional response surface (TDRS). Finally, the estimation results of the SOC and battery capacity are verified by different experimental validations over their entire lifespans. This study dedicates to the following two contributions: 1) A novel capacity estimation method based on a TDRS is proposed, and the model parameters are updated synchronously; and 2) based on the capacity and parameters revision, a co-estimation scheme is established for the SOC and capacity estimation simultaneously against different degradation statuses.

In the remainder of this study, Section 2 details the second-order RC model and the experiment test profiles. In Section 3, the co-estimation scheme of both capacity and SOC is elaborated. The validation results are exhibited and discussed in Section 4. Finally, Section 5 draws the main conclusions and looks to future works.

#### **2. The Lithium–Ion Battery Model and the Experimental Details**

#### *2.1. Battery Modeling and Analysis*

To better estimate battery states, various mathematical models have been established, including electrochemical models [26] and ECMs [27]. However, they differ greatly in accuracy, computation complexity and reliability. Considering the precision and complexity of models, a second-order RC-based ECM, as shown in Figure 1, is deployed in this work, thanks to its relatively satisfactory precision and acceptable computation intensity [28]. As can be seen, it contains two parallel RC networks connected in series topology to characterize the battery polarization. Based on Figure 1, the circuit equation can be built, as:

$$\dot{\mathcal{U}}\_s = -\frac{\mathcal{U}\_s}{R\_s \mathcal{C}\_s} + \frac{I}{\mathcal{C}\_s} \tag{1}$$

$$
\dot{\mathcal{U}}\_l = -\frac{\mathcal{U}\_l}{\mathcal{R}\_l \mathcal{C}\_l} + \frac{\mathcal{I}}{\mathcal{C}\_l} \tag{2}
$$

$$
\mathcal{U}\_l = \mathcal{U}\_{\text{occ}} - \mathcal{U}\_s - \mathcal{U}\_l - \mathcal{R}\_t I \tag{3}
$$

where *Rl* and *Cl* indicate the internal resistance and capacitance of electrochemical polarization, while *Rs* and *Cs* denote the internal resistance and capacitance of the concentration polarization.

**Figure 1.** The second-order resistance–capacitance (RC) equivalent circuit model.

*Us* and *Ul* denote the voltage drop across *RsCs* and *RlCl*, respectively; *Ut* indicates the terminal voltage; *I* is the loading current; *Uocv* stands for the open circuit voltage; and *Re* represents the internal ohmic resistance. In addition, the SOC denotes the ratio of available remaining capacity over the rated capacity (the maximum available capacity), as:

$$SOC(t) = SOC(t\_0) - \frac{\int\_{t\_0}^{t} \eta\_c I(t)dt}{Q\_N} \tag{4}$$

where *SOC*(*t*) indicates the SOC value at *t*, respectively; *SOC*(*t*0) stands for the SOC value at *t*0; *QN* is the rated capacity of battery; and η*<sup>c</sup>* represents the columbic efficiency.

#### *2.2. Experiments*

The basic specifications of the battery are illustrated in Table 1. The battery's nominal voltage is 3.6 V, and the nominal capacity is 2.55 Ampere hour (Ah). To characterize the battery's electrical performance, some prerequisite experiments are conducted, including an accelerated aging test, performance test, and dynamic test. The accelerated cycle life aging test is carried out at 25 ◦C, which is divided into seven stages. They are separated by the cycles 0, 30, 60, 90, 120, 150 and 180 (defined as cyc0, cyc30, cyc60, cyc90, cyc120, cyc150 and cyc180 hereinafter). The battery cell is charged by means of the constant current-constant voltage (CCCV) scheme with the current rate of 1 C, and discharged by means of constant current (CC), with the current rate of 2 C in each cycle. Here, C denotes the rated capacity value of the battery with the unit of Ah. The performance tests, including the capacity test and hybrid pulse power characterization (HPPC) test, are carried out periodically during the cycle life test. In addition, a typical dynamic test based on the urban dynamometer driving schedule (UDDS) is executed to verify the performance under dynamic operating conditions. Figure 2 shows the decay variation of the discharge capacity. It can be clearly observed that the discharge capacity basically remains unchanged from cyc0 to cyc30, and it tends to decline faster after 30 cycles, and the amount of electric energy decreases faster as the cycle number increases. After 180 cycles, the maximum discharge capacity decreases from the initial 2.614 Ah to 1.129 Ah, which remains only 44.3% of the nominal capacity. The discharge capacity decreases to 1.97 Ah at the fifth stage (cyc120), which is 77.25% of the nominal capacity. When the capacity drops to 80% of the rated value, the battery should be abandoned; and therefore, only the experimental results of the first five aging stages are applied to analyze and verify the performance of proposed algorithm in this work.

**Table 1.** The specifications of the battery cell.


**Figure 2.** The decay of battery cell discharge capacity under 180 cycles.

Furthermore, battery degradation not only features the capacity reduction, but also embodies the variation of the OCV-SOC relationship. The relationship curve between OCV and SOC at various aging status is presented in Figure 3. It is apparent that the OCV–SOC relationship curve changes gradually as the cycle number increases. When the SOC is more than 20%, the trend of the curves remains basically the same. However, when the SOC is less than 20%, especially under 10%, the OCV changes significantly. Generally, to protect and extend the battery life of EVs, the battery discharge cut-off SOC is usually set to 10% or 20%. When the SOC ranges from 20% to 60%, the OCV at different aging status differs obviously, and the maximum difference value is 20 mV. When the SOC is more than 60%, the OCV values at different aging statuses are relatively close, and the maximum difference is within 10 mV.

**Figure 3.** The relationship curve of the open circuit voltage (OCV) and state of charge (SOC) at different aging stages.

#### **3. The Joint Estimation of SOC and Battery Capacity**

The framework of the joint capacity and SOC co-estimation is shown in Figure 4. It mainly includes four parts: the strategy module, the modeling module, the capacity and parameters estimation module, as well as the SOC estimation module. First, the strategy module starts to accumulate the experimental data until the length of data is more than the preset threshold. Then, the capacity estimation module employs the GA to conduct the parameter identification and capacity estimation, based on the acquired data and the established model. After finding the model parameters and capacity, the SOC estimation module is triggered to estimate the SOC based on the SRCKF. Note that the modeling and parameters estimation modules are not invoked every time. The detailed co-estimation procedure is elaborated in the following.

#### *3.1. The Capacity Estimation Algorithm*

The battery capacity is deemed to be a significant parameter that needs to be identified. Firstly, based on the OCV–SOC curves illustrated in Figure 3, a TDRS with respect to the capacity, SOC and OCV, is constructed, as plotted in Figure 5. Next, the TDRS is imported into the established battery model. Finally, the usable battery capacity is incorporated into the battery model parameters, and is identified by the GA [29].

**Figure 4.** The scheme of the SOC and battery capacity co-estimation algorithm for lithium–ion batteries.

**Figure 5.** The three-dimensional response surface via capacity–SOC–OCV.

By this manner, the problem of usable battery capacity estimation can be transformed into the problem of searching the optimal OCV–SOC relationship match on the built TDRS by applying the optimization algorithm. It is worth noting that the ambient temperature plays an important influence on battery capacity. However, the temperature influence is not taken into account in this study, as the battery system onboard is generally equipped with a good thermal management system, thereby ensuring the temperature variation is within ±5 ◦C [30].

In consideration of accuracy and complexity, a fifth-order polynomial function is selected to describe the relationship among OCV, SOC and capacity, as:

$$\text{dL}\_{\text{OCV}}(\text{SOC}, \mathbb{C}\_{\mathfrak{u},i}) = a\_{1,i} \times \text{SOC}^5 + a\_{2,i} \times \text{SOC}^4 + a\_{3,i} \times \text{SOC}^3 + a\_{4,i} \times \text{SOC}^2 + a\_{5,i} \times \text{SOC} + a\_{6,i} \tag{5}$$

where *Ca*,*<sup>i</sup>* represents the available battery capacity at the *i*th capacity point.

α*j*,*i*(*j* = 1, ··· , 6) denotes the fitting coefficient of OCV and SOC at the *i*th capacity point, which is no longer a constant, and is herein defined as a quadratic function of *Ca*,*i*, as:

$$a\_{j,i} = b\_{2,i} \mathbb{C}\_{a,i}^2 + b\_{1,i} \mathbb{C}\_{a,i} + b\_{0,i} \tag{6}$$

where *bt*,*i*(*t* = 0, 1, 2) denotes the capacity coefficient. The above equation can be rewritten into a matrix form, as:

$$\left[\left[\alpha\_{1,i}, \alpha\_{2,i}, \alpha\_{3,i}, \alpha\_{4,i}, \alpha\_{5,i}, \alpha\_{6,i}\right]^T = \Gamma \times \left[\begin{array}{cc} \mathbb{C}^2\_{a,i} & \mathbb{C}\_{a,i} & 1 \end{array}\right]^T \tag{7}$$

where Γ refers to a 6 × 3 capacity coefficient matrix, which can be obtained by the polynomial fitting, and the results are described in Equation (8).

$$
\Gamma = \begin{bmatrix}
14.2473 & -57.3965 & 61.0568 \\
38.9053 & -151.7664 & 161.0539 \\
3.0020 & -10.7412 & 11.42.84 \\
\end{bmatrix} \tag{8}
$$

Now, according to Equations (5)–(8), a nonlinear relationship can be built among OCV, capacity and SOC. Once the SOC is determined, it will be a deterministic mapping function between OCV and SOC. Together with Equations (1)–(3), the capacity identification can be conducted simultaneously with other parameters, including those of the RC networks. During the parameters identification, the variation of model parameters is eventually reflected by the difference between the estimation result and the terminal voltage. Based on the OCV model, the discrete mathematical expression of the second-order RC-based ECM can be reformulated as:

$$\begin{split} \mathcal{U}l\_{l,k} &= \mathcal{U}\_{\text{occ},k}(\text{SOC}\_{k}, \mathbb{C}\_{\mathfrak{a}}) - \left\{ \exp(-\Delta t/\tau\_{s}) \mathcal{U}\_{s,k} + \mathcal{R}\_{\mathfrak{s}} [1 - \exp(-\Delta t/\tau\_{s})] I\_{k} \right\} \\ &- \left\{ \exp(-\Delta t/\tau\_{l}) \mathcal{U}\_{l,k} + \mathcal{R}\_{\mathfrak{l}} [1 - \exp(-\Delta t/\tau\_{l})] I\_{k} \right\} - I\_{k} \mathcal{R}\_{\mathfrak{c}} \end{split} \tag{9}$$

where Δ*t* indicates the time interval, and both τ*<sup>s</sup>* = *RsCs* and τ*<sup>l</sup>* = *RlCl* belong to time constants. As can be seen from Equation (9), the TDRS based on the OCV model in Equation (5) is imported into the ECM. Hence, the difference of relationship between SOC and OCV at different capacity levels is eventually highlighted by the estimation results of the terminal voltage. In this manner, the battery capacity can be added into the model parameter series for identification. To attain it, the GA is employed to find the optimal combination of model parameters and capacity, in which the optimal parameter group to be identified can be expressed as:

$$\boldsymbol{\Theta}\_{\text{optimal}} = \begin{bmatrix} \boldsymbol{R}\_{\text{c}\prime} \, \boldsymbol{R}\_{\text{s}\prime} \, \boldsymbol{R}\_{l\prime} \, \mathbf{C}\_{\text{s}\prime} \, \mathbf{C}\_{l\prime} \, \mathbf{C}\_{a} \end{bmatrix} \tag{10}$$

During the identification process, the minimum root mean squared error (RMSE) of the terminal voltage is taken as the fitness function, as:

$$\begin{cases} \,\,\rho(\boldsymbol{\theta}\_{\text{optimal}}) = \min[RMSE] \\\\ \,\,\,RMSE = \sqrt{\frac{\sum\_{k}^{N} \left[ \boldsymbol{\mu}\_{t,k} - \boldsymbol{\Omega}\_{t,k}(\boldsymbol{\theta}\_{\text{optimal}}) \right]^{2}}{N}} \end{cases} \tag{11}$$

where *Ut*,*<sup>k</sup>* and *U*ˆ *<sup>t</sup>*,*k*(θ*optimal*) represent the measured terminal voltage and the estimated terminal voltage at step *k*, respectively.

In addition, the constraints of optimization algorithm are subject to:

$$\begin{cases} \begin{cases} 20\% \le \text{SOC} \le 100\% \\\ 0.005 \text{ } \Omega \le R\_{\varepsilon}, R\_{s}, R\_{l} \le 0.1 \text{ } \Omega \end{cases} \\\ \begin{cases} 0.5 \text{ s} \le \tau\_{s} \le 1000 \text{ s} \\\ 0.5 \text{ s} \le \tau\_{l} \le 1000 \text{ s} \end{cases} \end{cases} \tag{12}$$

The setting of these constraints is explained as follows. In practical applications, to avoid the over-charge and over-discharge of the battery, the range of the SOC is generally set to 20% to 100% for guaranteeing proper operation and extending the service life of batteries [31]. The upper limits of internal ohmic resistance *Re*, internal resistance *Rl* of electrochemical polarization and internal resistance *Rs* of the concentration polarization are determined in terms of the specifications of batteries and the technical parameters supplied by the manufacturers. Their low limits are all determined to be 0.005 Ω, based on the parameter calculation of ECM introduced in [32], as well as the experimental analysis. Meanwhile, the range of *Cl* and *Cs* can be deduced to be 100 F to 10<sup>4</sup> F. In addition, τ*<sup>l</sup>* and τ*<sup>s</sup>* are time constants, where τ*<sup>l</sup>* = *RlCl* and τ*<sup>s</sup>* = *RsCs*. Hence, the range of τ*<sup>l</sup>* and τ*<sup>s</sup>* can be limited with 0.5 s to 1000 s. In summary, when the above-mentioned battery model parameters and capacity are identified, these parameters will be transmitted into the SOC estimation module. However, it is worth noting that the parameter identification based on the GA requires a certain amount of data to obtain an ideal identification result. Therefore, the capacity estimation method proposed in this paper only runs when the data length reaches a pre-set condition, and the determination of data length will be discussed in the next section.

#### *3.2. The SOC Estimation Algorithm*

After obtaining the model parameters and battery capacity, the SRCKF is proposed to attain the estimation of battery SOC with the cyclic recurrence based on the established second-order RC ECM. In comparison with the traditional cubature Kalman filter (CKF), the SRCKF can directly perform iterative update in the form of calculating the square-roots of the covariance matrices during the filtering process, which determines the non-negative definite value of the covariance matrix, and avoids the divergence of filter [33]. In general, a discrete nonlinear dynamic system with enhanced noise can be modeled, as:

$$\begin{cases} \mathbf{x}\_{k+1} = f(\mathbf{x}\_k, \mathbf{u}\_k) + w\_k \\ \mathbf{z}\_k = h(\mathbf{x}\_k, \mathbf{u}\_k) + v\_k \end{cases} \tag{13}$$

where *xk* ∈ *Rn* and *zk* indicate the system state vector and the system output at time *k*, respectively. *f*(·) and *h*(·) denote the nonlinear system state function and nonlinear measurement function, respectively. *wk* stands for random process noise indicating uncertain input. *vk* denotes the observation noise, which is generally employed to simulate sensor noise affecting the output measurement. Additionally, the corresponding covariance of *wk* and *vk* are *Qk* and *Rk*, respectively. Based on the established ECM, the time-discrete state equation and measurement equation can be respectively expressed, as:

$$
\begin{bmatrix}
\mathcal{U}\_{k,k+1} \\
\mathcal{U}\_{l,k+1} \\
\mathcal{SOC}\_{k+1}
\end{bmatrix} = \begin{bmatrix}
\exp(-\Delta t/\tau\_{s}) & 0 & 0 \\
0 & \exp(-\Delta t/\tau\_{l}) & 0 \\
0 & 0 & 1
\end{bmatrix} \begin{bmatrix}
\mathcal{U}\_{s,k} \\
\mathcal{U}\_{l,k} \\
\mathcal{SOC}\_{k}
\end{bmatrix} + \begin{bmatrix}
R\_{c}(1-\exp(-\Delta t/\tau\_{s})) \\
R\_{l}(1-\exp(-\Delta t/\tau\_{l})) \\
\end{bmatrix} \begin{bmatrix}
\mathcal{U}\_{k,k} \\
w\_{2,k} \\
w\_{3,k}
\end{bmatrix} \tag{14}
$$

$$
\mathcal{U}\_{l,k} = \mathcal{U}\_{\text{arc}\underline{k}}(\mathcal{SOC}\_{k}, \mathcal{C}\_{d}) + \begin{bmatrix}
\end{bmatrix} \begin{bmatrix}
\mathcal{U}\_{s,k} \\
\mathcal{U}\_{l,k} \\
\mathcal{SOC}\_{k}
\end{bmatrix} + \begin{bmatrix}
\end{bmatrix} \tag{15}
$$

where the system state variable *xk* <sup>=</sup> *Us*,*<sup>k</sup> Ul*,*<sup>k</sup> SOCk T* , input variable *uk* = *Ik* and system output *zk* = *Ut*,*k*. In this study, The SRCKF algorithm is adopted to estimate the SOC, of which the general process is summarized in Table 2, where *n* is the state dimension, and *m* denotes the total number of volume points, which number is twice those of the state dimension. The sample [1] indicates a complete set of fully symmetric points, of which the set of points is obtained through the complete permutation of elements of the n-dimensional unit vector *e* = [1, 0 ··· 0] *<sup>T</sup>* and the alteration of the element symbol. [1]*<sup>g</sup>* represents that the point is centered at the *g*th point of [1]. *x*ˆ*<sup>k</sup>* and *z*ˆ*<sup>k</sup>* are the predicted state and measurement, respectively. *SQ*,*k*−<sup>1</sup> and *SR*,*<sup>k</sup>* denote the square-roots of the process noise covariance matrix *Qk*−<sup>1</sup> and the measurement noise covariance matrix *Rk*, respectively.

**Table 2.** The process of SOC estimation based on the square root cubature Kalman filter (SRCKF) algorithm.

(a) Initialization:

$$\begin{cases} \begin{aligned} \hat{\mathbb{x}}\_{\text{00}} &= E[\mathbb{x}\_{\text{0}}] \\ P\_{\text{00}} &= E\left[ (\mathbb{x}\_{\text{0}} - \mathbb{x}\_{\text{00}}) (\mathbb{x}\_{\text{0}} - \mathbb{x}\_{\text{0}\text{0}})^{T} \right] \end{aligned} \tag{16}$$

Determine the initial value *S*0|<sup>0</sup> of the square roots of the error covariance matrix by the Cholesky decomposition:

$$S\_{\rm{\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tag}} = \left[\text{clol}(P\_{\rm{\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tiny\tag}})\right]^{T} \tag{17}$$

(b) Calculate the basic cubature points and weight:

$$\xi\_{\mathcal{S}} = \sqrt{\frac{m}{2}} [1]\_{\mathcal{S}'} (\mathcal{S} = 1, 2, \cdots, m) \tag{18}$$

(c) Iteration:

for *k* = 1, 2, ··· , *N*

Time update:

Step 1: calculate the cubature points:

$$X\_{\mathcal{G},k-1|k-1} = S\_{k-1|k-1} \xi\_l + \pounds\_{k-1|k-1} \tag{19}$$

Step 2: calculate the propagated cubature points:

$$X\_{\mathbf{g},k|k=1}^{\*} = f(X\_{\mathbf{g},k-1|k=1}, \mu\_k) \tag{20}$$

Step 3: calculate the predicted state:

$$\pounds\_{k|k-1} = \frac{1}{m} \sum\_{\mathcal{S}^{-1}}^{m} X^\*\_{\mathcal{S},k|k-1} \tag{21}$$

Step 4: calculate the state-weighted center matrix:

$$\chi\_{k|k-1}^\* = \frac{1}{\sqrt{m}} \left[ \begin{array}{cccc} X^\*\_{1,k|k-1} - \p\_{k|k-1} & \cdots & X^\*\_{m,k|k-1} - \p\_{k|k-1} \end{array} \right] \tag{22}$$

Step 5: calculate the square-root of the prediction error covariance matrix:

$$S\_{k|k-1} = \operatorname{Train}\left(\left[\chi^\*\_{k|k-1}, S\_{Q,k-1}\right]\right) \tag{23}$$

Measurement update:

Step 1: recalculate the cubature points:

$$X\_{\mathcal{G},k\mathbb{K}-1} = S\_{k\mathbb{K}-1}\xi\_{\mathcal{G}} + \mathfrak{A}\_{k\mathbb{K}-1} \tag{24}$$

Step 2: update the propagated measurement cubature points:

$$Z\_{\mathbb{S},k\mathbb{k}-1} = h(X\_{\mathbb{S},k\mathbb{k}-1}, \mu\_k) \tag{25}$$

Step 3: estimate the predicted measurement:

$$\mathfrak{A}\_{k|k-1} = \frac{1}{m} \sum\_{\mathcal{S}=1}^{m} Z\_{\mathcal{S},k|k-1} \tag{26}$$

Step 4: evaluate the measurement-weighted center matrix:

$$Z\_{k|k-1} = \frac{1}{\sqrt{m}} \left[ \begin{array}{cccc} Z\_{1,k|k-1} - \mathbb{A}\_{k|k-1} & \cdots & Z\_{m,k|k-1} - \mathbb{A}\_{k|k-1} \end{array} \right] \tag{27}$$

Step 5: estimate the square root of the innovation covariance matrix:

$$S\_{\mathbb{Z}\mathbb{Z},k|k-1} = \operatorname{Tri}\left(\left[\mathbb{Q}\_{k|k-1}, \mathbb{S}\_{\mathbb{R},k}\right]\right) \tag{28}$$

Step 6: update the state-weighted center matrix:

$$\chi\_{k|k-1} = \frac{1}{\sqrt{m}} \left[ \begin{array}{cccc} X\_{1,k|k-1} - \p\_{k|k-1} & \cdots & X\_{m,k|k-1} - \p\_{k|k-1} \end{array} \right] \tag{29}$$

Step 7: estimate the cross-covariance matrix:

$$P\_{\chi z \downarrow k \vert k-1} = \chi\_{k \vert k-1} \zeta\_{\chi\_{\parallel k-1}}^T \tag{30}$$

Moreover, update the Kalman gain, state and square root of the error covariance Step 1: estimate the Kalman gain matrix:

$$\mathcal{W}\_{k} = \frac{P\_{xz,k\mathbb{j}k-1}/S\_{zz,k\mathbb{j}k-1}^{T}}{S\_{zz,k\mathbb{j}k-1}} \tag{31}$$

Step 2: estimate the final updated state:

$$\mathbb{L}\pounds\_{k\nparallel} = \pounds\_{k\nparallel-1} + \mathcal{W}\_k(z\_k - \mathcal{Q}\_{k\nparallel-1})\tag{32}$$

Step 3: update the corresponding square-root of the error covariance matrix:

$$S\_{k\&k} = \operatorname{Trid}\Big(\Big[\chi\_{k\&-1} - \mathcal{W}\_k \zeta\_{k\&-1}, \mathcal{W}\_k S\_{R,k}\Big]\Big)\tag{33}$$

End

#### **4. Verification and Discussion**

#### *4.1. Verification Study on Di*ff*erent Data Lengths*

In this section, different data lengths are selected to investigate the effectiveness of the proposed co-estimation algorithm. Actually, the experimental data at different aging stages can be chosen to verify the SOC estimation of different data lengths. Nonetheless, based on the estimated capacity, the estimation error of battery capacity is maximum at cyc30. Hence, to better verify the performance of the proposed estimation algorithm, the battery after being cycled 30 times is chosen as the test target, and the current schedules acquired based on the UDDS experiment are repetitively operated until the terminal voltage reaches the cut-off voltage designated by the manufacturer. Note that when the data length is less than the pre-set threshold value, the SOC module still uses the previously identified capacity and parameters to conduct the estimation.

Figure 6 shows the current profiles under the UDDS experiment. It can be clearly found that the entire discharging process takes around 554 min. To evaluate the influence incurred by different data lengths when identifying the model parameters, the data with the duration of 65, 84, 130 and 200 min (defined as 65 min, 84 min, 130 min and 200 min, respectively) are randomly selected as the test target, and the remaining data are applied for SOC estimation. Note that when the data length is 554 min (total loading profile process), the SOC is estimated without the update of the capacity value and parameters. Table 3 compares the battery capacity estimation results with respect to different data lengths.

**Figure 6.** The urban dynamometer driving schedule (UDDS) current.

**Table 3.** The estimation results and errors of the battery capacity corresponding to various data lengths.


The actual capacity 2.5321 Ah is measured through the calibration test, and the estimated capacity ranges from 2.4366 Ah to 2.5065 Ah. It can be observed that data duration shows certain influence on the estimation results, and the estimated error of the battery capacity decreases by 2.7606%, from 3.7716% to 1.011%, after increasing the data length.

Based on the obtained capacity, the detailed results of the SOC and estimation error with different identification data lengths are presented in Figure 7, and the statistic results are provided in Figure 8. As demonstrated in Figure 7, when the date length increases from 54 min to 200 min, the estimated SOC can quickly converge to the reference value according to the updated parameters and capacity, and the maximum absolute error decreases from 3.643% to 0.989%. When the data length reaches 200 min, the estimation errors are restricted within a small range, less than 1%. Besides, Figure 7 also shows the SOC estimation results without considering the capacity's update. The maximum absolute error, the mean absolute error and the RMSE are 2.538%, 1.661% and 1.737%, respectively. It is apparent that the estimated SOC looks more divergent without the capacity update, thereby manifesting the advance of the joint estimation algorithm. From Figure 8, we can find that when the data length increases from 65 min to 200 min, the maximum absolute error, mean absolute error and RMSE decrease from 3.643%, 2.331%, 2.498% to 0.989%, 0.234%, 0.319%, respectively. The results demonstrate that as the calculated data length increases, the estimated SOC becomes closer to the reference value, and this is mainly because the GA shows a global optimization ability, and when more input data samples are referred, the prediction results will be more accurate. Hence, appropriately increasing the duration of data is beneficial for improving the accuracy of capacity identification and SOC estimation. Nonetheless, it is appreciably time-consuming when increasing the amount of data to estimate the battery capacity. To balance the relationship between error and calculation time, the data duration of 200 min is considered as the preferred length. In the following, the estimated results with the data length of 200 min are all adopted for SOC estimation and comparison.

**Figure 7.** The results of SOC estimation based on various data lengths: (**a**) SOC estimation results; (**b**) SOC error.

**Figure 8.** Comparison of the battery SOC estimation results of different data lengths.

#### *4.2. SOC Estimation under Various Degradation Stages*

To verify the feasibility of the proposed co-estimation scheme, the battery cells are experimentally and circularly tested with the UDDS current at different aging levels. According to the estimation algorithm of capacity addressed previously, the pre-set data length is 200 min. Table 4 and Figure 9 compare the estimated results of battery capacity at different aging stages (fresh, nearly fresh, slightly cycled, severely cycled and lifespan exceeded), of which the number of cycles ranges from 0 to 120, with 30 as the interval. As illustrated in Table 4, the battery capacity declines with the cycling operation. The proposed algorithm enables that the maximum relative and absolute errors are less than 1.011% and 0.026 Ah when the battery is cycled for 30 times, thereby indicating its preferable capability of estimating the battery capacity at different aging statuses. Furthermore, the estimated results can also commendably reflect the decay trend of battery capacity.


**Table 4.** The estimated results and errors of the battery during the entire lifespan.

After finding the model parameters including the capacity value, the estimated SOC results at different aging status are demonstrated in Figure 10, and the statistic results are summarized in Figure 11. As Figure 10 suggests, it is obvious that when the battery ages, the total discharging time gradually decreases under the same operating conditions. The initial SOC is 20%, with the error of 80%, and the estimated SOC at various aged status can all converge to the reference values. Figure 10 also reveals that the maximum absolute error of SOC is restricted within 1% after the correction of the initial

SOC error, even when the battery is aged, and the convergence time is less than 120 s. As discussed previously, the accuracy of SOC estimation is heavily influenced by the battery degradation. Without considering the update of battery capacity, the SOC estimation error increases towards higher numbers of cycles. In this study, the updated battery capacity is exploited to assist in improving the SOC estimation in the entire lifespan. As can be found in Figure 11, the maximum value of absolute error, mean absolute error and the RMSE are 0.987%, 0.484% and 0.566%, respectively, occurring in cyc30. The reason is that when the cycle number is 30, the estimation error of capacity reaches 1.011%, thus leading to the worst SOC estimation. Even so, the maximum absolute error of SOC is still restricted within 1.1% after the correction of initial SOC. As the number of battery cycles increases, the estimated SOC error does not increase obviously, manifesting that the updated capacity value contributes to the SOC estimation. From this point of view, regular updates of battery capacity in the aging process are imperative to improve the accuracy of SOC estimation.

**Figure 9.** Capacity value of measurement and estimation at various degradation extents. (**a**) Estimation results; (**b**) Relative error.

**Figure 10.** The results of SOC estimation with the aged battery cell: (**a**) SOC estimation results; (**b**) SOC error.

**Figure 11.** The estimated errors of the SOC corresponding to the aging battery cell.

#### **5. Conclusions**

In this study, a model-based adaptive joint estimation algorithm of SOC and capacity is proposed for lithium–ion batteries. The SOC estimation is implemented based on a second-order ECM, with the SRCKF algorithm considering the capacity degradation and parameters variation. The battery capacity is imported into the model parameter group, and it is jointly identified by the GA and the constructed TDRS. After obtaining the parameters and the capacity, the SOC is accurately estimated by the SRCKF. Through the experimental validations in terms of different degradation status, varying duration of recorded data and various dynamic operating conditions, the preferable performance of the proposed method is satisfactorily verified. The experimental results elucidate that the co-estimation approach can improve the SOC estimation accuracy in the entire battery lifespan cycle with the update of capacity, even in the cases of aged batteries and under complicated operating conditions.

In addition, this paper only investigates the SOC and capacity estimation for battery cells. However, the capacity and SOC of battery packs are also particularly critical in practical applications, and they will certainly be our research focus in the future.

**Author Contributions:** Conceptualization, Z.C. and J.X.; methodology, Z.C.; software, J.X.; validation, Z.C., J.X., X.S., S.S., and Y.L.; formal analysis, J.X.; investigation, Z.C.; resources, X.S.; data curation, J.X., X.S., and J.S.; writing—original draft preparation, Z.C., J.X., X.S. and J.S.; writing—review and editing, Z.C.; visualization, J.X.; supervision, Z.C.; project administration, Z.C.; funding acquisition, Z.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the National Science Foundation under Grant 61763021 and Grant 51775063, in part by the National Key R&D Program of China under Grant 2018YFB0104000, and in part by the EU-funded Marie Skłodowska-Curie Individual Fellowships Project under Grant 845102-HOEMEV-H2020-MSCA-IF-2018 in part.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Dual Nonlinear Kalman Filter-Based SoC and Remaining Capacity Estimation for an Electric Scooter Li-NMC Battery Pack**

#### **Filip Maleti´c \*, Mario Hrgeti´c and Joško Deur**

Department of Robotics and Automation of Manufacturing Systems, Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, 10000 Zagreb, Croatia; mario.hrgetic@fsb.hr (M.H.); josko.deur@fsb.hr (J.D.)

**\*** Correspondence: filip.maletic@fsb.hr; Tel.: +385-1-6168-555

Received: 23 December 2019; Accepted: 17 January 2020; Published: 22 January 2020

**Abstract:** Accurate, real-time estimation of battery state-of-charge (SoC) and state-of-health represents a crucial task of modern battery management systems. Due to nonlinear and battery degradation-dependent behavior of output voltage, the design of these estimation algorithms should be based on nonlinear parameter-varying models. The paper first describes the experimental setup that consists of commercially available electric scooter equipped with telemetry measurement equipment. Next, dual extended Kalman filter-based (DEKF) estimator of battery SoC, internal resistances, and parameters of open-circuit voltage (OCV) vs. SoC characteristic is presented under the assumption of fixed polarization time constant vs. SoC characteristic. The DEKF is upgraded with an adaptation mechanism to capture the battery OCV hysteresis without explicitly modelling it. Parameterization of an explicit hysteresis model and its inclusion in the DEKF is also considered. Finally, a slow time scale, sigma-point Kalman filter-based capacity estimator is designed and inter-coupled with the DEKF. A convergence detection algorithm is proposed to ensure that the two estimators are coupled automatically only after the capacity estimate has converged. The overall estimator performance is experimentally validated for real electric scooter driving cycles.

**Keywords:** electric vehicle; lithium-ion battery; estimation; Kalman filter; state-of-charge; state-of-health; resistance; open-circuit voltage; battery capacity

#### **1. Introduction**

Modern battery management systems (BMSs), among other functionalities, include a number of algorithms for estimating key battery state variables such as state-of-charge (SoC) and remaining available charge capacity, and model parameters such as internal resistance [1]. The SoC estimate can be used for predicting the current vehicle range, as well as for identification of current battery operating point which is important from the standpoint of ensuring battery safety. On the other hand, the internal resistance and capacity estimates are the main indicators used for tracking the battery degradation level, i.e., estimation of battery state-of-health (SoH) [2]. Furthermore, almost every battery model parameter is changing with battery degradation, so that for robust SoC and SoH estimation, those changes should be accurately tracked, as well.

Battery state and parameter estimation algorithms are often based on Kalman filters (KF), which in its basic linear version represent an optimal recursive solution for estimating hidden states of a linear, time-varying Gaussian system (i.e., probabilistic inference) [3]. While the Gaussian assumption holds in many cases based on the central limit theorem, the battery model is inherently nonlinear, which calls for application of nonlinear KF forms. Two of the most widely used nonlinear KFs are extended Kalman filter (EKF) and sigma-point Kalman filter (SPKF) [3]. The EKF relies on analytical linearization

of the model around a time-varying operating point (i.e., an expected value of the estimated random state), while SPKF statistically linearizes the model around several operating points (depending on the number of states that are estimated).

Topic of state and parameter estimation of Li-ion battery cells has been addressed by many previous studies. One of the first implementations of dual extended Kalman filter-based (DEKF) estimator of SoC and resistance parameters can be found in [4]. Researchers have been upgrading the estimators ever since, e.g., using an adaptation mechanism for process noise variance recalculation [5], or applying more advanced filters such as SPKF [6] or particle filter (PF) [7]. These approaches are based on the assumption of constant model parameters/characteristics such as the SoC-dependent open-circuit voltage (OCV) characteristic *Uoc*(*SoC*) or battery remaining capacity. Since those parameters are in fact dependent on SoH [8] and temperature [9], they should be estimated as well, for accurate and robust overall estimation.

There are several studies that account for *Uoc*(*SoC*) variation with SoH by implementing the offline identified response surface model of *Uoc* with respect to SoC and remaining capacity [10–12]. Authors in [13] use the model migration method to adapt an offline trained model. An obvious disadvantage of this approach is related to the need of having a large data set from previously conducted aging experiments on the same cell type, as well as lack of temperature dependency in the response surface model. This disadvantage is tackled in this paper by describing the characteristic *Uoc*(*SoC*) with a model whose parameters are estimated along the rest of model states and parameters within the DEKF structure. Moreover, this approach includes an adaptation mechanism of *Uoc*(*SoC*) which allows for identification of *Uoc*(*SoC*) hysteresis profile.

Remaining capacity estimators based on EKF and PF can be found in [14], while a recursive approximate least-squares approach is proposed in [15]. In both cases the characteristic *Uoc*(*SoC*) is again considered as a constant-parameter dependence. Dual estimation of SoC and capacity can be found in [16], where authors use multiscale estimation with the online identified model, which can be regarded as a next step towards complete estimator. Certain weaknesses of that approach include: (i) Still an offline identified *Uoc*(*SoC*) map is used, (ii) capacity estimate shows considerable variations in steady state, and (iii) the capacity estimator needs to be turned on manually after 25 min in order to ensure overall estimator stability. The multiscale estimator presented in this paper improves the capacity estimation accuracy and flexibility by using a more accurate SPKF and automated turning on the capacity estimator by means of applying a convergence detection algorithm.

Finally, a fully-electric scooter-based experimental verification of the proposed battery estimators is conducted, including consideration of different temperature operating points.

#### **2. Experimental Setup**

The experimental setup includes the fully-electric scooter Govecs S2.6+, powered by the 3.3 kW BLDC electric motor and the battery pack of 400 Li-Ni0.33Mn0.33Co0.33O2 (Li-NMC) cells, connected in the 20 × 20 matrix [17], with the total nominal voltage of 72 V, and the total energy capacity of 4.1 kWh. Battery pack is equipped with BMS which provides basic battery measurements and estimates accessible through the scooter CAN bus.

Electric scooters became an attractive transportation solution in urban areas with mild climate conditions, thus contributing to the current transport electrification effort aimed at reducing traffic congestion, and air and noise pollution. There are already several strong electric scooter manufacturers in the EU (and worldwide), e.g., Govecs, Ujet, Hrowin, Torrot, etc. NMC-type Li-ion batteries represent a preferred energy storage solution in scooter applications [17], because they offer favorable energy density, while not experiencing high loads (in terms of battery C-rate) and not operating in extreme, particularly low temperature conditions, in those applications.

For the research purposes, the scooter has been equipped with the measurement and telemetry system illustrated in Figure 1. The system is built around the Artronic SkyTrack telemetry module, custom-programmed for the acquisition and storage of measurement data, as well as for communication

with the server through GPRS connection in real time. The measurement system consists of voltage and current measurement on battery output nodes, and acquisition of data available from the scooter CAN bus. The battery current is measured by using a precise, low-offset current transducer (LEM CAB 300, [18]), while the battery voltage is measured through a 12-bit analogue input of the telemetry module. Those two measurement values are sampled every 0.1 s and stored in the module. Selected values from the scooter internal CAN bus, such as battery voltage, current and temperature, vehicle's distance travelled, motor on/off flag, as well as the vehicle's current GPS coordinates and longitudinal velocity are stored with the sample rate of 1 s. GPRS connection is used to send data relevant for real-time tracking of scooter, such as its GPS coordinates, battery SoC, and other diagnostic parameters. The whole measurement dataset, including the fast current and voltage measurements, is stored in the telemetry module memory card and can be occasionally downloaded through USB connection to a local PC.

**Figure 1.** Scooter measurement and telemetry system.

#### **3. Battery Pack Model**

This section first presents a battery mathematical model used as a basis for SoC estimator design. Next, models employed for estimation of battery internal parameters used by the SoC estimator are presented. Finally, two offline identification experiments are described, which have been conducted to determine battery model parameters that are considered as constant or used in estimator verification.

#### *3.1. Mathematical Model*

The battery model used in this research is based on the equivalent-circuit model (ECM) showed in Figure 2, which consists of (i) a voltage source dependent on the battery SoC, i.e., the OCV characteristic *Uoc*(*SoC*), (ii) an ohmic resistance *Rohm* which models voltage drops in the electrolyte and electrical contacts, and (iii) a single polarization RC term (*Rp* and *Cp*) which models the slow battery dynamics, i.e., diffusion process. It should be noted that the diffusion process is more accurately modelled with the Warburg element [19] which is here avoided due to the complexity, but it can be approximated by a single or more RC elements connected in series (a single RC element is usually used as a good trade-off between simplicity and accuracy [20]). Moreover, note that the polarization resistance *Rp* in Figure 2 models all voltage drops that are not related to the ohmic one, including that related to charge transfer.

**Figure 2.** Battery equivalent circuit model used in the DEKF.

The above ECM can be described by the following discrete-time time-varying state-space mathematical model [4]:

$$
\begin{bmatrix} \text{SoC}(k) \\ \boldsymbol{i}\_{p}(k) \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & e^{-\frac{T\_{\rm H}}{\tau\_{p}(\text{ScC}(k-1))}} \end{bmatrix} \begin{bmatrix} \text{SoC}(k-1) \\ \boldsymbol{i}\_{p}(k-1) \end{bmatrix} + \begin{bmatrix} -\frac{T\_{\rm H}}{\tilde{\mathbf{C}}\_{\tau\_{p}}} \\ 1 - e^{-\frac{T\_{\rm H}}{\tau\_{p}(\text{ScC}(k-1))}} \end{bmatrix} \dot{\boldsymbol{i}}\_{b}(k-1) \tag{1}
$$

$$\mathcal{L}l\_b(k) = \mathcal{L}l\_{\text{oc}}(\text{SoC}(k)) - \mathcal{R}\_{\text{olm}}(k)i\_b(k) - \mathcal{R}\_{\mathcal{P}}(k)i\_{\mathcal{P}}(k) \tag{2}$$

where *Tu* is the filter sampling time, *Cn* is the battery capacity, τ*<sup>p</sup>* = *RpCp* is the polarization term time constant, and *k* is discrete sample step.

#### *3.2. OCV Model*

Since the battery OCV is a nonlinear function of SoC, and to a lower extent temperature [9] and SoH [8], it is desirable to describe it using a parametric model such as the one used in [4]:

$$\text{LIL}\_{\text{G}}(\text{SoC}) = \begin{bmatrix} K\_0 & K\_1 & K\_2 & K\_3 & K\_4 \end{bmatrix} \begin{bmatrix} 1 & -\frac{1}{\text{S}\text{cC}} & -\text{SoC} & \ln(\text{SoC}) & \ln(1 - \text{SoC}) \end{bmatrix}^T = k\_{\text{0c}} \mathbf{x}\_{\text{0c}} \tag{3}$$

where vector *koc* contains *Uoc*-model parameters that need to be estimated.

#### *3.3. Model of Internal Resistance Parameters*

The presented ECM has two resistance parameters in its model. Both of those resistances are known to depend on SoH and temperature [21]. So, it is important to have them estimated along with the model states. Since there is no resistance model feasible for online estimator implementation, resistances are modelled as random-walk variables:

$$
\begin{bmatrix} R\_{\text{olm}}(k) \\ R\_{\text{P}}(k) \end{bmatrix} = \mathbf{I} \cdot \begin{bmatrix} R\_{\text{olm}}(k-1) \\ R\_{\text{P}}(k-1) \end{bmatrix} + \mathbf{r} \tag{4}
$$

where *I* is the identity matrix, and *r* is the vector containing variances of both resistances. Other variable model parameters, such as those from Equation (3), can be modelled using this approach, as well.

#### *3.4. Identification Experiments*

The battery model parameters that are assumed to be constant or used in estimator verification should be determined by means of specific (targeted) offline identification tests.

#### 3.4.1. Battery OCV Curve

The curve *Uoc*(*SoC*) has been identified during low- and constant-load experiments (during which the vehicle was in rest, while only considerable battery load was scooter headlight), in which case any voltage drop in the battery can be neglected due to the low current (~*C*/50), so that the measured voltage *Ub* can be taken as the OCV *Uoc*. The SoC was estimated by Coulomb counting, i.e., by integrating the measured current. The battery capacity was also identified in this experiment by integration of measured current during the process of full battery discharge, which gave *Cn* = 49.57 Ah. Graphical illustration of the identification experiment and related results are shown in Figure 3. The identified curve *Uoc*(*SoC*) has been used in validation of *Uoc* estimation results (see next sections).

**Figure 3.** (**a**) Low- and constant-load experiment: Current, voltage, and SoC responses, and illustration of (**b**) capacity and (**c**) *Uoc*(*SoC*) identification.

#### Polarization Time Constant

This parameter can be identified during the battery relaxation periods, i.e., parts of driving cycle where current has dropped to zero and remained equal to zero for at least 15 min. The relaxation transients to be identified were extracted from the voltage response (see Figure 4a,b) and approximated with the ECM model shown in Figure 4c.

The identified values of relaxation time constant τ*p*(*SoC*) are shown in Figure 4d. These values were then approximated with a 3rd-order polynomial in dependence on SoC, and that polynomial was later used for calculation of τ*<sup>p</sup>* at every estimator step based on the current, slowly changing SoC working point.

It is important to note that the polarization time constant can also vary with battery temperature and aging [22,23]. These effects are neglected in the estimator problem formulation in this paper, i.e., parameters of the characteristic τ*p*(*SoC*) are not estimated online. This is motivated by the following main reasons: (i) τ*<sup>p</sup>* is not directly involved in the ECM voltage equation (see Equation (2)), thus making it weakly observable in the proposed estimator design; ii) error in τ*<sup>p</sup>* will cause an error in voltage modelling during the transient periods (i.e., before voltage has relaxed), so that the polarization dynamics may influence estimator accuracy only in transient conditions. As needed, the slow temperature- and aging-influenced polarization dynamics can be accounted for in the estimator design either by extending the τ*<sup>p</sup>* characteristic with the temperature and SOH inputs or by considering τ*<sup>p</sup>* as an additional parameter to be estimated, which is a subject of future work.

**Figure 4.** Illustration of τ*p* identification procedure: (**a**)Measured battery voltage responses, (**b**) extracted voltage relaxation periods, (**c**) illustration of typical Li-ion cell voltage response after the current steps with the approximation equation, taken from [1], and (**d**) τ*p*(*SoC*) identification results.

#### **4. State and Parameter Estimator**

This section deals with design, parametrization, and verification of the SoC estimator. It is designed as a dual state and parameter estimator, thus allowing for accompanying estimation of selected ECM parameters (i.e., the battery internal resistances and OCV parameters). A special attention is devoted to estimation of battery OCV hysteresis based on two complementary approaches (adapting the OCV parameters to current sign change or using an explicit hysteresis model).

#### *4.1. DEKF-Based State and Parameter Estimator*

States and parameters of the ECM are estimated with the DEKF, as a well-known approach in the model-based estimation problems where model states and slowly varying model parameters are to be estimated simultaneously [1]. The DEKF equations are not listed here due to paper size constraints, and they can be found in [24]. DEKF consists of two filters operating in parallel based on the state and parameter models:


where *x* and *u* are the vectors of model states and inputs, respectively, *w* is the vector of state variances (with the corresponding covariance matrix *Qx*), θ is the vector of model parameters with their variances contained in vector *r* (with the corresponding covariance matrix *Q*θ), *h* is the model output function (the same output function is used in both state and parameter models), *y* is the measured model output vector with measurement noise and corresponding covariance matrix denoted by *v* and *R*, respectively. The complete, discrete-time state-space model for simultaneous state and parameter estimation then reads (cf. Equations (1)–(4)):

$$\hat{\mathfrak{X}} = \begin{bmatrix} \mathrm{SoC}(k) \\ i\_p(k) \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & e^{-\frac{T\_y}{\tau\_p(\mathrm{SoC}(k-1))}} \end{bmatrix} \begin{bmatrix} \mathrm{SoC}(k-1) \\ i\_p(k-1) \end{bmatrix} + \begin{bmatrix} -\frac{T\_y}{\tilde{\mathbf{C}}\_v} \\ 1 - e^{-\frac{T\_y}{\tau\_p(\mathrm{SoC}(k-1))}} \end{bmatrix} \langle i\_b(k-1) + w \rangle \tag{5}$$

$$\boldsymbol{\Theta} = \begin{bmatrix} \boldsymbol{R}\_{\text{ohm}}(k) \\ \boldsymbol{R}\_{p}(k) \\ \boldsymbol{k}\_{\text{oc}}^{T}(k) \end{bmatrix} = \boldsymbol{I} \cdot \begin{bmatrix} \boldsymbol{R}\_{\text{ohm}}(k-1) \\ \boldsymbol{R}\_{p}(k-1) \\ \boldsymbol{k}\_{\text{oc}}^{T}(k-1) \end{bmatrix} + \boldsymbol{r} \tag{6}$$

$$y(k) = \mathcal{U}\_b(k) = k\_{\text{oc}} \mathbf{x}\_{\text{oc}} - R\_{\text{olm}}(k) i\_b(k) - R\_p(k) i\_p(k) + \upsilon \tag{7}$$

The state estimator model, given by Equations (6) and (8), is considered linear in state equation under the assumption that the nonlinearity of function τ*p*(*SoC*) can be neglected. The only nonlinearity resides in the output equation of the state estimator, related to the *xoc* term (see Equation (3)), so that an EKF is finally used as a model state estimator. On the other hand, the parameter estimator model, given by Equations (7) and (8), is linear, so that the estimator reduces to KF.

#### *4.2. Estimator Parametrization*

The DEKF needs to be properly parametrized. For instance, appropriate statistic parameters such as process and output noise covariances *Q* and *R* should be determined offline. Polarization time-constant τ*<sup>p</sup>* was assumed to be degradation-invariant and used as the identified SoC-dependent profile (see previous section), while battery capacity was in this case taken as a constant value that was measured as described in the previous section. This section also describes an estimator adaptation mechanism that indirectly compensates for the influence of unmodelled hysteresis of curve *Uoc*(*SoC*).

#### 4.2.1. DEKF Covariance Matrices Parametrization

The measurement variable in the DEKF model is the battery output voltage *Ub* (see Equation (8)). Its measurement noise has been estimated by approximating the voltage measurement error histogram with normal distribution, as shown in Figure 5a. The parameter μ identified in Figure 5a is the voltage noise mean value (expectation), while σ is the standard deviation which, after being squared, yields the measurement covariance *<sup>R</sup>* = <sup>53</sup>·10−<sup>3</sup> 2 mV2. The parameter *Lstat* in Figure 5a is the result of Lilliefors normality test. Further in this paper, we calculate *Lstat* for estimator voltage residuals and compare it to the calculated *Lstat* = 0.0436 of voltage sensor noise (see Figure 5a) to check how similar they are, i.e., how accurate is the estimator.

**Figure 5.** (**a**) Estimated voltage sensor noise, (**b**) amplitude of current sensor noise with respect to measured current, taken from [18].

The process noise relates to the current sensor noise, as can be seen from Equation (6). The current is in this case measured with LEM CAB 300 sensor, whose datasheet specifies a linear relation between measured current and magnitude of its measurement error (see Figure 5b). The standard deviation of current sensor noise can be estimated as a value three times lower than the noise magnitude, and covariance matrix is then the diagonal matrix of current sensor noise variances:

$$\sigma\_x = \frac{1.75}{350} \cdot \frac{1}{3} \cdot i\_b \to \mathbf{Q}\_x = \text{diag}\{\sigma\_x^2, \sigma\_x^2\} = \text{diag}\{\left(\frac{i\_b}{600}\right)^2, \left(\frac{i\_b}{600}\right)^2\} \tag{8}$$

#### 4.2.2. Adaptation Mechanism

The relatively simple battery model given by Equations (1) and (2) does not take into account some secondary, but generally influential effect such as the hysteresis of OCV curve *Uoc*(*SoC*) [25]. Since the hysteresis cannot be directly measured in this case, an adaptation mechanism is introduced in the form of single-step increase of the elements of parameter covariance submatrix *Q*θ[3, 7; 3, 7] when the start or end of charging is detected. This approach allows faster convergence of the *Uoc* parameters (written in *koc*), which abruptly change when the sign of battery current (or SoC derivative) occurs due to the existence of hysteresis of *Uoc*(*SoC*) curve. Note that the battery current for the given scooter changes its sign only when the scooter is exposed to change from normal driving to charging or vice versa, because it does not incorporate regenerative braking.

#### *4.3. Estimation Results*

The presented DEKF was validated based on the recorded scooter real city driving cycle data consisting of seven load cycles (i.e., charge/discharge cycles) lasting for 150 h in total. The obtained estimation results are shown in Figure 6. Since the battery SoC cannot be measured, and there is no fully reliable SoC estimate available, the DEKF accuracy is evaluated by analyzing *a posteriori* voltage residual, i.e., difference between the recorded voltage *Ub* and the voltage calculated from output Equation (8) using *a posteriori* estimated states and parameters. The perfectly accurate filter would reduce the voltage residual to the voltage sensor noise, i.e., the residual mean value, standard deviation, and *Lstat* would be close to the values from Figure 5a.

**Figure 6.** DEKF verification results: (**a**) Voltage residual histogram including normal distribution fit, (**b**,**d**) estimated resistances *Rohm* and *Rp*, (**c**) estimated and recorded *Uoc*(*SoC*) curves for a long set of real-life discharging and charging cycles.

Figure 6a shows the voltage residuals histogram including the corresponding normal distribution fit and its parameters. Residual mean value is low, while standard deviation and *Lstat* are larger than those of the voltage sensor noise. The estimated values of resistances *Rohm* and *Rp* are shown in Figure 6b,d, respectively, vs. SoC and color-mapped with respect to battery temperature. These results point out that both resistances show negative correlation with respect to temperature (note: ρ*X*,*<sup>Y</sup>* stands for correlation coefficient between vectors *X* and *Y*, and are obtained by using the MATLAB function corrcoef), which is expected for the Li-ion cell resistances [21]. As of the correlation with respect to SoC (based on visual inspection of Figure 6b,d), both *Rohm* and *Rp* do not seem to be correlated with SoC, which is an expected result for the particular SoC range, based on the estimator results from the available literature [16,21,26] in which resistances more significantly depend on SoC only at the very low and very large SoC bands. The estimated *Uoc* curves during charging and discharging intervals are shown in Figure 6c, along with the "measured" one adopted from Figure 3c. Evidently, the estimated and "measured" curves are in good agreement, and a relatively small hysteresis is apparent (i.e., the charging and discharging curves do not overlap).

The two sets of estimated *Uoc*(*SoC*) curves from Figure 6c have been averaged and shown as dotted lines in Figure 7a. Half of the difference between those two curves yields the estimate of battery hysteresis voltage which is shown in Figure 7b. The estimated hysteresis voltage trend is in line with the results from the literature (e.g., [25]), except in the low-SoC region (*SoC* < 20%), where the estimated hysteresis is larger than what would be expected based on the literature.

**Figure 7.** (**a**) Replot of Figure 6c with added average values of estimated *Uoc*(*SoC*) curves for charge and discharge periods, (**b**) estimated hysteresis voltage.

Now, when the hysteresis voltage is known, the adaptation mechanism may be omitted, and the hysteresis can be accounted for directly through a proper *Uoc*(*SoC*) model extension. A complex, dynamic hysteresis model [27] is not necessary in this case, because the particular scooter does not support regenerative breaking (i.e., its battery is not exposed to often changes of current sign). A simple, instantaneous hysteresis model can be described by introducing an auxiliary variable *s* described as [4]:

$$s(k) = \begin{cases} 1, \ i\_b(k) > 3\sqrt{Q\_x} \\ -1, \ i\_b(k) < -3\sqrt{Q\_x} \\ s(k-1), \ i\_b(k) < \left| 3\sqrt{Q\_x} \right| \end{cases} \tag{9}$$

(where 3 <sup>√</sup> *Qx* is the current sensor noise amplitude calculated using the current sensor variance from Equation (9)) and using it to modify the output equation (cf. Equation (8)):

$$y(k) = \mathcal{U}\_b(k) = \mathcal{U}\_{\text{oc}}(\text{SoC}(k)) - R\_{\text{shm}}(k)\dot{\imath}\_b(k) - R\_p(k)\dot{\imath}\_p(k) + s(k)\mathcal{M}\_0(k) + \upsilon \tag{10}$$

where *M*<sup>0</sup> is the hysteresis voltage value obtained from data shown in Figure 7b by means of 10th-order approximation polynomial. Described hysteresis is used instead of the adaptation mechanism in the rest of the paper.

#### **5. Battery Capacity Estimation**

This section presents the battery remaining charge capacity estimator, and its integration into the overall SoC and capacity estimation algorithm. The capacity estimator is supplemented with a convergence detection algorithm to perform automatic coupling of the capacity estimator with the SoC estimator after the capacity estimate has converged. Finally, the complete estimation algorithm is verified for real driving battery load cycles.

#### *5.1. Capacity Estimation Model*

Since the battery capacity parameter is not directly involved in the model output equation (i.e., Equation (8)), it is not convenient to estimate it as another random-walk parameter in the DEKF [15]. Instead, the model for capacity estimation could be defined as [15]:

$$\mathcal{C}(k) = \mathcal{C}(k - L) + r\_{\mathcal{C}} \tag{11}$$

$$\text{SoC}(k - L + 1) - \text{SoC}(k) = \frac{T\_{\text{II}}}{\text{C}(k)} \sum\_{j=k-L+1}^{k} i\_{\text{b}}(j) + v\_{\text{SoC}}(k) \tag{12}$$

where *C*(*k*) is capacity, *rC* is random walk noise for capacity parameter model with the corresponding covariance *QC*, *L* is the number of basic (DEKF) sampling steps between two capacity estimates, and *vSoC* is measurement noise of SoC signal difference with the corresponding covariance *RSoC*.

The model output is the SoC difference between two capacity estimates, while its input is the cumulative sum of battery current between those time instances. The SoC, as an output term, cannot be measured, but can be estimated by using the previously designed DEKF (both, estimates of SoC mean value and its variance are available). By looking at Equation (13) it can be seen that capacity estimate cannot be updated at the same rate as DEKF, because the signal-to-noise ratio of SoC estimate would be too low for the SoC dynamics being much slower than the current dynamics. The capacity estimator is therefore executed every *L* time steps, where *L* is in the range of 600–6000, i.e., 1 to 10 min. The overall estimator, i.e., the previously discussed DEKF extended with the capacity estimator, is shown in Figure 8.

**Figure 8.** Overall algorithm for dual SoC and remaining charge capacity estimation.

#### *5.2. SPKF-Based Capacity Estimator*

Since the battery capacity model output equation is distinctively nonlinear, the EKF-based estimator application has been found to give too noisy estimates with slow convergence rate. This is an expected result since EKF uses analytic linearization through Taylor series expansion around the current operating point, i.e., around the state variable (in this case capacity *C*) mean value. Another, more coherent approach to this problem is statistical linearization which linearizes the model at multiple points drawn from prior distribution of *C*. The estimator derived using this approach is called SPKF [3]. There is a couple of SPKF versions which differ in calculation of sigma-points for

linearization; in this paper, the method called central-difference Kalman filter (CDKF) is used because it provides simple parametrization without compromising accuracy [3]. Comparison between EKFand SPKF-based capacity estimation, shown in Figure 9, clearly illustrates the benefits of using SPKF when compared to EKF.

**Figure 9.** Comparison between EKF- and SPKF-based capacity estimation, where the estimated capacity is not fed back to DEKF-based state and parameter estimator.

It is important to note that in the case shown in Figure 9 the capacity estimates were not fed back into the state model of the DEKF, i.e., into Equation (6). If this were the case, i.e., if the state model of DEKF was updated with capacity estimates every *L* time stamps, the estimator would not converge to correct estimates, as shown in Figure 10a–c. This is because every model parameter is estimated in a coupled manner, so there are multiple parameter combinations where output voltage residual would be minimized. For instance, Figure 10d shows an estimate of *Uoc*(*SoC*) which is narrower than the actual curve, because the capacity is estimated higher than the actual one.

**Figure 10.** SPKF-based capacity estimation with capacity adaptation of the DEKF from the start, i.e., *tstart* = 0: (**a**–**c**) estimated capacity vs. time with zooms, (**d**) estimated and measured *Uoc*(*SoC*) curves.

The capacity estimate feedback to the DEKF should be, therefore, turned on with some delay, i.e., until capacity estimate convergence is detected. For that purpose, capacity convergence detection algorithm has been designed, as presented in the next subsection.

#### *5.3. Capacity Convergence Detection Algorithm*

The capacity convergence detection algorithm is based on monitoring of the normalized estimation error (NEE) [28]:

$$\varepsilon\_{\hat{y}}(k) = (\hat{y}(k) - \hat{y}(k))P\_{\hat{y}}^{-1}(\hat{y}(k) - \hat{y}(k))^T \tag{13}$$

where *y*(*k*) is the SoC estimate generated by the DEKF, *y*ˆ(*k*) is the SoC calculated from the SPKF model output, and *Py* is the SPKF innovation matrix (which is regularly calculated as a part of SPKF; note that it is a scalar in the particular case of single estimated parameter—the capacity). The convergence algorithm monitors the NEE, and when it is lower than some predefined value during some predefined number of consecutive time steps, the convergence is claimed.

#### *5.4. Capacity Estimation Results*

Results of SPKF-based capacity estimation algorithm with delayed and automatically calculated (through capacity convergence detection algorithm) start of capacity update (i.e., *tstart*) within the DEKF state model (version with hysteresis model included was used) are shown in Figure 11. The capacity estimates plotted versus time are shown in Figure 11a along with the "measured" capacity (see Figure 3b for details about capacity identification). Capacity convergence has automatically been detected after 2.9 h and from that point on, SPKF has been coupled to the DEKF. Figure 11b shows capacity estimates during the discharge periods plotted versus SoC and color-mapped with respect to temperature. Capacity shows expected (based on the [29]) positive correlation with the temperature.

**Figure 11.** SPKF-based capacity estimation with automatic convergence detection: (**a**) Capacity estimates vs. time, (**b**) capacity estimates vs. SoC and temperature.

Figure 12 shows the same plots as in the case of Figure 6, but instead of using the adaptation mechanism the estimator relies on the explicit hysteresis model and has the capacity estimation included. The voltage residual is shown in Figure 12a together with the usual statistics. This residual has higher *Lstat* value than the one from Figure 6a, which may be explained by the influence of added capacity estimation. The estimates of *Rohm* and *Rp*, plotted in Figure 12b,d with respect to SoC and temperature, respectively, are similar to those from Figure 6b,d, but with slightly higher correlation with temperature for both resistances. Finally, it should be noted that there are no distinguishable sets of estimated *Uoc* curves in Figure 12c (unlike in Figure 6c), because estimated *Uoc*(*SoC*) now describes the central curve while the hysteresis is accounted for in the model (see Figure 7 and Equations (10)

and (11)). Estimated *Uoc*(*SoC*) is slightly larger than the recorded one (see Figure 3c for details about *Uoc*(*SoC*) identification) because the latter is discharge *Uoc*(*SoC*) curve while we estimate the average *Uoc*(*SoC*) since hysteresis is explicitly modelled in this case. The overall estimation algorithm is parametrized as given in Appendix A.

**Figure 12.** DEKF verification results with added hysteresis model and capacity estimation: (**a**) Voltage residual histogram including normal distribution fit, (**b**,**d**) estimated resistances *Rohm* and *Rp*, (**c**) estimated and recorded *Uoc*(*SoC*) curves.

#### **6. Conclusions**

An algorithm for dual estimation of battery state-of-charge (SoC) and remaining charge capacity has been proposed, which is aimed to be accurate over the whole battery lifetime and real-driving conditions including varying ambient temperatures. This was achieved by simultaneous estimation of relevant battery degradation-dependent parameters such as internal resistances and parameters of open-circuit voltage vs. SoC characteristic, *Uoc*(*SoC*).

To this end, the dual extended Kalman filter-based SoC estimation algorithm has been extended to estimate parameters of the characteristic *Uoc*(*SoC*) along with the resistance parameters. This extension allows the DEKF to adapt for *Uoc*(*SoC*) variations and capture its hysteresis without explicitly modelling it. The latter can be useful in cases when the exact hysteresis profile is not known in advance or when it needs to be updated at the given state-of-health level without a specific identification experiment.

Next, a battery capacity estimator has been designed as a separate estimator, as it is based on a different model than the one that has been used in the DEKF design. Moreover, capacity estimation is meant to be executed on a significantly slower time scale than the DEKF. It has been shown that the EKF-based capacity estimator gives rather inconsistent estimates with a slow convergence rate, which is explained by a distinctively nonlinear capacity model. Capacity estimator has, therefore, been designed by using a sigma-point Kalman filter (SPKF). Furthermore, it has been demonstrated that SoC and capacity estimators (i.e., DEKF and SPKF, respectively) cannot be started in a coupled manner, unless it is ensured that both estimators have converged. A capacity convergence detection algorithm has, therefore, been designed to automatically couple the two estimators.

Finally, the overall estimator has been successfully verified based on real driving cycle data acquired by using a fully electric scooter equipped with a telemetry measurement system. The DEKF output voltage estimation residual distribution was confirmed to be close to the voltage measurement noise, while resistance estimates showed expected correlations with temperature. The estimated capacity was shown to be close to the measured one and expectedly correlated with temperature, as well.

Future work will be directed towards further extensions and verifications of the proposed estimator to account for temperature- and aging-dependent variations of the polarization time constant τ*<sup>p</sup>* and further analyze the sensitivity of estimator for broader operating conditions (e.g., wider temperature range), respectively. The emphasis will be on using the estimator to track battery degradation features in support of modelling the battery degradation process.

**Author Contributions:** Conceptualization, J.D., F.M., and M.H.; Methodology, F.M., M.H., and J.D.; Software, F.M.; Validation, F.M., J.D., and M.H.; Writing—Original Draft Preparation, F.M.; Writing—Review and Editing, J.D. and M.H.; Supervision: J.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** It is gratefully acknowledged that this work has been supported through a small-scale grant program of the University of Zagreb, related to the project entitled "Experimental Characterization of Electric Scooter Battery Pack Aging", as well as through project "Advanced Methods and Technologies in Data Science and Cooperative Systems" (DATACROSS) of the Centre of Research Excellence within the University of Zagreb.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **Appendix A. Estimator Parameters**

The overall estimation algorithm is parametrized as given in Table A1.



#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **SOC and SOH Joint Estimation of the Power Batteries Based on Fuzzy Unscented Kalman Filtering Algorithm**

#### **Miaomiao Zeng, Peng Zhang, Yang Yang, Changjun Xie \* and Ying Shi \***

School of Automation, Wuhan University of Technology, Wuhan 430070, China **\*** Correspondence: jackxie@whut.edu.cn (C.X.); yingsh@whut.edu.cn (Y.S.)

Received: 5 July 2019; Accepted: 12 August 2019; Published: 14 August 2019

**Abstract:** In order to improve the convergence time and stabilization accuracy of the real-time state estimation of the power batteries for electric vehicles, a fuzzy unscented Kalman filtering algorithm (F-UKF) of a new type is proposed in this paper, with an improved second-order resistor-capacitor (RC) equivalent circuit model established and an online parameter identification used by Bayes. Ohmic resistance is treated as a battery state of health (SOH) characteristic parameter, F-UKF algorithms are used for the joint estimation of battery state of charge (SOC) and SOH. The experimental data obtained from the ITS5300-based battery test platform are adopted for the simulation verification under discharge conditions with constant-current pulses and urban dynamometer driving schedule (UDDS) conditions in the MATLAB environment. The experimental results show that the F-UKF algorithm is insensitive to the initial value of the SOC under discharge conditions with constant-current pulses, and the SOC and SOH estimation accuracy under UDDS conditions reaches 1.76% and 1.61%, respectively, with the corresponding convergence time of 120 and 140 s, which proves the superiority of the joint estimation algorithm.

**Keywords:** power batteries; improved second-order RC equivalent circuit; fuzzy unscented Kalman filtering algorithm; joint estimation

#### **1. Introduction**

The power batteries serving as the power supply for electric vehicles (EVs) have direct effects on the overall performance of EVs, and the battery overcharge may cause overheating or even an explosion, while the battery over-discharge may result in accelerated aging and permanently reduced capacity [1]. Concerning the issues of safety usage, the state estimation of batteries available for safety precautions can facilitate the elimination of safety hazards, which means the state of charge (SOC) and state of health (SOH) joint estimation is of great significance for the research on power batteries [2].

Lots of scholars have proposed many SOC estimation methods, such as the open circuit voltage method [3,4], the Coulomb counting method [5], the neural network method [6] and the Kalman filtering algorithm [7]. Among them, the open circuit voltage method was to first establish a corresponding function of the open circuit voltage and the SOC and then obtain the SOC by measuring the open circuit voltage after the battery was stationary [8]; the Coulomb integral method, which discretizes the current flowing through the battery and sums it up, and obtains the SOC value by simple division [9]; the neural network method optimizes the relevant parameters of the SOC estimation algorithm and solves complex abstract problems through autonomous learning [7]; a series of Kalman filtering algorithms based on the extended Kalman filtering algorithm optimize autoregressive data processing, which can make the optimal estimation in the minimum variance sense for the state of the dynamic system [10,11].

The estimation methods in the SOH are mainly divided into two categories: One is to start with the characteristic parameters of the battery, and the other is to analyze the aging characteristics and electrochemical reaction characteristics of the battery [12]. The former mainly uses the direct measurement method, obtaining the current SOH by obtaining aging characteristic parameters such as capacity and ohmic internal resistance [8,13]. There are also methods such as neural networks [14] and fuzzy logic [15], which can directly estimate the SOH of the battery through data training without an a priori model. The latter uses an electrochemical model method [12] that models the internal physical and chemical reactions during the charging process and designs an estimator for SOH estimation. There is also a method based on an equivalent circuit model [16] that establishes a circuit that reflects internal variables for SOH estimation.

All of the above algorithms are only a single estimate for the SOC or the SOH, ignoring the close relationship between the SOC and the SOH. The SOC estimate is affected by battery aging—as the battery ages, inaccurate SOC estimates can affect the SOH correction. Therefore, a joint estimate of the SOC and the SOH is necessary. The literature proposes an online SOH estimation method for the lithium battery using the constant-voltage (CV) charge current, as proposed in reference document [17], which can ensure the estimation error of less than 2.5%. However, it is difficult to accurately estimate the true state of the lithium battery by merely estimating the value of the SOH. Another SOC and SOH joint estimation method applicable to the cycle life of lithium-ion batteries for EVs, as proposed in reference document [18], involves an SOC and SOH identification using offline state estimators with different time scales; this requires substantial data to ensure asymptotic convergence without the real-time update.

Reference document [19] analyzed the error sources from the four angles of measurement, model, algorithm and state parameters for the SOC estimation. Finally, the author put forward new concerns in the practical application of SOC estimation. A multi-time-scale observer of the SOC and the SOH for a lithium-ion battery with coupled fast and slow dynamics was proposed in the reference document [16]. The authors used a deterministic transformation of the extended Kalman filter. The paper made an effective estimation of the SOC and the SOH by strictly characterizing the stability of estimation error. Three model-based filtering algorithms [20] were used to estimate the SOC, and the tracking accuracy, calculation time, robustness, etc., were analyzed and compared. Experimental results showed the advantages of three algorithms; the unscented Kalman filtering (UKF) algorithm has a good stability and the Particle filter (PF) algorithm, in the early stage has extreme rapidity. This article gave a combination of the two algorithms to improve the accuracy of the research direction.

In this paper, full consideration was given to the estimation error caused by the change in ohmic resistance during the service of power batteries, and the constant ohmic resistance was replaced by that of gentle variations resistance so as to propose a joint estimation algorithm of the power battery SOC and SOH based on a fuzzy control trace-free Kalman filter. This algorithm uses two complete fuzzy unscented Kalman filtering (F-UKF) algorithms to estimate the SOC and ohmic resistance of the battery at the same time. First of all, the use of a fuzzy controller can effectively reduce the impact of observation noise under complex conditions and to further improve the accuracy of battery SOC estimation. Secondly, the fuzzy controller is used to make a real-time correction of the variance matrix of the observed noise so as to finally realize the estimation of the battery ohmic internal resistance; experiments show that the joint estimation algorithm is not affected by the initial value of SOC, and it still has good convergence speed and tracking accuracy under complex conditions.

The rest of this paper is organized as follows: Section 2 introduces the model of lithium battery, open circuit voltage, SOC calibration experiment, and parameter identification. Section 3 reviews the implementation method of traceless Kalman filtering, fully considers the intrinsic coupling relationship between the SOC and the SOH, puts forward the fuzzy and traceless Kalman filter algorithm on the basis of traceless Kalman filtering, and uses two F-UKF algorithms to estimate the SOC and ohmic internal resistance at the same time. Section 4 discusses the relevant experimental process and conclusions, and Section 5 summarizes the full text.

#### **2. Model for the Lithium Battery**

#### *2.1. Setup of Equivalent Circuit Model for the Lithium Battery*

An accurate battery model can effectively describe the external features and characteristics of internal electrochemical reactions, which is of great significance for the SOC and SOH evaluation of power batteries [21]. In this paper, the SOC is defined as the ratio between remaining battery capacity and nominal battery capacity under the same environmental conditions and specified discharge rate [22]:

$$\text{SOC} = \frac{Q\_{\text{res}}}{Q\_N} \times 100\% \tag{1}$$

In which *Qres* is the remaining battery capacity after the discharge of partial electric quantity and *QN* is the nominal battery capacity.

Through the comparative analysis of differences between old and new batteries, the researchers found that the ohmic resistance and actual maximum battery capacity have more significant changes due to the SOH variations, and SOH is defined as follows from the perspective of ohmic resistance [23]:

$$SOH\_R = \frac{R\_0(end) - R\_0(t)}{R\_0(end) - R\_0(0)} \times 100\% \tag{2}$$

In which *SOHR* is the battery SOH, which defined based on the ohmic resistance *R*0; *R*0(*end*) is the ohmic resistance when the actual maximum battery capacity drops to 80% of the nominal battery capacity; *R*0(*t*) is the ohmic resistance of the battery at *t*; and *R*0(0) is the ohmic resistance upon the battery delivery from the factory.

The proposed improved second-order Resistor-capacitor (RC) equivalent circuit model based on the equivalent circuit model [24,25] is shown in Figure 1. The high capacitance *Cp* and current-controlled current source (CCCS) on the left side characterize the battery capacity, SOC and running time. The second-order RC circuit on the right side simulates the internal polarization characteristics of the battery, R1 and *C*<sup>1</sup> describe the concentration polarization characteristics of the battery, while *R*<sup>2</sup> and *C*<sup>2</sup> describe the electrochemical polarization characteristics of the battery. The voltage-controlled voltage source (VCVS) simulates the nonlinear relationship between the open-circuit voltage and *Usoc*, which links the circuit parts on both sides.

**Figure 1.** Improved second-order resistor-capacitor (RC) equivalent circuit model.

Based on the improved second-order RC equivalent circuit model for the lithium battery, select *x* = [SOC *U*<sup>1</sup> *U*2] *<sup>T</sup>* as the state variable to obtain the following continuous state space equation:

$$
\begin{bmatrix}
\dot{\mathbf{SOC}} \\
\dot{\boldsymbol{\mathcal{U}}}\_{1} \\
\dot{\boldsymbol{\mathcal{U}}}\_{2}
\end{bmatrix} = \begin{bmatrix}
0 & 0 & 0 \\
0 & -\frac{1}{\overline{\mathcal{R}\_{1}\mathbf{C}\_{1}}} & 0 \\
0 & 0 & -\frac{1}{\overline{\mathcal{R}\_{2}\mathbf{C}\_{2}}}
\end{bmatrix} \cdot \begin{bmatrix}
\text{SOC} \\
\boldsymbol{\mathcal{U}}\_{1} \\
\boldsymbol{\mathcal{U}}\_{2}
\end{bmatrix} + \begin{bmatrix}
\end{bmatrix} \cdot \boldsymbol{i} \tag{3}
$$

In Equation (3), *QN* is the nominal battery capacity and η is the charge–discharge efficiency of the battery. The discretized state equation and observation equation are as follows:

$$\begin{cases} \mathbf{x}(k+1) = A \cdot \mathbf{x}(k) + B \cdot \dot{\mathbf{z}}(k) \\ \mathcal{U}\_L(k) = \mathcal{U}\_{0\mathbf{c}}(\text{SOC}) - \mathcal{U}\_1(k) - \mathcal{U}\_2(k) - R\_0 \cdot \dot{\mathbf{z}}(k) \end{cases} \tag{4}$$

In which *T* is the sampling period of the system, with *A* and *B* expressed as follows:

$$A = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \exp(-\frac{T}{R\_1 C\_1}) & 0 \\ 0 & 0 & \exp(-\frac{T}{R\_2 C\_2}) \end{bmatrix}, B = \begin{bmatrix} -\frac{\eta \cdot T}{Q \circ i} \\ R\_1 (1 - \exp(-\frac{T}{R\_1 C\_1})) \\ R\_2 (1 - \exp(-\frac{T}{R\_2 C\_2})) \end{bmatrix} \tag{5}$$

#### *2.2. Open-Circuit Voltage and SOC Setting Experiments*

The procedures of open-circuit voltage and SOC setting experiments for the lithium battery at a normal temperature (25 ◦C) based on the ITS5300 battery test platform (ITECH ELECTRONIC CO., LTD., Nanjing, China) are as follows: Charge the battery until the full-load capacity is reached before the 3-hour standing and record the open-circuit voltage of the battery, discharge the battery for 6 minutes at a discharge rate of 1 C (40 A), and repeat the above steps until the cutoff voltage is reached. The fitting of open-circuit voltage curve corresponding to the SOC variations of lithium battery was completed via the MATLAB software (2017a, The MathWorks, Inc, Natick, MA, USA), which showed that the fitting curve had the minimum root-mean-square error when the polynomial order was 5. The fitting curve is shown in Figure 2, and the function expression is as follows.

#### *Uoc*(*SOC*) = 3.2821 · *SOC*<sup>5</sup> <sup>−</sup> 10.3004 · *SOC*<sup>4</sup> <sup>+</sup> 13.0068 · *SOC*<sup>3</sup> <sup>−</sup> 7.9724 · *SOC*<sup>2</sup> <sup>+</sup> 2.4054 · *SOC* <sup>+</sup> 2.9752 (6)

**Figure 2.** Fitting curves of the open-circuit voltage and state of charge.

#### *2.3. Parameter Identification of the Lithium Battery Model*

In accordance with the improved second-order RC equivalent circuit model, the Bayesian identification algorithm based on the least-square equation was adopted for the identification of resistance and capacitance parameters of the equivalent circuit model. Taking the parameters to be estimated as random variables, the Bayesian identification algorithm achieved the optimal estimation indirectly through the observation on other related parameters [26,27].

Kirchhoff's law should be adopted to obtain the following Laplace's equation of the improved second-order RC equivalent circuit model:

$$\text{lL}\_{\text{oc}}(s) - \text{lL}\_{\text{L}}(s) = i(s) \cdot \left( \frac{R\_1}{R\_1 C\_1 s + 1} + \frac{R\_2}{R\_2 C\_2 s + 1} + R\_0 \right) \tag{7}$$

The equation obtained using the bilinear transformation method is as follows:

$$d(k) = -k\operatorname{id}(k-1) - k\operatorname{zd}(k-2) + k\operatorname{si}(k) + k\operatorname{si}(k-1) + k\operatorname{si}(k-2) \tag{8}$$

In which *i*(*k*) is the system input, *d*(*k*) = *UOC*(*k*) − *UL*(*k*), and the final derivation of the Bayesian identification algorithm is as follows.

$$\begin{cases} \begin{aligned} \boldsymbol{\theta}(k) &= \boldsymbol{\theta}(k-1) + \boldsymbol{K}(k) \cdot \left[ \boldsymbol{z}(k) - \boldsymbol{H}^{T}(k) \boldsymbol{\theta}(k-1) \right] \\ \boldsymbol{K}(k) &= \boldsymbol{P}\_{\boldsymbol{\theta}}(k-1) \boldsymbol{H}^{T}(k) \cdot \left[ \boldsymbol{H}^{T}(k) \boldsymbol{P}\_{\boldsymbol{\theta}}(k-1) \boldsymbol{H}(k) + \frac{1}{\sigma\_{r}^{2}} \right]^{-1} \\ \boldsymbol{P}\_{\boldsymbol{\theta}}(k) &= \left[ \boldsymbol{I} - \boldsymbol{K}(k) \boldsymbol{H}^{T}(k) \right] \cdot \boldsymbol{P}\_{\boldsymbol{\theta}}(k-1) \end{aligned} \tag{9}$$

The initial value of θ(0) is 0, and the initial value of covariance matrix *P*<sup>θ</sup> (0) is *a*·*I*, among which *a* is a small positive number and *I* is a 5-order unit matrix. Use the recursion Formula (9) of the Bayesian identification algorithm to estimate the model parameters and then calculate the resistance and capacitance values of the model via Equation (8). In practical applications, it is necessary to consider the amount of calculation and the length of time for parameter identification. The joint estimation algorithm designed in this paper has a large amount of computation. Therefore, the mean value of the online identification result was selected as the parameter identification result. The results are shown in Table 1.

**Table 1.** Online identification result based on the Bayesian identification algorithm.


#### **3. SOC and SOH Joint Estimation Based on F-UKF**

#### *3.1. Unscented Kalman Filtering Algorithm*

The unscented Kalman filtering algorithm adopts the linear Kalman filter framework instead of the traditional linearization for nonlinear functions, with the nonlinear transfer of mean value and covariance completed via unscented transformation in the one-step prediction equation [7,28]. The unscented Kalman filtering algorithm is applicable to the nonlinear dynamic systems described with the following state-space equation:

$$\begin{cases} \mathbf{x}(k+1) = f[\mathbf{x}(k), \mathbf{u}(k)] + \mathbf{c}(k) \\ \mathbf{y}(k) = \mathbf{g}[\mathbf{x}(k), \mathbf{u}(k)] + \mathbf{v}(k) \end{cases} \tag{10}$$

In which *f* is the function of nonlinear state equation and *g* is the function of nonlinear observation equation. Assume that *e*(*k*) has the covariance matrix *Q* and *v*(*k*) has the covariance matrix *R*; thus, the essential operation steps of the unscented Kalman filtering algorithm for a random variable *X* at the different time *K* are shown in Figure 3.

**Figure 3.** Essential operation steps of the unscented Kalman filtering algorithm.

#### *3.2. Fuzzy Unscented Kalman Filtering Algorithm*

The application of the unscented Kalman filtering algorithm should have been based on the already known statistical characteristics of process noises and observation noises; however, the insufficient estimation accuracy caused by the difficulty in noise determination during the use of power batteries required the introduction of adaptive filtering technique for algorithm optimization [29,30]. With reference to the unscented Kalman filtering algorithm, the covariance matching technique based on the fuzzy inference system was adopted in this section to effectively improve the accuracy of real-time observation noise estimation.

Assume that the statistical characteristics of process noises are already known, implement the recursive correction of observation noise variance based on the calculation of real-time ratio between the theoretical and actual covariances of observation errors.

Calculate the theoretical covariance *N*(*k*) and actual covariance *M*(*k*) of observation errors first, among which *i* = *k* − *n* + 1.

$$N(k) = \sum\_{i=0}^{2n} a\_c^i \cdot \varepsilon\_y(k|k-1) \cdot \varepsilon\_y^{\top}(k|k-1) + V(k) \tag{11}$$

$$M(k) = \frac{1}{n} \sum\_{i}^{k} \left[ y(i) - y(i|i-1) \right] \cdot \left[ y(i) - y(i|i-1) \right]^T \tag{12}$$

The input value *G*(*k*) of the fuzzy controller is a ratio between the theoretical and actual covariances of observation errors, while its output value α(*k*) is the adjustment factor of observation noise variance. Take the adjusted observation noise variance *V*ˆ (*k*) into the unscented Kalman filtering algorithm to calculate the new observation-error covariance matrix *<sup>P</sup>*ˆ*y*(*k*|*<sup>k</sup>* <sup>−</sup> 1) and then update the Kalman filter gain and state-error covariance matrix with the updated *<sup>P</sup>*ˆ*y*(*k*|*<sup>k</sup>* <sup>−</sup> 1).

$$G(k) = \frac{M(k)}{N(k)}\tag{13}$$

$$
\hat{V}(k) = a(k) \cdot V(k) \tag{14}
$$

As a kind of uncertainty reasoning method, the fuzzy controller is composed of three parts. First, initiate the fuzzy processing in accordance with the input membership function shown in Figure 4 for the input value *G*(*k*) of the fuzzy controller to obtain the corresponding fuzzy index.

**Figure 4.** Input membership function.

Second, concerning the fuzzy controller with a single input/output, the correspondent fuzzy rules are relatively simple. The greater observation noise will result in the greater actual covariance *M*(*k*) and *G*(*k*), while the change in theoretical covariance *N*(*k*) is subject to the variation of observation noise variance *V*(*k*). In order to keep the variation consistency between *N*(*k*) and *M*(*k*), adjust α(*k*) to enlarge *V*(*k*) when the observation noise becomes greater, so that the decreased *G*(*k*) will approach 1. The decreased observation noise will result in the decreased actual covariance *M*(*k*) and *G*(*k*). Therefore, α(*k*) should be adjusted accordingly to decrease *V*(*k*) so that the enlarged *G*(*k*) will approach 1. The fuzzy rules established in accordance with the above derivation process is shown in Table 2.



Finally, initiate the anti-fuzzy processing in accordance with the output membership function shown in Figure 5 for the output fuzziness to obtain the output value α(*k*) of the fuzzy controller.

**Figure 5.** Output membership function.

The diagram of working principle using the fuzzy controller for the adjustment of observation noise variance matrix is shown in Figure 6.

**Figure 6.** Structure diagram of the fuzzy controller.

#### *3.3. Design and Implementation of the SOC and SOH Joint Estimation Algorithm*

The following state-space equation can be used to express the improved second-order RC equivalent circuit model:

$$\begin{cases} x(k+1) = A\_x(k) \cdot x(k) + B\_x(k) \cdot i(k) + \varepsilon\_x(k) \\ y(k) = \mathcal{U}\_{\text{oc}}(\text{SOC}(k)) - \mathcal{U}\_1(k) - \mathcal{U}\_2(k) - \mathcal{R}\_0 \cdot i(k) + v\_x(k) \end{cases} \tag{15}$$

In which *x*(*k*) is the state variable; *y*(*k*) is the predicted terminal voltage of the battery; *ex*(*k*) is the process noise, with the mean value of zero and variance *Ex*(*k*); *vx*(*k*) is the observation noise, with the mean value of zero and variance *Vx*(*k*); and *Ex*(*k*) and *Vx*(*k*) are irrelevant.

The gradual increase of ohmic resistance in a non-linear way is unnoticeable within a short period, which means the ohmic resistance of the battery at two adjacent moments can be taken as the constant value. Therefore, the state equation and observation equation for the estimation of ohmic resistance can be expressed as follows:

$$\begin{cases} \mathcal{R}\_0(k+1) = R\_0(k) + \varepsilon\_R(k) \\ \mathcal{Y}(k) = \mathcal{U}\_{\rm oc}(\text{SOC}(k)) - \mathcal{U}\_1(k) - \mathcal{U}\_2(k) - R\_0(k) \cdot i(k) + \upsilon\_R(k) \end{cases} \tag{16}$$

In which *R*0(*k*) is the state variable; *y*(*k*) is the predicted terminal voltage of the battery; *eR*(*k*) is the process noise, with the mean value of zero and variance *ER*(*k*); *VR*(*k*) is the observation noise, with the mean value of zero and variance *VR*(*k*); and *ER*(*k*)and *VR*(*k*) are irrelevant.

The optimal estimated value of the battery SOC was adopted in the joint estimation algorithm for the one-step-ahead prediction of ohmic resistance; meanwhile, the optimal estimated value of ohmic resistance was also available for the one-step-ahead prediction of the battery SOC, and the above mutual application facilitated the acquisition of the estimated battery SOC and ohmic resistance closer to the actual values.

The flowchart of the SOC and SOH joint estimation algorithm based on F-UKF is shown in Figure 7.

**Figure 7.** Flowchart of the joint estimation algorithm.


#### **4. Experimental Verification and Result Analysis**

The ITS5300-based battery test platform available to verify the proposed SOC and SOH joint estimation algorithm is shown in Figure 8. The nominal capacity of a single lithium iron phosphate battery is 40 Ah, and the corresponding performance parameters are shown in Table 3. In order to measure the terminal voltage and working current of the battery, the software of IT9320 battery test system (ITECH ELECTRONIC CO., LTD, Nanjing, China) was used to simulate the discharge conditions with constant-current pulses and the urban dynamometer driving schedule (UDDS) driving cycles, and the MATLAB software was adopted for the simulation verification and analysis of the joint estimation algorithm proposed in this paper.

**Figure 8.** Battery test platform.


**Table 3.** Parameters of the lithium iron phosphate battery.

#### *4.1. Sensitivity Verification of the F-UKF Algorithm against Initial Values*

Concerning the estimation of the battery SOC and ohmic resistance using the F-UKF algorithm, it was difficult to obtain the accurate initial values of battery SOC, but the values of ohmic resistance were relatively stable without violent fluctuations. The lithium iron phosphate battery was charged until the battery SOC reached 85% of the initial state before the experiment. Under the discharge conditions with constant-current pulses, the different initial values of the battery SOC were set to verify the F-UKF sensitivity against initial values. In the MATLAB software, the respective initial values of the SOC were set to 40%/0% and 85% for the F-UKF algorithm and the ampere-hour integration

method, with the sampling period and discharge rate set to 1 s and 0.5 C (20 A) for the 500-second simulation experiment.

The simulations under discharge conditions with constant-current pulses shown in Figures 9 and 10 indicate that the different initial values of battery SOC converged to the vicinity of reference value after a period of filtering iteration. Though the greater deviation of SOC initial values resulted in a longer convergence time, the stabilized values could follow the reference value well, and the estimation error was extremely small. Therefore, the F-UKF algorithm proposed in this paper is insensitive to the initial values.

**Figure 9.** State of charge (SOC) estimation curve based on the fuzzy unscented Kalman filtering algorithm (F-UKF) under discharge conditions with constant-current pulses.

**Figure 10.** SOC estimation error curve based on F-UKF under discharge conditions with constant-current pulses.

#### *4.2. Joint Simulation Verification of UDDS Driving Cycles*

The lithium iron phosphate battery was charged until the battery SOC reached 85% of the initial state before the experiment.

The software of IT9320 battery test system was used to compile the pulse current driving cycles before the acquisition of terminal voltage and working current from UDDS driving cycles. In the MATLAB software, the initial values of SOC were set to 80% for the F-UKF and joint estimation algorithms and 85% for the ampere-hour integration method, respectively. The initial value of ohmic resistance was set to 1.50 mΩ for the joint estimation algorithm, which was greater than the reference value of 1.278 mΩ.

The simulations of UDDS driving cycles shown in Figure 11 indicate that the convergence time of SOC from the initial 80% to the vicinity of reference value based on the F-UKF algorithm was about 170 s, while the UKF algorithm was about 185 s. The SOC estimation error of the F-UKF algorithm after convergence could be controlled within 2.82%, while the estimation error of the UKF algorithm was about 2.93%. Since the noise of UKF was random white noise, the F-UKF algorithm that introduces adaptive technology was to adjust the noise instead of eliminating the noise. It can be seen that the convergence performance of the F-UKF algorithm was not only better than the UKF algorithm, but the estimation accuracy was also relatively improved in complex conditions.

**Figure 11.** (**a**) SOC estimation curve of urban dynamometer driving schedule (UDDS) driving cycles based on the F-UKF and UKF algorithms. (**b**) SOC estimation error curve of UDDS driving cycles based on the F-UKF and UKF algorithms.

The simulations of UDDS driving cycles shown in Figure 12a indicate that the convergence time of SOC from the initial 80% to the vicinity of reference value based on the F-UKF algorithm was about 170 s, while the corresponding convergence time with an increased rising velocity based on the joint estimation algorithm was about 120 s. Therefore, the convergence performance of the joint estimation algorithm was better than that of the F-UKF algorithm under complex conditions. Figure 12b shows that the respective SOC estimation errors of the F-UKF algorithm and the joint estimation algorithm after convergence were less than 2.82% and 1.76%. Therefore, the tracking performance of the joint estimation algorithm was better than that of the F-UKF algorithm in terms of the SOC estimation.

**Figure 12.** (**a**) SOC estimation curve of UDDS driving cycles based on the F-UKF and joint estimation algorithms. (**b**) SOC estimation error curve of UDDS driving cycles based on the F-UKF and joint estimation algorithms.

The simulations of UDDS driving cycles shown in Figures 13 and 14 indicate that the convergence time of ohmic resistance from 1.50 mΩ to the vicinity of reference value (1.278 mΩ) based on the joint estimation algorithm was about 140 s, and the corresponding battery SOH was about 89.87% of the reference value after stabilization. Figure 15 shows that the maximum SOH estimation error based on the joint estimation algorithm was 1.61%, and the SOH estimation error was less than 1.20% over time.

**Figure 13.** Ohmic resistance estimation curve of UDDS driving cycles based on the joint estimation algorithm.

**Figure 14.** State of health (SOH) estimation curve of UDDS driving cycles based on the joint estimation algorithm.

**Figure 15.** SOH estimation error curve of UDDS driving cycles based on the joint estimation algorithm.

#### **5. Conclusions**

In order to implement the real-time state estimation of power batteries for EVs, taking account observation noises and gradually changed ohmic resistance, an improved second-order RC equivalent circuit model was established in this paper for the SOC and SOH joint estimation using the fuzzy unscented Kalman filtering algorithm (F-UKF). The experimental data obtained from a test bench was adopted for the simulation to verify the convergence and stability of the F-UKF algorithm and to achieve the required design effects. The experimental data were obtained by the ITS5300 battery test platform, and the proposed joint estimation algorithm considered the influence of ohmic internal resistance change and noise error.

**Author Contributions:** Conceptualization, Y.S.; data curation, M.Z.; formal analysis, Y.Y.; funding acquisition, C.X.; investigation, M.Z.; methodology, P.Z.; project administration, Y.Y.; software, P.Z.; supervision, C.X.; validation, Y.S.; writing—original draft, M.Z.; writing—review & editing, C.X. and Y.S.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 51477125", "the Hubei Science Fund for Distinguished Young Scholars, grant number 2017CFA049" and "the Hubei province technological innovation major project, grant number 2018AAA059".

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

### **High Reynold's Number Turbulent Model for Micro-Channel Cold Plate Using Reverse Engineering Approach for Water-Cooled Battery in Electric Vehicles**

#### **Satyam Panchal 1,2,\*, Krishna Gudlanarva 2, Manh-Kien Tran 3, Roydon Fraser <sup>1</sup> and Michael Fowler <sup>3</sup>**


Received: 4 March 2020; Accepted: 30 March 2020; Published: 2 April 2020

**Abstract:** The investigation and improvement of the cooling process of lithium-ion batteries (LIBs) used in battery electric vehicles (BEVs) and hybrid electric vehicles (HEVs) are required in order to achieve better performance and longer lifespan. In this manuscript, the temperature and velocity profiles of cooling plates used to cool down the large prismatic Graphite/LiFePO4 battery are presented using both laboratory testing and modeling techniques. Computed tomography (CT) scanning was utilized for the cooling plate, Detroit Engineering Products (DEP) MeshWorks 8.0 was used for meshing of the cooling plate, and STAR CCM+ was used for simulation. The numerical investigation was conducted for higher C-rates of 3C and 4C with different ambient temperatures. For the experimental work, three heat flux sensors were attached to the battery surface. Water was used as a coolant inside the cooling plate to cool down the battery. The mass flow rate at each channel was 0.000277677 kg/s. The k-ε model was then utilized to simulate the turbulent behaviour of the fluid in the cooling plate, and the thermal behaviour under constant current (CC) discharge was studied and validated with the experimental data. This study provides insight into thermal and flow characteristics of the coolant inside a cooing plate, which can be used for designing more efficient cooling plates.

**Keywords:** heat and mass transfer; thermal analysis; Lithium-ion battery; micro-channel cooling plate; battery thermal management; MeshWorks; CFD

#### **1. Introduction**

The collective effects of global warming, environmental degradation, and energy crisis have prompted attention towards clean and sustainable energy [1]. However, there is inconsistency of renewable energy harvesting, since it depends on the effects of climate, which could result in complications in providing sufficient electricity in contrast with traditional nonrenewable energy sources. This has led to an interest in developing large-scale energy storage systems (ESS), for which batteries show promise [2]. Among the available secondary batteries, lithium-ion and lead-acid are broadly considered as effective candidates for energy storage systems. In the automotive industry, plug-in hybrid electric vehicles (PHEVs), hybrid electric vehicles (HEVs), and battery electric vehicles (BEVs) commonly utilize lithium-ion batteries (LIBs) [3]. The widespread use of LIBs is the result of 1) high specific power and energy densities [4]; 2) high nominal voltage and low self-discharge rate [5];

and 3) long cycle-life and no memory influence [6,7]. All these characteristics are required for electric vehicles (EVs) to achieve desirable driving range and vehicle speed [8]. In addition to the driving range and vehicle speed, the life cycle of the LIBs is also a critical factor in EVs. Some important factors in determining the allowable discharging and charging currents (also known as C-rates) and the batteries' life cycle are battery materials, working temperature, and assembling process. Current research on enhancing the life cycle of a battery has mainly been focused on the improvement in the materials and assembling technology, with a specific goal to obtain desired energy density. In addition, little attention has been directed toward the advancement and change of battery cooling systems (BCS) [9].

Pouch-type LIBs are being commonly utilized in new EVs. However, many problems in their safety and life span remain. First, when using the pouch-type lithium-ion battery, particularly in cases of high discharge rate, high heat generation may occur [10]. This type of battery expands due to overheating when the heat is not removed promptly. In some cases, the LIB may even burst and explode. In addition, these pouch-type batteries are connected either in parallel or in series within the LIB packs or modules, which also generate high amount of heat during both discharging and charging. Therefore, a good battery thermal management system (BTMS) is necessary [11]. The heat of LIBs increases when EVs accelerate and experience fast charging. If this generated heat is not adjusted or if it is overtaken by the rate of heat production, the battery pack temperature drastically increases. A high operating temperature of lithium-ion batteries can lead to capacity fade of the battery [12]. The impact of high working temperatures on the execution of a LIB, particularly for cylindrical battery cells (Sony 18650), was researched by Ramadass et al. [13] and the prismatic LIB cell (A123 20 Ah) was explored by Panchal et al. [14]. The authors found that the capacity fading is not the main negative impact related to high working temperatures of the battery; it also may lead to the explosion of the electrolyte. In addition, thermal runaway of a LIB cell can cause the whole LIB pack to fail [15]. In addition, there is a major effect on the electrochemical behaviors in terms of the degradation of electrolyte, electrodes, separator, and the life cycle cost [14,16]. Hence, a robust BTMS is required, to achieve better LIB pack performance in low-temperature conditions and a better lifespan in high-temperature conditions [12,17,18]. A typical operating range is between 20 ◦C and 40 ◦C [19], and an extended range is between −10 ◦C and +50 ◦C for certain applications [20].

Based on the coolant utilized in the BTMS, the BTMS can be divided into (1) air-, (2) fluid- or liquid-, and (3) phase change material (PCM)-based. For the air-based BTMS, cooling of the batteries is done by airflow passing between the batteries in a module or pack. The stream of air can come from the motion of the vehicle, and these types of systems are called natural air-cooled BTMS. If the airflow is generated by power-operated equipment, then these systems are called forced-air BTMS. The air streaming direction for air-based BCSs is usually 90 degrees to the axes for cylindrical batteries. Fluid-based BTMS utilizes coolants in the liquid phase, and such systems typically use power-operated equipment to move the coolant flow. In the PCM cooling method, the cooling is provided by the latent heat of the PCM, while in air- and liquid-cooling methods, the cooling is provided by sensible cooling. The advantages and disadvantages of air, fluid and PCM cooling are: (1) the air-cooling method is simple and lightweight [21]; (2) the water-cooling method is more efficient since it absorbs more heat, and takes less volume, yet more complexities are involved, including high weight and cost [22]; (3) in comparison with the air-cooling strategy, the liquid-cooling strategy provides better cooling performance owing to the high thermal conductivity of water compared to air [23]; and (4) because of the low thermal conductivity of air [24], a high speed of air is required in order to provide adequate cooling for LIBs [25,26]. The advantage of PCM over fluid and air-based strategies is better temperature consistency throughout the battery pack. Most recent examples of EVs and HEVs that utilize air-based cooling systems are the Toyota Prius, the Nissan Leaf, and the Honda Insight [27]. A very noted use of fluid-cooling systems for thermal management of LIBs is in Tesla vehicles, including the Roadster. These cooling techniques utilize a 50/50 mixture of water and glycol to keep the battery pack temperature within the appropriate limits [28]. Both the Chevrolet Bolt and BMW i3 utilize a bottom-cooling plate in their battery packs, with the cooling medium being a water-glycol solution

for the Chevrolet and a refrigerant for the BMW [29]. The objective of this study is to conduct a reverse engineering study to experimentally investigate the design and heat generation of the 20 Ah lithium-ion battery with the cold plate. Subsequently, heat flux data and developed CAD modeling will be used to numerically characterize the thermal and flow behaviors of coolant in the cold plate at different inlet coolant temperatures and higher C-rates. The study provides insight into thermal and flow characteristics of the coolant inside a cooing plate which can be used for future improvements of the cooling plate.

#### **2. State of the Art**

According to the research by several scholars from all over the world, there are numerous papers available on BTMS with air, water, and PCM cooling, and battery modeling in the open literature [28–33].

Al-Hallaj et al. [34] were the first to utilize a PCM experiencing a solid to fluid change, and the authors used the PCM to cool down the 18,650 cylindrical-shaped lithium-ion batteries. In the BCS, they utilized a paraffin blend of pentacosane and hexacosane as the PCM. Moreover, for a passive cooling strategy, PCM-based BTMS is less expensive, requires smaller volume, and accomplishes preferred temperature consistency over air and fluid-based BTMS. The authors additionally utilized a commercial finite-element (FE) software, PDEase2D™, to simulate the thermal behavior of EV battery modules with a PCM BTMS. The authors claimed that the research was important for EV performance under cold conditions, or in space applications where the battery working temperature drops greatly when an orbiting satellite moves from the light to the dark side of the earth.

Zhang et al. [35] worked on the simulation of pouch-type LIB with thermal management using a cooling plate approach. The authors used a 22 Ah battery and developed a computational fluid dynamics (CFD) model using ANSYS Fluent. The authors studied the effect of inlet mass flow rate, difference in temperature, and pressure drop at 4C discharge rate, with the end goal of improving the efficiency and economy of the cooling plate. Their results demonstrated that the increase in mass flow rate of coolant reduced the maximum temperature and temperature difference of cells, but when mass flow rate exceeded 0.003 kg/s, then the economy of the cooling plate worsened.

Omkar et al. [36] developed a PCM/cooling-plate-coupled BTMS using CFD. The authors used LiFePO4/C as a battery and determined the heat generation in simulation. They varied several factors in simulating battery and cold plate, such as the inlet mass flow rate, PCM, thermal conductivity, direction of flow, water cooling, to know the impact on the cooling execution of module. The authors concluded that as the space between adjoining batteries increased, the most extreme temperature had little change, yet the temperature field was uniform.

Chen et al. [37] worked on the comparison of four diverse cooling strategies, including air cooling, indirect fluid cooling, direct fluid cooling, and fin cooling for LIB cells. The authors evaluated the effectiveness on the basis of coolant parasitic power utilization, temperature difference in a cell, the most extreme temperature rise, and extra weight utilized for the cooling strategy. The authors found that (1) to keep the same average temperature, an air cooling strategy needed two to three times more energy than other techniques; (2) an indirect liquid cooling strategy had the lowest maximum temperature rise; (3) a fin cooling system includes an approximately 40% higher weight of cell, and (4) indirect fluid cooling is more efficient than direct fluid cooling.

Lu et al. [38] worked on BTMS of thickly pressed EV batteries with forced-air cooling systems to investigate the air cooling capacity on the temperature consistency and hotspots alleviation of a smaller battery pack subject to different airflow rates as well as airflow paths. The authors used 252 cylindrical lithium-ion batteries (32650). Their numerical outcomes demonstrated that the effective improvement in heat transfer exchange zones between air-coolant and battery surfaces could bring down the highest temperature and improve the most extreme temperature variation in the thickly pressed battery box.

In another study, Qian et al. [39] worked on thermal performance of a thermal management system (TMS) of LIBs by using a minichannel cooling method. The authors utilized a fluid cooling technique based on mini-channel cold plate and later a 3D numerical model was developed. They investigated the effects of the direction of flow, the mass flow rate at the inlet, the impact of number of channels, and the channel width on the thermal performance of the pack. Their outcomes demonstrated that at 5C discharge, the TMS based on minichannel cooling plates was effective in cooling productivity and controlling the battery temperature. They also found that a five-channel cold-plate was sufficient to control the temperature by increasing the mass flow rate at the inlet.

Jarrett et al. [40] developed a CFD model of a battery cooling plate while considering the impact of working conditions on the ideal design of electric vehicle battery cooling plates. The authors considered three important performance measurements: (1) average temperature, (2) temperature consistency, and (3) drop in pressure. They identified that out of these three, temperature consistency was the most sensitive to the working conditions, particularly the circulation of the heat flux input and the flow rate of the coolant.

Zou et al. [41] worked on an experimental study on multiwalled carbon nanotube (MWCNT)-based, graphene-based and MWCNT/graphene-based PCM to enhance the thermal performance of lithium-ion BTMS. Their results demonstrated that a composite PCM mass ratio of 3:7 of the MWCNT/graphene could display the best synergistic improvement for the heat transfer effect, for which the thermal conductivity was increased by 31.8%, 55.4% and 124% compared to graphene-based composite PCM, MWCNT-based composite PCM and pure PCM respectively.

Greco et al. [42] developed a heat-pipe-based BTMS arranged in a sandwiched pattern to improve the cooling for EVs. The authors also built a 1D model utilizing the thermal circuit technique. The proposed model was contrasted to an analytical solution in view of variable partition and CFD simulations in 3D. They found that the higher surface contact of the heat pipes allowed a better cooling management compared to forced convection cooling.

Liang et al. [43] researched the thermal execution of a BTMS under various ambient temperatures using heat pipes. The authors examined impacts of environment temperature, coolant flow rate, coolant temperature, and start-up time on the thermal execution of BTMS. They also claimed that the power utilization can be minimized by diminishing run time of HP-BTMS.

Lastly, Wang et al. [44] experimentally examined a high capacity LiFePO4 battery pack at high temperatures and quick discharge using a novel fluid cooling method. They designed and developed thermal silica plate-based BTMS. Their test results demonstrated that adding the thermal silica plates significantly improved the cooling limit. This can enable the most extreme temperature distinction to be controlled at 6.1 ◦C and decrease the highest temperature by 11.3 ◦C in the battery module.

In the above paragraphs, various methods of lithium-ion BTMS have been studied, and the cooling plates cooling method demonstrates various practical application prospects. The main research contents of the existing literature on the cooling plate cooling methods include the impact of flow direction and mass flow rate. However, as a key part of the cooling system, the effect of external temperature (or boundary conditions) within mini channels was seldom studied. Therefore, in this paper, a cooling plate design and development was done in such a way that it gives maximum cooling close to the anode and the cathode, as the greatest heat production is close to electrodes of LIBs during high acceleration of EVs. A comprehensive examination and modeling was conducted on the cooling plate in the LIB system. From this, the investigation under the rates of 3C and 4C (constant current discharge), and operating conditions of 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C was assessed in detail. In our previous CFD studies (Panchal et al. [45,46]), we designed and developed a cooling plate with only one channel composed of a single inlet and outlet, and put one on both the top and bottom of the battery to cool it down during the discharging rates of 1C, 2C, 3C, and 4C, and diverse cooling temperatures. In this paper, STAR CCM+ was used for CFD simulation and then the simulated results were validated with experimental data of various temperature and velocity profiles. The results of this research

can assist in the design, development and optimization of a cooling-plate cooling system. The data generated is also helpful in battery thermal modeling and EVs development.

The rest of the manuscript is organized as follows. Section 3 introduces the experimental studies including laser scanning, heat flux locations on the surface of the battery, experimental plan and procedure. Section 4 explains the cooling plate physical model in detail, provides geometry and boundary conditions as well as meshing. Section 5 analyzes the results of the numerical calculation, specifically, the effect of inlet mass flow rate on the temperature and velocity profiles of the cooling plate mini channels. Section 6 presents a summary of the conclusions.

#### **3. Experimental Studies**

Here, the lab testing details are given through the laser scanning, test set-up, heat flux distributions, and testing plan and procedure.

#### *3.1. Reverse Engineering*

For reverse engineering, we used a 25S2P battery pack in an EV. The battery pack and reverse engineering are shown in Figure 1.

**Figure 1.** Reverse engineering approach: taking out battery cell from pack.

#### *3.2. Laser Scanning*

The 3D laser-scanning machine used for this work is shown in Figure 2. It consists of probe, Light Emitting Diode (LED) indicators, wide stable joints, dampener, stable base, etc. The laser probe has an accuracy of ±25 μm (±0.001 in); field depth: 115 mm (4.5 in); width of effective scan: near field 3.1 in (80 mm), far field 5.9 in (150 mm); minimum point spacing: 40 μm (±0.0015 in); scan rate: 280 fps (frames/second), 280 fps × 2000, point/line = 560,000 points/sec; and laser class: 2 M. Arm specifications were measuring range: 1.8 m (6 ft); volumetric accuracy: ±0.034 mm (±0.0013 in); single point repeatability: 0.024 mm (0.0009 in); and seven-axis movement. After scanning the microchannel cooling plate, the image was transferred into CAD using DEP MeshWorks 8.0 software.

**Figure 2.** Laser scanning set-up.

#### *3.3. Battery Description, Experimental Set-up and Heat Flux Locations*

The test set-up utilized for this work is explained in detail in our previously published paper [47]. In this work, a different cooling plate configuration was used. Three heat flux sensors were attached on the principle surface of the battery (one near the anode, one near the cathode, and one near the mid body) and the sensor measurements were used for simulation. The heat flux sensor location is shown Figure 3. A pouch-type 20-Ah-capacity LIB cell was utilized for the testing and model validation. Table 1 organizes the LIB cell technical details. The battery cell was placed between two cooling plates to form a sandwich structure.

**Figure 3.** Heat flux sensors locations near (**a**) cathode and (**b**) anode.


**Table 1.** Technical details of 20 Ah lithium-ion battery (LIB) cell.

#### *3.4. Test Plan*

In the experiment, four different fluid inlet temperatures were chosen for the water-cooling strategy: 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C. Two distinctive discharge currents were selected: 60 A and 80 A (3C and 4C). The charge current is 20 A (1C). The testing sequence is presented in Table 2. The test layout and experimental uncertainty are available in our previously published paper [47].


#### **4. Cold Plate Cooling System Modeling**

#### *4.1. Governing Equations*

The fluid stream in this test was considered turbulent because the Reynold's number was 8700. As such, the stream was modeled by Reynolds-averaged Navier–Stokes Equations (RANS). STAR CCM+ software was used for the CFD simulation. In this investigation, the realizable *k-*ε turbulence model was utilized because of the strengths of the model, which include reasonable precision for an extensive variety of flows and its demonstrated capacity in heat transfer and stream examination. The equations used in STAR CCM+ for turbulent kinetic energy and eddy viscosity were:

$$\frac{\partial}{\partial t}(\rho k) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho k \mu\_i) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ (\mu + \frac{\mu\_t}{\sigma\_k}) \frac{\partial k}{\partial \mathbf{x}\_j} \right] + G\_k + G\_b - \rho \varepsilon - \mathbf{Y}\_M + S\_k \tag{1}$$

$$\frac{\partial}{\partial t}(\rho \varepsilon) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho \varepsilon u\_i) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ (\mu + \frac{\mu\_l}{\sigma\_\ell}) \frac{\partial \varepsilon}{\partial \mathbf{x}\_j} \right] + \mathcal{C}\_{1\varepsilon} \frac{\varepsilon}{k} (\mathcal{G}\_k + \mathcal{C}\_{3\varepsilon} \mathcal{G}\_b) - \mathcal{C}\_{2\varepsilon} \rho \frac{\varepsilon^2}{k} + \mathcal{S}\_{\varepsilon} \tag{2}$$

In the above equations, *Sk* and *S*<sup>ε</sup> are user-defined source terms. *YM* is the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. *Gk* is the production of turbulence kinetic energy due to the average speed gradients. *Gb* is the generation of turbulence kinetic energy due to buoyancy. *C*1ε, *C*2ε, and *C*3<sup>ε</sup> represents the model constants, σ*<sup>k</sup>* and σε are the turbulent Prandtl numbers for *k* and ε, respectively. The turbulent (or eddy) viscosity was computed by combining *k* and ε as follows:

$$
\mu\_t = \rho \mathbf{C}\_\mu \frac{k^2}{\varepsilon} \tag{3}
$$

where *C*<sup>μ</sup> is a constant. The model constants *C*1ε, *C*2ε, *C*μ, σ*<sup>k</sup>* and σε have default values of: *C*1<sup>ε</sup> = 1.44, *C*2<sup>ε</sup> = 1.92, *C*<sup>μ</sup> = 0.09, σ*<sup>k</sup>* = 1.0 and σε = 1.3.

#### *4.2. CFD Modeling Details Using STAR CCM*+

In a CFD simulation, the boundary condition "wall" is considered at locations where the stream cannot penetrate and includes walls, the ceiling, and the floor. The accompanying parameters were chosen for the model development: (1) the stream is incompressible, turbulent, and steady state; (2) water is selected as the working medium with 997.56 kg/m3 density; (3) The mass flow rate at each channel is 0.000277677 kg/s, while the aggregate mass flow rate at all nine channels is 0.002499003 kg/s; (4) the area at each channel is 5.272 <sup>×</sup> 10−<sup>7</sup> m2; (5) the dynamic viscosity of 0.00088871 Pa s; (6) the specific heat is 4181.72 J/kg K; (7) the thermal conductivity is 0.62 W/m K; and (8) the turbulent Prandtl number is 0.9. In addition to this, the thermal conductivity of the outlet aluminum cover is 237 W/m K, the density of cover is 2702 kg/m3, and specific heat is 903 J/kg K. The selected parameters for model set-up were (1) flow: turbulent; (2) fluid: incompressible; (3) time: steady state; (4) realizable K-epsilon (RANS); (5) two-layer wall: y+ wall treatment (y+ ≈ 5); (6) solver: segregated; (7) convection: second order; (8) turbulence intensity: 0.01 (default); and (9) turbulent viscosity ratio: 10.0 (default).

#### *4.3. Meshing in DEP MeshWorks 8.0*

The meshing of the area was a crucial step because different lattice parameters, such as quality criteria, mesh size, the shape of the elements, and the number of nodes have a significant impact on the result accuracy and the numerical solution. Here, the meshing was done using DEP MeshWorks 8.0 software. The screenshot of DEP MeshWorks during meshing of cooling plate is shown in Figure 4. Meshing in all nine-inlet channels of the cooling plate and meshing in the top portion of the cooling plate, which is specifically designed for this prismatic battery cooling, is shown in Figure 5. This design provides maximum cooling in this region because the heat production is the highest near the electrodes. In order to accurately represent the heat and flow transfer characteristics, the mesh was refined at regions of high geometrical deviations. Furthermore, P1, P2, and P3 areas for all cases for the heat flux sensor (HFS) is shown in Figure 6.

**Figure 4.** MeshWorks 8.0 screenshot during meshing of cooling plate.

**Figure 6.** Heat flux sensor (HFS) positions for model development.

#### **5. Results and Discussion**

This section presents on the results obtained from investigations for a specific prismatic LIB at various higher discharge currents of 60 A (3C) and 80 A (4C) for water cooling at working conditions of 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C.

#### *5.1. Temperature Contours at 3C (60 A) Discharge*

The temperature contours determined from STAR CCM+ CFD at 60 A discharge current and 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C working temperatures are shown in Figures 7–10. As previously mentioned, three heat flux sensors were put on the main surface of the battery: one was situated close to the positive terminal or cathode, the second was situated close to the negative terminal or anode, and the third was situated at the center of the cell. Figure 7 demonstrates the simulation results at 60 A discharge current

and 5 ◦C coolant inlet temperature with heat flux near cathode = 2347.7 W/m2, anode = 2259.5 W/m2, and mid surface = 539.3 W/m2. Similarly, Figure 8 shows temperature contours at 60 A discharge current and 15 ◦C coolant inlet temperature with heat flux values near cathode = 1711.8 W/m2, anode= 2351.6 W/m2, and mid surface = 548.4 W/m2. It is observed that the temperature contours and trends are similar with the inlet being cold and outlet being hot. During the battery operation at high C-rates, the generated heat from the battery is conducted to the cooling plate. As the coolant flows inside from the inlet, the heat is absorbed continuously with increments in coolant temperature along the flow path. The maximum temperature of coolant is observed at the outlet surface, as expected. Figure 9 shows temperature contours at 60 A discharge current and 25 ◦C coolant inlet temperature with heat flux values near cathode = 1597.3 W/m2, anode = 1851.6 W/m2, and mid surface = 413.0 W/m2. Figure 10 demonstrates temperature contours at 3C discharge rate and 35 ◦C coolant inlet temperature with heat flux values near cathode = 1468.4 W/m2, anode = 1579.9 W/m2, and mid surface = 340.6 W/m2.

**Figure 7.** Temperature profile at 60 A and 5 ◦C with heat flux values near cathode = 2347.7 W/m2, anode = 2259.5 W/m2, and mid surface = 539.3 W/m2. (**a**) Top view and (**b**) bottom view.

(**b**)

**Figure 8.** Temperature profile at 60 A and 15 ◦C with heat flux values near cathode = 1711.8 W/m2, anode = 2351.6 W/m2, and mid surface = 548.4 W/m2. (**a**) Top view and (**b**) bottom view.

(a) **Figure 9.** *Cont*.

**Figure 9.** Temperature profile at 60 A and 25 ◦C with heat flux values near cathode = 1597.3 W/m2, anode = 1851.6 W/m2, and mid surface = 413.0 W/m2. (**a**) Top view and (**b**) bottom view.

The result for temperature from simulation for 4C discharge rate and ambient temperature of 35 ◦C showed 38.12 ◦C, which is 2.85% higher or lower than experimental values for similar boundary conditions. The result for temperature from simulation for 3C discharge rate and ambient temperature of 5 ◦C showed 7.93 ◦C, which is 2.19% higher or lower than experimental values for similar boundary conditions. It is also noted that working temperature had a great impact on battery discharge capacity. As the working temperature rose from 5 ◦C to 35 ◦C, the temperature contour values also increased for a particular C-rate. The general cooling patterns are identical, demonstrating more noteworthy contrasts at the inlet of the cooling plate where the water is coldest. The temperatures differ with the inlet working condition temperature, yet the overall pattern remains generally consistent. Table 3 gives the summary of inlet and outlet water temperatures at 60 A discharge current and working temperature conditions of 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C. Table 3 also provides the difference in experimental and simulated values obtained from STAR CCM+ software. It is additionally seen that the simulated values are higher than experimental values because of the assumption of ideal wall conditions of adiabatic wall at non-heat-transferring boundaries. In the actual experiment, a small amount of heat transfer could be observed at some places where insulation was provided.



**Figure 10.** Temperature profile at 60 A and 35 ◦C with heat flux values near cathode = 1468.4 W/m2, anode = 1579.9 W/m2, and mid surface = 340.6 W/m2. (**a**) Top view and (**b**) bottom view.

#### *5.2. Temperature Contours at 4C (80 A) Discharge*

Figure 11 shows temperature contours at 80 A discharge current and 5 ◦C water inlet temperature with heat flux values near cathode = 3112.2 W/m2, anode = 3072.8 W/m2, and mid surface = 764.1 W/m2. Figure 12 shows temperature contours at 80 A discharge current and 15 ◦C water inlet temperature with heat flux values near cathode = 2419.0 W/m2, anode = 2887.1 W/m2, and mid surface = 697.3 W/m2. Figure 13 demonstrates temperature contours at 80 A discharge current and 25 ◦C water inlet temperature with heat flux values near cathode = 2309.3 W/m2, anode = 2648.2 W/m2, and mid surface = 611.1 W/m2. Figure 14 shows temperature contours at 80 A discharge current and 35 ◦C water inlet temperature with heat flux values near cathode = 2160.2 W/m2, anode = 2101.5 W/m2, and mid surface = 471.8 W/m2. It is noted that, as the battery discharged, the flowing water inside the cooling plate got heated because the heat was first conducted to the cooling plate and subsequently transferred to the coolant by convection. The joule heating is the dominant factor for heat generation. As the discharge current changed from 60 A to 80 A, there was an increase in temperature values as well. The pattern observed was that increased discharge currents and increased ambient temperatures

resulted in higher temperatures in the cooling plate. Table 3 gives the outline of inlet and outlet water temperatures at 80 A discharge current and different working temperature conditions of 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C. Table 3 also provides the difference in experimental and simulated values obtained from STAR CCM+ software. It was also found that the simulated values were higher than experimental values. In addition, the general cooling patterns were identical, similar to the results discussed in Section 5.1. There were noteworthy temperature differences at the inlet of the cooling plate where the water was coldest.

**Figure 11.** Temperature profile at 80 A and 5 ◦C with heat flux values near cathode = 3112.2 W/m2, anode = 3072.8 W/m2, and mid surface = 764.1 W/m2. (**a**) Top view and (**b**) bottom view.

**Figure 12.** Temperature profile at 80 A and 15 ◦C with heat flux values near cathode = 2419.0 W/m2, anode = 2887.1 W/m2, and mid surface = 697.3 W/m2. (**a**) Top view and (**b**) bottom view.

**Figure 13.** Temperature profile at 80 A and 25 ◦C with heat flux values near cathode = 2309.3 W/m2, anode = 2648.2 W/m2, and mid surface = 611.1 W/m2. (**a**) Top view and (**b**) bottom view.

**Figure 14.** Temperature profile at 80 A and 35 ◦C with heat flux values near cathode = 2160.2 W/m2, anode = 2101.5 W/m2, and mid surface = 471.8 W/m2. (**a**) Top view and (**b**) bottom view.

#### *5.3. Velocity Contours at 3C (60 A) and 4C (80 A) Discharge*

The investigations of velocity contours can provide insights into future design considerations by comparing contour results of velocity and comparing with respective contour results of temperatures. The velocity contours at 60 A and 80 A discharge currents and 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C working temperatures appear in Figures 15 and 16. The velocity contours were identical in all the cases, in accordance with general trends, given the low temperatures associated in the modeling that would have had a minimal impact on the water density. These results might be influenced by the lower y+ value, wall functions and turbulence model utilized. It was also observed that the velocity distribution

at the inlet to the cooling plate and the outlet from the cooling plate was curved with relatively higher velocity gradients.

**Figure 15.** Velocity profile at 60 A with 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C.

**Velocity profile at 80 A\_5 °C Velocity profile at 80 A\_15 °C** 

**Figure 16.** *Cont*.

**Figure 16.** Velocity profile at 80 A with 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C.

#### *5.4. Transient Temperature Profiles of Water Flow and Voltage Distributions*

Figures 17 and 18 show the transient behavior of water flowing inside the cooling plates at 60 A and 80 A constant current discharges with various working temperatures of 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C. As discussed earlier, the increase in temperature was due to the joule heating *I 2R* from the LIB during discharge.

**Water temperature profile at 60 A\_25 °C Water temperature profile at 60 A\_35 °C** 

**Figure 17.** Transient temperature profile of water flow at 60 A with 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C.

*Energies* **2020**, *13*, 1638

**Figure 18.** Transient temperature profile of water flow at 80 A with 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C.

It was discovered that the working temperature had a great impact on the battery performance. At lower discharge currents, the battery capacity was close to the manufacturer's supplied data sheet, but as discharge current increased, there was a decrease in the discharge capacity. Further, when the working temperature changed from 35 ◦C to 5 ◦C, there was a more prominent decrease in the discharge capacity. Consequently, it is clear that as the working temperature decreased, the battery discharge capacity also decreased. These effects (reduction in battery discharge capacity) can be seen in Figures 19 and 20, which present the discharge/charge profiles at 60 A and 80 A constant current discharges (and charge current being 20 A) with various working temperatures of 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C.

**Discharge/charge voltage profile at 60 A\_5 °C Discharge/charge voltage profile at 60 A\_15 °C** 

*Energies* **2020**, *13*, 1638

**Discharge/charge voltage profile at 60 A\_25 °C Discharge/charge voltage profile at 60 A\_35 °C** 

**Figure 19.** Discharge/charge voltage profile at 60 A with 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C.

**Discharge/charge voltage profile at 80 A\_25 °C Discharge/charge voltage profile at 80 A\_35 °C** 

**Figure 20.** Discharge/charge voltage profile at 60 A with 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C.

#### **6. Conclusions**

This paper presented a numerical model using STAR CCM+ for CFD simulations at high C-rates and diverse working temperatures of the liquid (water) with 5 ◦C, 15 ◦C, 25 ◦C, and 35 ◦C. We discovered that the temperature distributions within cooling plate channels increased with C-rates (3C to 4C). As C-rate increased, the heat flux values measured near the anode, the cathode, and the middle surface also increased. The cooling patterns obtained from simulation were similar to the experimental values with slightly higher values. The velocity plots were identical for all cases. There results provide valuable information on the design considerations that must be made for battery cooling systems in EVs.

**Author Contributions:** Conceptualization, S.P., Formal analysis, M.-K.T.; Supervision, R.F. and M.F.; Validation, K.G.; Writing—original draft, S.P.; Writing—review and editing, M.-K.T. and S.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This research work (reverse engineering, CT scanning, MeshWorks and STAR CMM+ software) was done at Detroit Engineering Products (DEP), Michigan, USA, during year 2018–2019. The experimental data was used from University of Waterloo for the model development and validation. This work was supported by equipment and manpower from the Department of Chemical Engineering at the University of Waterloo. Special thanks to Danielle Skeba for the contribution in editing the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**



#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Review*
