**Small-Signal Modeling and Analysis for a Wirelessly Distributed and Enabled Battery Energy Storage System of Electric Vehicles**

#### **Yuan Cao**

Department of Electrical and Computer Engineering, College of Engineering, The University of Alabama, Tuscaloosa, AL 35487, USA; ycao19@crimson.ua.edu

Received: 31 August 2019; Accepted: 7 October 2019; Published: 11 October 2019

#### **Featured Application: The major object of this research is to provide a small-signal modeling method and controller design guidelines in wireless distributed and enabled battery energy storage system (WEDES) battery system for electric vehicles applications.**

**Abstract:** This paper presents small-signal modeling, analysis, and control design for wireless distributed and enabled battery energy storage system (WEDES) for electric vehicles (EVs), which can realize the active state-of-charge (SOC) balancing between each WEDES battery module and maintain operation with a regulated bus voltage. The derived small-signal models of the WEDES system consist of several sub-models, such as the DC-DC boost converter model, wireless power transfer model, and the models of control compensators. The small-signal models are able to provide deep insight analysis of the steady-state and dynamics of the WEDES battery system and provide design guidelines or criteria of the WEDES controller. The derived small-signal models and controller design are evaluated and validated by both MATLAB®/SIMULINK simulation and hardware experimental prototype.

**Keywords:** small-signal modeling; battery energy storage system; battery management system; control; stability; dynamic response; wireless power; state-of-charge; electric vehicle

#### **1. Introduction**

Battery energy storage systems (BESS) have been widely used in various applications, such as electric vehicles (EVs), consumer electronics, medical devices, smart grid, energy backup in data centers, and among others [1–10]. For EV applications, what is referred to be as "range anxiety" is one of the major reasons that prohibit or slows down the adoption of EVs [9–14].

To eliminate range anxiety in and extend the driving range of EVs, different methods have been discussed in the literature [9–17], such as increasing battery pack capacity, utilizing a faster charging method, utilizing a battery pack swapping method, achieving dynamic wireless charging, etc. While these methods can be effective to extend the driving range of EVs, some design challenges or drawbacks cannot be ignored.

When increasing capacity, the weight, size, and cost of the battery pack inevitably increase with the increase in battery capacity [10]. Further, the needed recharge time is also increased. For faster charging, the battery state-of-health (SOH) degrades at a higher rate if faster charging is applied [11,12]. In addition, the fast charger requires a high-power infrastructure that increases the cost of the overall system. For conventional battery swapping, specialized equipment, as well as the experienced personnel, are required to realize battery swapping [13]. For dynamic wireless charging, a large number of transmitter (Tx) coils are required with corresponding power supply units, which increases the infrastructure cost. In addition, this method might not be practical in all locations [14].

Among these methods, the battery swapping concept is a good candidate to reduce recharging time and extend the driving range with low infrastructure cost. To deal with the challenges in the conventional battery swapping concept, a new distributed and enabled battery energy storage (WEDES) system and WEDES controller for EVs are presented in [1], which allows for fast and safe exchange/swapping of smaller and lighter battery modules with wireless power transfer (WPT) technology [18,19].

An illustration of the WEDES system for EVs is shown in Figure 1b, and its example circuit diagram is shown in Figure 1c. Each of the battery modules consists of multiple battery cells, a dedicated electronics circuit, wireless power transmitter coil (Tx coil), wireless communication circuit, and client controller. While the on-board-unit (OBU) consists of a wireless power receiver coil (Rx coil), wireless communication circuit, and host controller. Different from the conventional battery swapping concept where the battery pack as a whole is exchanged at one time, in the WEDES battery system, the conventional single battery pack is divided into multiple small battery modules, which can deliver power through wireless power transfer (WPT) technology to the OBU. The distributed nature of the WEDES system combined with wireless power transfer (no physical connection between battery modules and OBU) makes the battery exchange/swapping easier, safer, and faster.

The distributed WEDES battery system with the WEDES controller addresses state-of-charge (SOC) balancing, bus voltage regulation, and battery module current/voltage regulation at the same time inside the system. Therefore, an SOC balancing control loop, a bus voltage regulation control loop, and a battery module current/voltage control loop are coupled with each other within one battery module as well as between multiple battery modules. These couplings make the analysis and design of the WEDES controller complex and critical. While the initial concept of the distributed WEDES battery system is discussed in [1], the small-signal modeling and controller design analysis are not focused on.

The main contributions of this paper can be summarized as follows:


The next Section discusses the detailed small-signal derivation of the WEDES system. Section 3 presents the design of compensators for each control loop. Simulation models and experimental results are presented and discussed in Section 4 to validate the derived small-signal models. Section 5 is the additional comments, and Section 6 concludes the paper.

**Figure 1.** *Cont.*

**Figure 1.** Illustration diagrams of battery system for electric vehicle (EV) application. (**a**) The conventional battery pack and electrics drive system in EVs, (**b**) the wireless distributed and enabled battery energy storage (WEDES) battery system in EVs, and (**c**) example circuit diagram of the WEDES system [1].

#### **2. Small-Signal Modeling of the Distributed WEDES System**

#### *2.1. Overview of the WEDES System and Controller Operation Principle*

Figure 1c shows the illustration of an example circuit diagram of the WEDES system, which consists of two major parts: battery modules and on-board-unit (OBU).

Inside each battery module, multiple battery cells are connected in series and/or in parallel to form a battery bank, which can provide voltage/current/power to the rest of the system. The output of the battery bank is then connected to the input of a DC-DC power converter, which is used to achieve bus voltage regulation, battery module current/voltage regulation as well as SOC balancing at the same time as described later in this section. The output of the power converter is connected to an inverter stage for DC-AC power conversion. At the end of the battery module, the AC power from the inverter is applied to the transmitter (Tx) for inductive wireless power transfer (I-WPT) to the OBU.

The OBU mainly consists of receiver coils (Rx) followed by an AC-DC power conversion/rectification stage (rectifier). The outputs of each rectifier (*V*o1 through *V*oN) are connected in series to the bus/output (*V*bus = *V*o1 + *V*o2 + ... + *V*oN).

To realize the functionalities of SOC balancing, bus voltage regulation, and battery module current/voltage regulation, the WEDES controller consists of three different control loops: the SOC balancing control loop (referred to be by the SOC balancing loop), battery module voltage control loop (referred to be by the module voltage loop) and bus voltage control loop (referred to be by the bus voltage loop). Figure 2 shows the diagram of the wirelessly distributed WEDES controller, where *V*bus\_ref is the desired value of bus voltage, VMN-total is an intermediate value for voltage regulation, *C*\_re\_MX1 though *C*\_re\_MXN are the remaining capacities of battery modules for SOC calculation, *SOC*MX1 though *SOC*MXN are the SOC values of battery modules, λDC1 through λDCN are weighting factors to generate the reference values of *V*MX1-DC-ref through *V*MXN-DC-ref for each battery module, αMX1 through αMXN are the SOC multipliers and δMX1 through δMXN are enable/disable values.

**Figure 2.** Illustration of the WEDES controller presented in [1].

In the WEDES controller, the SOC balancing loop is used to generate multipliers αMX1 through αMXN to realize SOC balancing between multiple battery modules, as given by Equation (1).

⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩ α*MX*<sup>1</sup> = (*SOCMX*−*re f* − *SOCMX*1) × *GvSOC* + 1 α*MX*<sup>2</sup> = (*SOCMX*−*re f* − *SOCMX*2) × *GvSOC* + 1 ...... . α*MXN* = (*SOCMX*−*re f* − *SOCMXN*) × *GvSOC* + 1 , (1)

where *SOC*MX-ref is the average SOC value of all battery modules, as given by Equation (2). When all battery modules are inserted and active, the sum of all δMX1 through δMXN equals to *N* (δ*MX*<sup>1</sup> + δ*MX*<sup>2</sup> + ... + δ*MXN* = *N*).

$$\text{SOC}\_{\text{MX-ref}} = \frac{(\delta\_{\text{MX1}} \times \text{SOC}\_{\text{MX1}}) + (\delta\_{\text{MX2}} \times \text{SOC}\_{\text{MX2}}) + \dots + (\delta\_{\text{MXN}} \times \text{SOC}\_{\text{MXN}})}{\delta\_{\text{MX1}} + \delta\_{\text{MX2}} + \dots + \delta\_{\text{MXN}}},\tag{2}$$

$$
\lambda\_{\rm DC} = \frac{\delta\_{\rm MXr} \times a\_{\rm MXr}}{(\delta\_{\rm MX1} \times a\_{\rm MX1}) + (\delta\_{\rm MX2} \times a\_{\rm MX2}) + \dots + (\delta\_{\rm MXN} \times a\_{\rm MXN})} \tag{3}
$$

If the SOC value of *r*th battery module is larger than others, the corresponding multiplier αMXr will be set larger, and vice versa. These multipliers αMX1 through αMXN are then multiplied by enabled/disable values δMX1 through δMXN to further generate the weighting factors λDC1 through λDCN, as given by Equation (3). The sum of weighting factors λDC1 through λDCN always equals to one. These weighting factors are then used in the battery module voltage loop to regulate the output voltage of the battery modules *V*MX1 through *V*MXN at the primary side (Tx side), and as a result, achieve bus voltage regulation at the second side (Rx/OBU side).

It should be emphasized that due to the inevitable power loss during wireless power transfer (i.e., transmission efficiency is less than 100%), the bus voltage control loop is important to adaptively adjust *V*MN-total value to compensate the conversion ratios and losses of multiple power conversion stages (DC-AC-AC-DC) and realize bus voltage regulation. The relationship between different voltages can be calculated as given by Equation (4).

$$\begin{cases} \begin{aligned} V\_{MX-total} &= \left(V\_{bus-ref} - V\_{bus}\right) \times G\_{v-bus} \\ V\_{MX-total} &= \left.V\_{MX1} + \left.V\_{MX2} + \ldots + \left.V\_{MXN}\right| \right. \end{aligned} \right. \\\begin{aligned} \begin{cases} V\_{MXr-DC-ref} = V\_{MX-total} \times \lambda\_{DCr} \\ V\_{bus} = V\_{o1} + \left.V\_{o2} + \ldots + \left.V\_{oN}\right| \end{cases} \end{cases} \end{cases} \tag{4}$$

To summarize, the presented WEDES controller can dynamically control SOC multipliers αMX1 through αMXN to adjust the discharging rate for each battery module to achieve SOC balancing, while keeping λDC1 + λDC2<sup>+</sup> ... +λDCN = 1 such that the bus voltage is always regulated as *V*bus-ref.

#### *2.2. Small-signal Modeling*

Based on the block diagram of the WEDES controller shown in Figure 2, the small-signal of the distributed WEDES system with controller is shown in Figure 3. The transfer functions and symbols in Figure 3 are summarized as follows. For simplicity, the *r*th battery module is used for illustration.

*Lbus*(*s*): Bus voltage control loop gain;

*LMXr*(*s*): Battery module voltage control loop gain;

*LSOCr*(*s*): SOC balancing control loop gain;

*Gvdr*(*s*): Duty cycle to DC-DC converter output voltage *V*MXN transfer function;

*GidN*(*s*): Duty cycle to cell current transfer function;

*Gsocir*(*s*): Cell current to cell SOC transfer function;

*GiTR*(*s*): Gain of the input current of the half-bridge inverter to the output current of rectifier;

*GvTR*(*s*): Gain of the input voltage of the half-bridge inverter to the output voltage of rectifier; *GPWM*: PWM module gain;

*Kdivr*(*s*): Output voltage sensing gain (including voltage divider gain and ADC conversion gain);

*Delaywr*(*s*): Delay of wireless communication;

*Delaydr*(*s*): Delay of digital computation;

*ZOHvr*(*s*): Zero order hold model of voltage sampling;

*ZOHir*(*s*): Zero order hold model of current sampling;

*ZOHSOCr*(*s*): Zero order hold model of SOC sampling.

**Figure 3.** Small-signal of distributed the WEDES battery system with the WEDES controller.

#### *2.3. Derivation of Transfer Functions*

Since the design parameters and equilibrium operation point of all WEDES battery modules are the same under steady-state operation (when SOC balancing is achieved), the derivation of transfer functions of all battery modules follows the same procedure. The detailed derivation of *r*th battery module is discussed as follows:

#### 2.3.1. Transfer Function of DC-DC Boost Stage

The circuit diagram of the DC-DC boost stage is shown in Figure 4. When the lower side switch SL-r is on and the upper side switch SU-r is off, the differential equation of the boost converter can be derived as follows:

$$\begin{cases} \begin{array}{c} L\_{M\mathcal{X}r} \frac{di\_{\rm irr}}{dt} = v\_{\rm irr} \\ \mathcal{C}\_{M\mathcal{X}r} \frac{d\mathcal{V}\_{M\mathcal{X}r}}{dt} = -i\_{M\mathcal{X}r} \end{array} \end{cases} \tag{5}$$

where *iinr* and *vinr* are the input current and input voltage, respectively, and *iMXr* and *vMXr* are the output current and output voltage of the boost converter, respectively. *L*MXr is the inductor value, and *C*MXr is the output capacitor.

**ͲŽŽƐƚ^ƚĂŐĞ**

**Figure 4.** Circuit diagram of DC-DC boost stage in the WEDES system.

The state-space form of Equation (5) can be rewritten as Equation (6).

$$K\begin{bmatrix} \dot{\mathbf{x}}\_1\\ \dot{\mathbf{x}}\_2 \end{bmatrix} = A\_1 \begin{bmatrix} \mathbf{x}\_1\\ \mathbf{x}\_2 \end{bmatrix} + B\_1 \begin{bmatrix} \boldsymbol{u}\_1\\ \boldsymbol{u}\_2 \end{bmatrix} = \begin{bmatrix} 0 & 0\\ 0 & 0 \end{bmatrix} \begin{bmatrix} \mathbf{x}\_1\\ \mathbf{x}\_2 \end{bmatrix} + \begin{bmatrix} 1 & 0\\ 0 & -1 \end{bmatrix} \begin{bmatrix} \boldsymbol{u}\_1\\ \boldsymbol{u}\_2 \end{bmatrix}.\tag{6}$$

where the variables are *x*<sup>1</sup> = *iinr*, *x*<sup>2</sup> = *vMXr*, *u*<sup>1</sup> = *vinr*, *u*<sup>2</sup> = *iMXr*, respectively. *K* = *L* 0 0 *C* .

Similarly, when the lower side switch SL-r is off and the upper side switch SU-r is on, the differential equation of the boost converter can be derived as follows:

$$\begin{cases} \begin{array}{c} L\_{M\mathbf{X}r} \frac{d\dot{l}\_{\mathrm{inv}}}{dt} = V\_{\mathrm{inv}} - V\_{M\mathbf{X}r} \\ \mathbf{C}\_{M\mathbf{X}r} \frac{d\dot{V}\_{M\mathbf{X}r}}{dt} = \dot{\mathbf{i}}\_{\mathrm{inv}} - I\_{M\mathbf{X}r} \end{array} . \end{cases} \tag{7}$$

The state-space form of Equation (7) can be rewritten as Equation (8).

$$K\begin{bmatrix} \dot{\mathbf{x}}\_1\\ \dot{\mathbf{x}}\_2 \end{bmatrix} = A\_2 \begin{bmatrix} \mathbf{x}\_1\\ \mathbf{x}\_2 \end{bmatrix} + B\_2 \begin{bmatrix} \boldsymbol{u}\_1\\ \boldsymbol{u}\_2 \end{bmatrix} = \begin{bmatrix} 0 & -1\\ -1 & 0 \end{bmatrix} \begin{bmatrix} \mathbf{x}\_1\\ \mathbf{x}\_2 \end{bmatrix} + \begin{bmatrix} 1 & 0\\ 0 & -1 \end{bmatrix} \begin{bmatrix} \boldsymbol{u}\_1\\ \boldsymbol{u}\_2 \end{bmatrix} \tag{8}$$

By using the state-space averaging method [20], the average matrix A and B are calculated as given by Equations (9) and (10).

$$A = A\_1 D\_7 + A\_2 (1 - D\_r),\tag{9}$$

$$B = B\_1 D\_{\mathcal{T}} + B\_2 (1 - D\_{\mathcal{T}}),\tag{10}$$

where *D*r is the duty cycle of the boost converter. The steady-state X is calculated as follows:

$$X = A^{-1} B I I = \begin{bmatrix} \frac{I \text{M} \text{X}r}{1 - D\_r} & V\_{inv} \end{bmatrix}^{-1}. \tag{11}$$

The small-signal equation becomes as given by Equation (12):

$$
\hat{X} = A\hat{\mathbf{x}} + B\hat{\boldsymbol{\mu}} + \left[ (A\_1 - A\_2)X + (B\_1 - B\_2)\mathbb{U} \right] \hat{d}\_{r\_2} \tag{12}
$$

where *A*<sup>1</sup> − *A*<sup>2</sup> = 0 1 <sup>−</sup>1 0 , *<sup>B</sup>*<sup>1</sup> <sup>−</sup> *<sup>B</sup>*<sup>2</sup> <sup>=</sup> 0, and <sup>ˆ</sup> *dr* is the small signal variation of duty cycle around its steady state operation point.

Equation (12) can be rewritten as Equation (13):

$$\begin{cases} L\_{M\mathbf{X}r} \frac{d\hat{i}\_{\mathrm{nr}r}}{dt} = -(1 - D\_r) \upsilon\_{M\mathbf{X}r} + \upsilon\_{\mathrm{nr}r}^\circ + \frac{V\_{\mathrm{nr}r}}{1 - D\_r} \hat{d}\_r\\ \quad \mathbb{C}\_{M\mathbf{X}r} \frac{d\hat{V}\_{M\mathbf{X}r}}{dt} = (1 - D\_r) \hat{i}\_{\mathrm{nr}r}^\circ - \dot{\mathbf{i}}\_{M\mathbf{X}r}^\circ - \frac{I\_{M\mathbf{X}r}}{1 - D\_r} \hat{d}\_r \end{cases} \tag{13}$$

To simply the analysis, the AC small-signal variation of *v* ˆ*inN* and ˆ *iMXN* is assumed to be 0 (negligible) because the dynamic variation of battery voltage and battery module output current is very slow compared to the dynamic variation of the control signal ˆ *dr* (duty cycle) of the power converter. Therefore, by performing the Laplace transformation, Equation (14) can be derived as:

$$\begin{cases} s\mathcal{L}\_{\rm MXr} i\_{\rm irr}(s) = -(1 - D\_r)\upsilon\_{\rm MXr}(s) + \frac{V\_{\rm irr}}{1 - D\_r}d\_r(s) \\\ s\mathcal{C}\_{\rm MXr}\upsilon\_{\rm MXr}(s) = (1 - D\_r)i\_{\rm irr}(s) - \frac{I\_{\rm MXN}}{1 - D\_r}d\_r(s) \end{cases} . \tag{14}$$

Based on Equation (15), the output voltage to the control signal transfer function of the power converter can be derived as

$$\mathbf{G}\_{\rm rdr} = \frac{v\_{\rm MXr}(s)}{d\_r(s)} = \frac{\frac{1}{\overline{\mathbf{C}\_{\rm MXr}}} \left(-\frac{I\_{\rm MXr}}{1 - D\_r}s + \frac{V\_{\rm inr}}{\overline{\mathbf{C}\_{\rm MXr}}}\right)}{s^2 + \frac{\left(1 - D\_r\right)^2}{L\_{\rm MXr}\overline{\mathbf{C}\_{\rm MXr}}}}. \tag{15}$$

Similarly, the input current to control signal transfer function can be derived as

$$\mathbf{G}\_{\rm idr} = \frac{i\_{\rm irr}(\mathbf{s})}{d\_r(\mathbf{s})} = \frac{\frac{1}{L\_{\rm MKr}} \left(\frac{I\_{\rm MKr}}{C\_{\rm MKr}} + \frac{V\_{\rm irr}}{1 - D\_r}\mathbf{s}\right)}{\mathbf{s}^2 + \frac{\left(1 - D\_r\right)^2}{L\_{\rm MKr}C\_{\rm MKr}}}.\tag{16}$$

#### 2.3.2. Transfer Function of WPT Stage (Half-Bridge Inverter, WPT Coils, and Half-Bridge Rectifier)

Figure 5 shows the circuit diagram of the WPT stage. By writing the Kirchhoff's voltage law (KVL) equations as given by Equation (17), the ratio between the output voltage *V*Rr and the input voltage *V*Tr at the resonance frequency can be calculated as given by Equation (18).

$$
\begin{bmatrix}
j\omega L\_{Tr} + \frac{1}{j\omega C\_{Tr}} + R\_{pTr} & -j\omega M\_{TR} \\
\end{bmatrix}
\begin{bmatrix}
I\_{Tr} \\
I\_{Rr}
\end{bmatrix} = \begin{bmatrix}
V\_{Tr} \\
0
\end{bmatrix},
\tag{17}
$$

$$\left. \text{G}\_{\text{r}TR} \right|\_{\omega = \omega\_{\text{b}}} = \frac{V\_{Rr}}{V\_{\text{Tr}}} = \frac{j\omega M\_{\text{TR}} R\_{\text{Lr}}}{Z\_{\text{Tr}} Z\_{\text{R}r} + \left(\alpha M\_{\text{TR}}\right)^2} = \frac{j\alpha k\_{\text{TR}} \sqrt{\mathcal{L}\_{\text{Tr}} L\_{\text{R}}} R\_{\text{Lr}}}{Z\_{\text{Tr}} Z\_{\text{R}r} + \left(\alpha k\_{\text{TR}}\right)^2 L\_{\text{Tr}} L\_{\text{R}r}},\tag{18}$$

where *ZTr* = *j*ω*LTr* + <sup>1</sup> *<sup>j</sup>*ω*CTr* + *RpTr* is the equivalent impedance of Tx side and *ZRr* = *<sup>j</sup>*ω*LRr* + 1 *<sup>j</sup>*ω*CRr* + *RpTr* + *RLr* is the equivalent impedance of Tx side. By substituting *<sup>s</sup>* = *<sup>j</sup>*<sup>ω</sup> into Equation (19), the following transfer function can be obtained.

$$\mathbf{G}\_{\rm vTRr}(\mathbf{s}) = \frac{V\_{\rm Rr}(\mathbf{s})}{V\_{\rm Tr}(\mathbf{s})} = -\frac{s k\_{\rm TRx} \sqrt{L\_{\rm Tr} L\_{\rm Rr}} R\_{\rm Lr}}{Z\_{\rm Tr} Z\_{\rm Rr} + s^2 k\_{\rm TRx}^2 L\_{\rm Tr} L\_{\rm Rr}}.\tag{19}$$

**Figure 5.** Circuit diagram of the wireless power transfer (WPT) stage in the WEDES system (**a**) equivalent circuit and (**b**) T-model.

Similarly, the ratio between the Tx current and Rx current can be derived as given by Equation (20).

$$\mathbf{G}\_{\rm ITRr}(\mathbf{s}) = \frac{I\_{\rm Rr}(\mathbf{s})}{I\_{\rm Tr}(\mathbf{s})} = -\frac{\mathbf{s}^2 k\_{\rm TRr} \mathbf{C}\_{\rm Rr} \sqrt{\mathbf{L}\_{\rm Tr} \mathbf{L}\_{\rm Rr}}}{\mathbf{s}^2 \mathbf{L}\_{\rm Rr} \mathbf{C}\_{\rm Rr} + \mathbf{s} \mathbf{C}\_{\rm Rr} (\mathbf{R}\_{\rm gTr} + \mathbf{R}\_{\rm Lr}) + 1}. \tag{20}$$

#### **3. Compensator Design**

#### *3.1. Battery Module Voltage Control Loop Compensator Design*

The WEDES system design parameters are shown in Table 1. Based on the small-signal model in Figure 3, the uncompensated discrete-time transfer function of the battery module voltage loop for the *r*th battery module *G*busr (Z) consists of *G*PWM, *G*vdr, *G*vTRr, *k*divr, *ZOH*vr, and digital computation delay *Delay*dr. *G*busN (Z) is calculated as given by Equation (21), and its bode plot is shown as the dashed curve in Figure 6.

$$L\_{\rm MXr\_{unmp}}(z) = Z \{ \mathbf{G}\_{\rm PWM}(\mathbf{s}) \cdot \mathbf{G}\_{\rm vdr}(\mathbf{s}) \cdot \mathbf{G}\_{\rm vTR}(\mathbf{s}) \cdot \mathbf{ZOH}\_{\rm vN}(\mathbf{s}) \cdot \mathbf{D} \rm lay\_{dN} \} = \frac{-0.0001999 \mathbf{z} - 0.000272}{x^3 - 1.977x^2 + 0.9789z}, \tag{21}$$

where *ZOHvr*(*s*) = <sup>1</sup>−*e*−*s*·*Ts <sup>s</sup>* ; *Delaydr*(*s*) = *e* <sup>−</sup>*sTdelay* ; *Tdelay* is the digital controller computation delay and it is equal to 10 μs in the experimental implementation; *GPWM* = 1/1024 ; *Kdiv* = 11 with 1 kΩ and 10 kΩ resistors as the voltage divider.


**Table 1.** Design parameters of the wireless distributed and enabled battery energy storage (WEDES) system.

**Figure 6.** The bode plot of the uncompensated (red-dashed curve) and compensated (blue-solid curve) battery module voltage control loop gain.

In the battery module voltage compensator design, there is a right-hand-plane-zero (RHPZ), which is located at 4.86 kHz. This RHPZ is introduced due to the existing boost converter. To guarantee the stability of the system, the compensated control bandwidth should be smaller than RHPZ. With a compensator Gv-bus(z) given by (22), the compensated battery module output voltage control loop gain (*GMXr*\_*comp* = *GMXr*\_*uncomp*(z)·*GMXr*(z)) achieves a control bandwidth of 1.62 kHz and a phase margin of 45◦, as shown on the solid curve in Figure 6.

$$\mathbf{G}\_{\rm MXr}(\mathbf{z}) = \frac{950.8\mathbf{z}^2 - 1850\mathbf{z} + 899.7}{\mathbf{z}^2 - \mathbf{z}}\tag{22}$$

#### *3.2. SOC Balancing Loop Compensator Design*

According to the small-signal model shown in Figure 3, the uncompensated SOC loop gain (i.e., with unity SOC loop compensator gain) is given by Equation (23), and its bode plot is represented by the dashed curve in Figure 7.

$$\mathcal{L}\_{\text{sacr-uncomp}}(z) = \mathcal{Z}\{\mathcal{G}\_{\text{kIN}}(s) \cdot \mathcal{Z}OH\_{\text{kellIN}}(s) \cdot \mathcal{G}\_{\text{saclN}}(s) \cdot \mathcal{Z}OH\_{\text{sacN}}(s) \cdot \mathcal{G}\_{\text{sacN}}(-1) \cdot \delta\_{\text{MXN}} \cdot \mathcal{G}\_{\text{a1}} \cdot \mathcal{D}\text{elay}\_{\text{b1}} \cdot \mathcal{G}\_{\text{b1-lus}}\},\tag{23}$$

where *ZOHicellN*(*s*) <sup>=</sup> <sup>1</sup>−*e*−*s*·*Ts <sup>s</sup>* , *ZOHsocN*(*s*) <sup>=</sup> <sup>1</sup>−*e*−*s*·*Tsoc <sup>s</sup>* , *T*soc is the sampling period for the SOC value in the SOC balancing loop. Since the SOC value of a battery cell varies very slowly compared to the switching period of the power converter, the sampling rate of the SOC balancing loop does not have to be very fast. *T*soc =1 s is found to be a good trade-off between the hardware resource consumption, system stability and SOC balancing speed.

**Figure 7.** The bode plot of uncompensated (red-dashed curve) and compensated (blue-solid curve) SOC balancing control loop.

With a compensator given by Equation (24), the compensated SOC balancing loop gain achieves a control bandwidth of 0.057 Hz and phase margin of 59.2◦, as shown on the solid curve in Figure 7. Due to the slow sampling rate of SOC value (1 Hz), it is expected that the control bandwidth of the SOC balancing loop is much lower than that of the voltage control loop.

$$G\_{\rm SOC}(z) = \frac{2.5 \times 10^8}{z - 1} \tag{24}$$

#### *3.3. Bus Voltage Control Loop Compensator Design*

Based on the small-signal model shown in Figure 3, the uncompensated bus voltage control loop gain (i.e., with gain = 1) is given by Equation (25), and its bode plot is represented by the dashed curve in Figure 8. With a compensator given by Equation (26), the compensated bus voltage loop gain achieves a control bandwidth of 8.12 kHz and a phase margin of 52.3◦, as shown on the solid curve in Figure 8.

$$L\_{\text{bus-uncomp}}(z) = Z \{ ZOH\_{v-bus}(s) \cdot Delay\_{\text{bw}} \} = \frac{0.2585z - 0.1798}{z^2 - 1.977z + 0.9789},\tag{25}$$

$$G\_{\rm bus}(z) = \frac{1.2z - 1.8}{z - 1},\tag{26}$$

where *ZOHv*−*bus*(*s*) <sup>=</sup> <sup>1</sup>−*e*−*s*·*Ts <sup>s</sup>* .

**Figure 8.** The bode plot of the uncompensated (red-dashed curve) and compensated (blue-solid curve) bus voltage control loop gain.

#### **4. Simulation and Model Experiment Validation**

In this section, the small-signal model was evaluated and validated in both simulation and hardware experiments. The simulation model was built in MATLAB®/SIMULINK software (2018a, MatchWorks, Natick, Massachusetts, USA) using the derived transfer functions in Sections 2 and 3. The hardware control compensator was implemented with Texas Instrumental microcontroller TMS320S28335. The design parameters are shown in Table 1. In the hardware experiment, three WEDES battery modules were implemented and utilized.

For the verification of derived small-signal models of the WEDES system, the dynamic response of both simulation and hardware experiments under different operation conditions was compared. If the dynamic performance from simulation results and experimental results are in good agreement, it can be implied that the developed model is valid.

As discussed in previous sections, there are three different control loops (battery module voltage control loop, bus voltage control loop, and SOC balancing control loop) in the WEDES system. During the test for a specific control loop, one of the operation parameters was changed while the rest of the operation parameters were set constant.

#### *4.1. Experimental Results for Battery Module Voltage Control Loop*

For the test of battery module voltage control loop, the reference output voltage for each battery module was changed suddenly from 10 V (VMX1-DC-ref = VMX2-DC-ref = VMX3-DC-ref = 10 V under steady-state operation) to VMX1-DC-ref = 13 V, VMX2-DC-ref = 9 V and VMX3-DC-ref = 8 V. The simulation results and experimental results are shown in Figure 9.

**Figure 9.** Waveforms for three battery module output voltage when the reference voltages were changed from VMX1-DC-ref = VMX2-DC-ref = VMX3-DC-ref = 10 V to VMX1-DC-ref = 13 V, VMX2-DC-ref = 9 V and VMX3-DC-ref = 8 V, (**a**) simulation results and (**b**) hardware experimental results.

From Figure 9, it can be observed that once the reference voltage values were changed, the output voltage of battery module 1 was controlled to increase while the output voltages of battery modules 2 and 3 were controlled to decrease. In Figure 9, the simulation model results and the hardware experimental results are in good agreement. In other words, the shape, magnitude, and overshoot/undershoot of the waveforms of simulation and hardware experiments match each other, which validates the small-signal model for the battery module voltage loop.

#### *4.2. Experimental Results for SOC Balancing Control Loop*

For the test of SOC balancing control loop, the SOC value of three battery modules were changed suddenly from SOC1 = SOC2 = SOC3 = 70% under balanced conditions to SOC1 = 75%, SOC2 = 70%, and SOC3 = 65%. The experimental results are shown in Figure 10.

It can be observed from Figure 10 that once there was a change in SOC values, Vo1 increased while Vo3 decreased. This is because the value of SOC1 was larger than the average SOC value of three battery modules (i.e., (75% + 70% + 65%)/3 = 70%), therefore resulted in a larger value of the battery module output voltage Vo1 for faster discharge. The variation of Vo3 for battery module 3 followed the same behavior but in the opposite direction. The value of SOC2 was equal to the average SOC value, and therefore, its corresponding battery module output voltage remained constant. Figure 10 also shows that the simulation results and hardware experimental results are in good agreement under the variation of SOC values test condition.

**Figure 10.** Waveforms for three battery module output voltage when SOC values of three modules changed from SOC1 = SOC2 = SOC3 = 70% to SOC1 = 75%, SOC2 = 70%, and SOC3 = 65%, (**a**) simulation results and (**b**) hardware experimental results.

#### *4.3. Experimental Results for Bus Voltage Control Loop*

The last test was implemented for the bus voltage control loop. In this test, the bus voltage reference value was changed suddenly from 30 V to 37.5 V, as shown in Figure 11. It can be observed from Figure 11 that as the bus voltage reference value increased, the battery module voltage values changed from 10 V to 12.5 V correspondingly, which yielded a total bus voltage equal to its new reference value. The waveforms shown in Figure 11 demonstrate the consistency between simulation results and hardware experimental results.

**Figure 11.** Waveforms for three battery module output voltage when bus reference voltage changes from Vbus-ref = 30 V to Vbus-ref = 37.5 V, (**a**) simulation results and (**b**) hardware experimental results.

#### **5. Additional Comment**

It should be noted that the presented WEDES system in this paper is different from a traditional Inductor-Capacitor (LC) compensation wireless power transfer (WPT) system or a traditional battery energy storage system (BESS). Instead, it is a combination of both systems. Table 2 shows a comparison between the presented WEDES system in this paper, conventional WPT system [18,21,22], and conventional BESS system [23–25]. It is shown that the presented WEDES system combines the advantages of both the conventional WPT system and BESS system, and can achieve bus voltage regulation and SOC balancing, and allow for fast and safe exchange/swapping of smaller and lighter battery modules with WPT technology.

**Table 2.** Comparison summary between the WEDES system, conventional Inductor-Capacitor (LC) compensation wireless power transfer (WPT) system, and conventional battery energy storage system (BESS).


#### **6. Conclusions**

In this paper, the small-signal modeling of a distributed WEDES battery system was derived to analyze the steady-state stability and dynamic response of the entire system, as well as provide guidelines for the controller design of multiple interacted control loops. Based on the small-signal models and associated transfer functions, all three control loops, including battery module output voltage control loop, SOC balancing control loop, and bus voltage control loop with compensators, were evaluated and validated by both simulations and a 3-module WEDES battery system. It was shown that the experimental results from simulation and hardware prototype were in good agreement, which validates the accuracy and effectiveness of the derived small-signal model and designed compensators.

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

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**

1. Abu Qahouq, J.; Cao, Y. Control scheme and power electronics architecture for a wirelessly distributed and enabled battery energy storage system. *Energies* **2018**, *11*, 1887. [CrossRef]


© 2019 by the author. 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 Estimation with an Adaptive Unscented Kalman Filter Based on Model Parameter Optimization**

#### **Xiangwei Guo 1,\*, Xiaozhuo Xu 1,\*, Jiahao Geng 1, Xian Hua 2, Yan Gao <sup>1</sup> and Zhen Liu <sup>1</sup>**


Received: 27 August 2019; Accepted: 2 October 2019; Published: 6 October 2019

#### **Featured Application: The method of the research can be applied to the estimation of the remaining battery charge of electric vehicles.**

**Abstract:** State of charge (SOC) estimation is generally acknowledged to be one of the most important functions of the battery management system (BMS) and is thus widely studied in academia and industry. Based on an accurate SOC estimation, the BMS can optimize energy efficiency and protect the battery from being over-charged or over-discharged. The accurate online estimation of the SOC is studied in this paper. First, it is proved that the second-order resistance capacitance (RC) model is the most suitable equivalent circuit model compared with the Thevenin and multi-order models. The second-order RC equivalent circuit model is established, and the model parameters are identified. Second, the reasonable optimization of model parameters is studied, and a reasonable optimization method is proposed to improve the accuracy of SOC estimation. Finally, the SOC is estimated online based on the adaptive unscented Kalman filter (AUKF) with optimized model parameters, and the results are compared with the results of an estimation based on pre-optimization model parameters. Simulation experiments show that, without affecting the convergence of the initial error of the AUKF, the model after parameter optimization has a higher online SOC estimation accuracy.

**Keywords:** SOC; second-order RC model; model parameter optimization; AUKF

#### **1. Introduction**

The power battery State of charge (SOC) is generally defined as the ratio between the available capacity and the reference capacity. An accurate SOC estimation is an important basis for the formulation of an optimal energy management strategy for the whole vehicle control system of an electric vehicle (EV). It is of great significance for extending the life of the battery pack and improving the safety of the battery system [1,2]. Due to the nonlinear nature of the power battery, the SOC cannot be directly acquired by sensors; rather, it must be estimated by measuring physical quantities, such as the battery voltage, operating current, and internal resistance of the battery, and by using certain mathematical methods [3,4].

Commonly used estimation methods include the open-circuit voltage method, the ampere-hour integration method, the neural network method and the Kalman filter method [5–9]. The open-circuit voltage (OCV) can represent the discharge capacity of a lithium battery in its current state. It has a good linear relationship with the SOC when the SOC is greater than 0.1. The OCV cannot be directly measured during the working state of the battery; it can only be approximated when the battery is not working. Therefore, this method is only applicable to EVs that are parked, and the OCV method is used to provide an initial value of the SOC for other estimation methods. The ampere-hour (Ah) integral method is a relatively common, simple and reliable method for estimating the SOC. When the

discharge current is positive and the charging current is negative, the calculation formula can be expressed as:

$$SOC\_t = SOC\_0 - \frac{1}{C\_N} \int\_0^t \eta i dt \tag{1}$$

where *SOC*<sup>0</sup> is the initial SOC value, *CN* is the maximum available capacity of the battery, which is almost invariant if the time scale is small and battery ageing is accordingly ignored [8], and η is the coulomb efficiency. The ampere-hour integral method has the advantages of low cost and convenient measurement. However, there are several problems in the application of this method in EVs: (1) Other methods are needed to obtain the initial value of the SOC. (2) The accuracy of the current measurement has a decisive influence on the accuracy of the SOC estimation. (3) The cumulative error of the integration process cannot be eliminated, and if the charging or discharging time is too long during one calculation, cumulative errors can cause estimates to be unreliable.

A neural network (NN) is an intelligent mathematical tool [10,11]. A NN has the adaptability and self-learning skills to demonstrate a complex nonlinear model. NNs use trained data to estimate the SOC without knowing any information about the internal structure of the battery and the initial SOC. Theoretically, the nonlinear characteristics of the power battery can be better mapped by a NN. The advantage of this method is that it is capable of working in nonlinear battery conditions while the battery is charging/discharging. Nevertheless, the algorithm needs to store a large amount of data for training, which not only requires large memory storage but also overloads the entire system.

The Kalman filter method [12–14] is the current research hotspot for SOC estimation. The core idea of the Kalman filter is the optimal estimate in the sense of minimum variance, including the two stages of prediction and updating. In the prediction stage, the filter applies the value of the previous state to estimate the current state, and in the updating stage, the filter optimizes the predicted value in the prediction stage by using the observed value in the current state to obtain a more accurate estimate of the current state. It should be noted that the basic Kalman filter is mainly applied to linear systems. However, the estimation of the SOC is related to many factors, such as charge and discharge current, and cut-off voltage, and the influence of these factors on the SOC is nonlinear. For this reason, some people have used the extended Kalman filter (EKF) to estimate the battery SOC. Although some achievements have been made, there is a linearization error in actual use, and the Jacobian matrix is difficult to estimate. In recent years, a novel derivative Kalman filter algorithm, the unscented Kalman filter (UKF), has emerged to realize nonlinear filtering. The UKF directly uses nonlinear unscented transform techniques without the need for Taylor approximations of nonlinear equations. This process allows the mean and variance of the state of the nonlinear system to propagate directly according to the nonlinear mapping, thereby achieving a higher estimation accuracy.

In the traditional UKF algorithm [15–19], the covariance is a constant and cannot satisfy the real-time dynamic characteristics of the noise, which has a certain impact on the accuracy [20,21]. To eliminate this effect, the traditional UKF algorithm is improved by updating the covariance in real-time, which thus improves the accuracy of the UKF in this paper. This type of algorithm is called the adaptive unscented Kalman filter (AUKF) algorithm. On this basis, a reasonable optimization of the model parameters is studied, and a reasonable optimization method is proposed to further improve the accuracy of SOC estimation.

The structure of this paper is arranged as follows. In Section 1, the most commonly used methods for SOC estimation are introduced, and the proposed method of this paper is briefly described. In Section 2, the battery model is established, the model parameters are identified, and the model accuracy is verified. In Section 3, the AUKF based on a second-order resistance capacitance (RC) equivalent circuit model is presented. In Section 4, the basis and method of the reasonable optimization of the model parameters are put forward. In Section 5, the accuracy and initial error convergence of the AUKF before and after model parameter optimization are compared by simulation experiments, and the advantages of the proposed method are demonstrated. Finally, in Section 6, the research results of this paper are summarized, and future research directions are provided.

#### **2. Establishment of Battery Model and Parameter Identification**

The open-circuit voltage and internal resistance are the most basic components of an equivalent circuit model. All equivalent circuit models add other components on the basis of these two components to improve the model accuracy. Typical equivalent circuit models of lithium-ion batteries include the Rint model, the Thevenin model, the PNGV (Partnership for a New Generation of Vehicles) model and the multi-order RC loop model [22–25]. For the complicated polarization characteristics of a battery, some studies deduce and suggest that the models with more parallel RC networks connected in series should have a much higher accuracy; however, at the same time, memory is consumed, waste calculations are performed and bad real-time applications are caused.

#### *2.1. Establishment of Lithium-Ion Battery Model*

The measured terminal voltage response curve of the Sanyo lithium battery with a rated capacity of 2.6 AH at a certain state (SOC = 0.9 with a constant current discharge of 260 mAh at 0.5 C) at the end of discharge is shown in Figure 1.

**Figure 1.** Terminal voltage response at the end of discharge.

The polarization effect of a lithium battery is usually expressed equivalently by an RC loop and is equivalent to the first-order zero input response after the battery is discharged. Its terminal voltage can be expressed as an exponential term (*IRe* <sup>−</sup>*<sup>t</sup>* <sup>τ</sup> ). Region 2 in Figure 1 is the disappearance process of the polarization effect. Exponential curves are used to fit the single index, double index and triple index coefficients of region 2, and the fitting results are shown in Figure 2.

**Figure 2.** Coefficient fitting of exponential curves.

It can be seen from Figure 2 that the reduced chi-squared statistics of the double and triple exponential fittings are smaller than that of the single exponential fitting. That is, the random error effect of the double and triple exponential fittings are smaller than that of the single exponential fitting. At the same time, the correction decision coefficient of the double and triple exponential fittings is larger

than that of the single exponential fitting and is closer to 1, and the fitting effect is better. Therefore, the double and triple exponential fittings can more accurately reflect the polarization effect of the battery. By comparing the results of the double and triple exponential fittings, the sum of the squares of the residuals of the double exponential fitting is less than that of the triple exponential fitting, and the correction decision coefficient of the double exponential fitting is also closer to 1 than that of the triple exponential fitting. The results show that the double exponential fitting effect is better than the triple exponential fitting effect. The reason for this fitting result is that although the third-order RC loop can theoretically better reflect the dynamic characteristics of the battery, the third-order RC loop has one more RC loop than the second-order RC loop, which means that the third-order RC loop has two more unknowns than the second-order RC loop in the process of data fitting by computer. Therefore, the fitting effect of the third-order RC loop is not as good as that of the second-order RC loop, and it is proved that the second-order RC model is better than the others for a lithium battery. On this basis, considered comprehensively, the second-order RC equivalent circuit model, as shown in Figure 3, is adopted in this paper.

**Figure 3.** Second-order resistance capacitance (RC)RC battery model.

The model consists of three parts:

Voltage source: *V*oc represents the open circuit voltage of the power battery. In this paper, the influences of temperature and state of health (SOH) on the OCV are not considered, and the functional relationship between *V*oc and the battery SOC is studied under the same temperature and SOH conditions. Additionally, the effect of corrosion on the battery is neglected because the focus of this paper is on the equivalent circuit model.

Ohmic resistance: *R* is used to represent the internal resistance of the battery, which can be determined by the sudden change in voltage after the end of discharge, as shown in region 1 of Figure 1.

RC loop circuit: Two links of a resistor and a capacitor are superposed to simulate battery polarization, which is used to simulate the process of voltage stabilization after discharge. Region 2 of Figure 1 shows the change in voltage influenced by the RC loop circuit.

As shown in Figure 3, the functional relationship of the equivalent circuit model is as follows:

$$\begin{cases} E(t) = i\mathbb{R} + \mathbb{u}\_s + \mathbb{u}\_p + \mathbb{U}(t) = F(\text{SOC}(t)) \\ \qquad i = \frac{\mathbb{u}\_s}{\mathbb{R}\_s} + \mathbb{C}\_s \frac{\text{du}\_k}{\text{d}t} \\ \qquad i = \frac{\mathbb{u}\_p}{\mathbb{R}\_p} + \mathbb{C}\_p \frac{\text{du}\_p}{\text{d}t} \end{cases} \tag{2}$$

For the discretization of Equation (2), the state equation is solved as:

$$
\begin{bmatrix} u\_{s,k} \\ u\_{p,k} \end{bmatrix} = \begin{bmatrix} a\_s & 0 \\ 0 & a\_p \end{bmatrix} \begin{bmatrix} u\_{s,k-1} \\ u\_{p,k-1} \end{bmatrix} + \begin{bmatrix} b\_s \\ b\_p \end{bmatrix} \mathbf{l}\_{k-1} + \begin{bmatrix} w\_3(k) \\ w\_5(k) \end{bmatrix} \tag{3}
$$

$$\mathcal{U}\_k = E\_k - I\_k \mathcal{R} - \mathcal{U}\_{s\mathcal{k}} - \mathcal{U}\_{p\mathcal{k}} + \upsilon(k) = F(\text{SOC}\_k) - I\_k \mathcal{R} - \mathcal{U}\_{s\mathcal{k}} - \mathcal{U}\_{p\mathcal{k}} + \upsilon(k) \tag{4}$$

of which:

$$\begin{cases} \ a\_{\mathfrak{s}} = \mathfrak{e}^{\frac{-\overline{T}}{R\_{\mathfrak{e}}C\_{\mathfrak{s}}}}, \ b\_{\mathfrak{s}} = R\_{\mathfrak{s}} - R\_{\mathfrak{s}} \mathfrak{e}^{\frac{-\overline{T}}{R\_{\mathfrak{e}}C\_{\mathfrak{s}}}}\\\ a\_{\mathfrak{p}} = \mathfrak{e}^{\frac{-\overline{T}}{R\_{\mathfrak{p}}C\_{\mathfrak{p}}}}, \ b\_{\mathfrak{p}} = R\_{\mathfrak{p}} - R\_{\mathfrak{p}} \mathfrak{e}^{\frac{-\overline{T}}{R\_{\mathfrak{p}}C\_{\mathfrak{p}}}} \end{cases} \tag{5}$$

#### *2.2. Identification and Verification of Battery Model Parameters*

In the process of establishing the battery model, it is necessary to use the approximate linear relationship between the OCV and the SOC. In this section, based on the OCV–SOC relationship curve, the parameters of the second-order RC model are identified, and then the accuracy of the model is verified.

#### 2.2.1. Obtaining the OCV–SOC Relationship Curve

The OCV–SOC curve was obtained from a 18,650 cylindrical lithium-ion battery manufactured by the Sanyo Corporation of Japan with a rated capacity of 2.6 Ah and a rated voltage of 3.7 V, as before. Charge and discharge experiments were carried out on the new battery at a constant temperature. The battery test platform is shown in Figure 4.

**Figure 4.** Battery test platform.

The OCV–SOC curves are calibrated under constant current and constant capacity intermittent discharge conditions of 0.2 C, 0.3 C, 0.4 C, 0.5 C, 0.6 C, 0.75 C and 1 C. The specific calibration process of each group is as follows [3]:


The calibration experiment results are shown in Figure 5. When the SOC is greater than 10%, the curves almost coincide. This shows that the OCV–SOC curves corresponding to different discharge rates are approximate under the same temperature (25 ◦C) and health conditions. In this paper, the OCV–SOC curve under the condition of 0.2 C constant current intermittent discharge is selected as the reference curve.

**Figure 5.** The calibration curve results.

Using MATLAB to fit the polynomial coefficients, the approximate linear relationship between the OCV and the SOC can be obtained, as shown in Equation (6):

$$\text{Voc} = a1 \times \text{SOC}^6 + a2 \times \text{SOC}^5 + a3 \times \text{SOC}^4 + a4 \times \text{SOC}^3 + a5 \times \text{SOC}^2 + a6 \times \text{SOC} + d7 \tag{6}$$

where *a*<sup>1</sup> = −34.72, *a*<sup>2</sup> = 120.7, *a*<sup>3</sup> = −165.9, *a*<sup>4</sup> = 114.5, *a*<sup>5</sup> = −40.9, *a*<sup>6</sup> = 7.31, and *a*<sup>7</sup> = 3.231.

#### 2.2.2. Model Parameter Identification

Figure 6 is a schematic diagram of the terminal voltage response curve at the end of lithium battery discharge.

**Figure 6.** Terminal voltage response curve at the end of battery discharge.

(*V*<sup>1</sup> − *V*0) is the process in which the voltage drop generated on the ohmic resistance inside the battery disappears at the end of discharge. Thus, the ohmic resistance of the battery can be obtained as *R* = *<sup>V</sup>*1−*V*<sup>0</sup> *<sup>I</sup>* . The battery model adopted in this paper simulates the polarization process of the battery by superposing two RC links. An RC parallel circuit with a small time constant is used to simulate the process of rapid voltage change (*V*<sup>2</sup> − *V*1) when the current of the battery changes abruptly. An RC parallel circuit with a large time constant is used to simulate the process of slow voltage change (*V*<sup>3</sup> − *V*2) when the current of the battery changes abruptly.

The battery is assumed to discharge for a period of time during (*t*<sup>0</sup> − *tr*) and then stay in a resting state for the remaining time. Here, *t*0, *td*, and *tr* are the discharge start time, the discharge stop time and the static stop time, respectively. In this process, the RC network voltages are:

$$MI\_s = \begin{cases} R\_s i(t) \left[ 1 - e^{\frac{-(t-t\_0)}{\tau\_s}} \right] & t\_0 < t \le t\_d \\\ ulI\_s(t\_d) e^{\frac{-(t-t\_d)}{\tau\_s}} & t\_d < t < t\_r \end{cases} \tag{7}$$

*t*

$$\mathcal{U}\_p = \begin{cases} R\_p i(t) \Big[ 1 - e^{-\frac{-(t-t\_0)}{\tau\_p}} \Big] & t\_0 < t \le t\_d\\ \mathcal{U}\_p(t\_d) e^{\frac{-(t-t\_d)}{\tau\_p}} & t\_d < t < t\_r \end{cases} \tag{8}$$

Let τ*<sup>s</sup>* = *RsCs* and τ*<sup>p</sup>* = *RpCp* be the time constants of the two parallel circuits. The voltage change (*V*<sup>3</sup> − *V*1) is caused by the disappearance of the polarization effect of the battery. In this process, the battery voltage output is *V* = *E* − *IRse* −*t* <sup>τ</sup>*<sup>s</sup>* − *IRpe* −*t* <sup>τ</sup>*<sup>p</sup>* , which can be simply written as *V* = *E* − *ae*−*ct* − *be*−*dt*. This form can be used to fit the coefficients of the double exponential terms with MATLAB. After calculating *a*, *b*, *c* and *d*, according to *Rs* = *<sup>a</sup> <sup>I</sup>* , *Rp* <sup>=</sup> *<sup>b</sup> <sup>I</sup>* , *Cs* <sup>=</sup> <sup>1</sup> *Rsc* , and *Cp* <sup>=</sup> <sup>1</sup> *Rpd* , the values of *Rs*, *Rp*, *Cs*, and *Cp* can be identified.

#### 2.2.3. Model Accuracy Verification

After the model parameters have been identified, the accuracy of the model needs to be verified to determine whether the model parameters can be used for SOC estimation and parameter optimization in this paper. The precision of the model is verified by referring to the hybrid pulse power characterization (HPPC) experiment [26]. The simulated initial SOC is set to 0.5. The comparison between the experimentally measured voltage and the model simulation output voltage is shown in Figure 7. To better distinguish the difference, the error is presented in Figure 8.

**Figure 7.** The voltage response comparison curve.

**Figure 8.** The voltage response error.

As seen from Figures 7 and 8, when the output voltage of the battery changes due to abrupt current changes, the output voltage of the simulation model can better track the measured voltage. The maximum error is 12 mV, which indicates that the model parameters can be used for the online SOC estimation of the Kalman filter and reasonable parameter optimization.

#### **3. The AUKF Based on the Second-Order RC Equivalent Circuit Model**

Based on the equivalent circuit model, the corresponding SOC estimation algorithm can be established. A power lithium battery is a typical nonlinear system. The nonlinear system is governed by the equations of state, and the observation equations are shown in Equation (9).

$$\begin{cases} \mathbf{x}\_k = f(\mathbf{x}\_{k-1}, \boldsymbol{\mu}\_{k-1}) + \boldsymbol{\omega}\_{k-1} \\ \boldsymbol{y}\_k = h(\mathbf{x}\_k, \boldsymbol{\mu}\_k) + \boldsymbol{v}\_k \end{cases} \tag{9}$$

where the random variables *w*<sup>k</sup> and *v*<sup>k</sup> represent the process and measurement noise, respectively. For the UKF, the iteration equation is based on a certain set of sample points, which are chosen to make their mean value and variance consistent with the mean value and variance of the state variables. Then, these points will recycle the equation of the discrete-time process model to produce a set of predicted points. After, the mean value and the variance of the predicted points will be calculated to modify the results, and the mean value and the variance will be estimated. Before the UKF recursion, the state variables must be modified in a superposition of the process noise and the measurement noise of the original states. The SOC of the lithium battery pack can be calculated using the ampere-hour integral method:

$$\text{SOC}(t) = \text{SOC}(t') - \frac{1}{\text{C}\_{N}} \int\_{t'}^{t} \eta \dot{u} dt \tag{10}$$

In Equation (10), η = *kikt kc* , where *ki* is the compensation coefficient of the charge and discharge rate, *kt* is the temperature compensation coefficient, *kc* is the cycle compensation coefficient, and *CN* is the actual available battery capacity. From the equivalent circuit model shown in Figure 3, the lithium battery equation of state can be obtained from Equations (2) and (10):

$$
\begin{bmatrix} SOC\_k \\ \
u\_{s,k} \\ \
u\_{p,k} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & a\_s & 0 \\ 0 & 0 & a\_p \end{bmatrix} \begin{bmatrix} SOC\_{k-1} \\ \
u\_{s,k-1} \\ \
u\_{p,k-1} \end{bmatrix} + \begin{bmatrix} \frac{\eta T}{C\_N} \\ \
b\_s \\ \
b\_p \end{bmatrix} \mathbb{I}\_{k-1} + \begin{bmatrix} \omega\_1(k) \\ \omega\_3(k) \\ \omega\_5(k) \end{bmatrix} \tag{11}
$$

$$\begin{aligned} lL\_k &= E\_k - I\_k R - lL\_{s,k} - lL\_{p,k} + v(k) \\ &= F(SOC\_k) - I\_k R - lL\_{s,k} - lL\_{p,k} + v(k) \end{aligned} \tag{12}$$

where the values of each coefficient are shown in Equation (5). For the circuit model shown in Equations (11) and (12):

$$\begin{bmatrix} X\_k = \begin{bmatrix} x\_k^T & u\_k^T & v\_k^T \end{bmatrix}^T = \begin{bmatrix} \text{SOC}\_k & u\_{s,k} & u\_{p,k} & u\_{\text{suc},k} & u\_{s,k} & u\_{p,k} & v\_k \end{bmatrix}^T \tag{13}$$

To facilitate the distinction, we can take *xk* = *SOCk*, *Us*,*k*, *Up*,*<sup>k</sup>* as the initial state of the system, *yk* as the raw output (its corresponding symbol is *Uk* in the circuit model of the lithium battery), and *uk* as the control variable (its corresponding symbol is *Ik*), and we can make Ψ = [*y*1, *y*2, ··· , *yk*]. The operations of an ordinary UKF are as follows [27–29]:

(1) Time update of state estimation

The mean and variance of the extended state are obtained based on the optimal state estimation at the previous time. Select (2*L* + 1) sampling points (*L* is the dimension of the extended state, where *L* = 7). Finally, the sampling points are transformed by the state equation, and the state prediction is completed.

(1) Initialization, initial state determination:

$$\begin{cases}
\dot{\mathfrak{x}}\_0 = E[\mathfrak{x}\_0] \\
 P\_0 = E\left[ (\mathfrak{x}\_0 - \mathfrak{x}\_0) (\mathfrak{x}\_0 - \mathfrak{x}\_0)^T \right]
\end{cases}
\tag{14}$$

(2) State expansion:

⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩

$$\begin{aligned} \dot{X}\_0 &= E[\mathbf{X}\_0] = E[\mathbf{\hat{x}}\_0 \ 0 \ 0 \ 0] \\ P\_{\mathbf{X},0} &= E\left[ \begin{pmatrix} \mathbf{X}\_0 - \mathbf{\hat{X}}\_0 \end{pmatrix} \begin{pmatrix} \mathbf{X}\_0 - \mathbf{\hat{X}}\_0 \end{pmatrix}^T \right] = \begin{bmatrix} P\_0 \\ & Q \\ & R \end{bmatrix} \end{aligned} \tag{15}$$

$$E[w\_{m\prime} \le\_n] = \begin{cases} \ Q \ m = n \\ \ 0 \ m \ne n \end{cases} \tag{16}$$

$$E[\upsilon\_{m\_\prime}, \upsilon\_n] = \begin{cases} \ R \ m = n \\ 0 \ m \neq n \end{cases} \tag{17}$$

*Q* and *R* are covariance matrices, which are symmetric, and the diagonal is the variance in each dimension.

State mean expansion:

$$\boldsymbol{\mathfrak{X}}\_{k-1} = \begin{bmatrix} (\boldsymbol{\mathfrak{x}}\_{k-1})^T, \boldsymbol{\overline{w}}\_{k'}^T \boldsymbol{\overline{w}}\_k^T \end{bmatrix}^T \tag{18}$$

State variance expansion:

$$P\_{X,k-1} = \begin{bmatrix} P\_{x,k-1} & & \\ & Q & \\ & & R \end{bmatrix} \tag{19}$$

(3) Sample point selection

Sample = {*zi*, *Xk*−1, *i*}, *i* = 0, 1, 2, ... , 2*L* + 1, where *Xk*−1,*<sup>i</sup>* represents the selected points, and *zi* is the corresponding weighting value. Then, select the points in the following manner:

$$\begin{cases} X\_{k-1,0} = \check{X}\_{k-1} \\ X\_{k-1,i} = \check{X}\_{k-1} + \left(\sqrt{(L+\lambda)P\_{X,k-1}}\right)\_i \\ X\_{k-1,i} = \check{X}\_{k-1} - \left(\sqrt{(L+\lambda)P\_{X,k-1}}\right)\_i \end{cases} \quad i=1\sim L \tag{20}$$

The corresponding weighting values are:

$$\begin{cases} z\_0^{(m)} = \frac{\lambda}{L+\lambda} \\ z\_0^{(c)} = \frac{\lambda}{L+\lambda} + \left(1 + \alpha^2 + \beta\right) \\ z\_i^{(m)} = z\_i^{(c)} = \frac{1}{2(L+\lambda)} \quad i = 1 \sim 2L \end{cases} \tag{21}$$

where <sup>λ</sup> = <sup>α</sup>2(*<sup>L</sup>* + *<sup>t</sup>*)−*L*, *<sup>z</sup>*(*m*) is the corresponding weighting value of the mean, *<sup>z</sup>*(*c*) is the corresponding weighting value of the variance, and ( (*L* + λ)*PX*,*k*−<sup>1</sup> *i* denotes the values of column *i* in the square-root matrix (*L* + λ)*PX*,*k*−1. To ensure that the covariance matrix is definitely positive, we must take *t* ≥ 0; <sup>α</sup> controls the distance of the selected points, with 10−<sup>2</sup> <sup>≤</sup> <sup>α</sup> <sup>≤</sup> 1, and <sup>β</sup> is used to reduce the error of the higher-order terms. For a Gaussian [30], the optimal choice is β = 2 in this paper, along with *t* = 0 and <sup>α</sup> <sup>=</sup> 1. *<sup>X</sup>*<sup>ˆ</sup> *<sup>k</sup>*−1,*<sup>i</sup>* consists of *<sup>X</sup>*<sup>ˆ</sup> *<sup>x</sup> k*−1,*i* , *X*ˆ *<sup>w</sup> <sup>k</sup>*−1,*<sup>i</sup>* and *<sup>X</sup>*<sup>ˆ</sup> *<sup>v</sup> k*−1,*i* .

Based on this, the time update of the state estimation is performed:

$$\begin{array}{lcl} \mathfrak{X}\_{\mathbb{k}|\mathbb{k}-1} &= E[\{f(\mathbf{x}\_{k-1}, \boldsymbol{\mu}\_{k-1}) + \boldsymbol{\nu}\_{k-1}\} | \boldsymbol{\Psi}\_{k-1}] \\ &= \sum\_{i=0}^{2L} z\_i^{(m)} \Big[ A\_{k-1} \boldsymbol{X}\_{k-1,i}^{\mathbf{x}} + B\_{k-1} \boldsymbol{\mu}\_{k-1} + \boldsymbol{X}\_{k-1,i}^{\mathbf{w}} \Big] \\ &= \sum\_{i=0}^{2L} z\_i^{(m)} \boldsymbol{X}\_{k|\mathbb{k}-1,i}^{\mathbf{x}} \end{array} \tag{22}$$

(2) Time update of mean square error:

$$\begin{array}{lcl}P\_{\mathbf{x},\mathbf{k}|\mathbf{k}-1} &= E\Big[ (\mathbf{x}\_{\mathbf{k}} - \hat{\mathbf{x}}\_{\mathbf{k}|\mathbf{k}-1}) (\mathbf{x}\_{\mathbf{k}} - \hat{\mathbf{x}}\_{\mathbf{k}|\mathbf{k}-1})^T \Big] \\ &= \sum\_{i=0}^{2L} z\_i^{(c)} \Big( \mathbf{X}\_{\mathbf{k}|\mathbf{k}-1,i}^{\mathbf{x}} - \hat{\mathbf{x}}\_{\mathbf{k}|\mathbf{k}-1} \Big) \Big( \mathbf{X}\_{\mathbf{k}|\mathbf{k}-1,i}^{\mathbf{x}} - \hat{\mathbf{x}}\_{\mathbf{k}|\mathbf{k}-1} \Big)^T \end{array} \tag{23}$$

(3) A priori estimate of system output:

$$\begin{array}{lcl}\mathcal{Y}\_{k} &= E\{ [h(\mathbf{x}\_{k'}, \boldsymbol{\mu}\_{k}) + \boldsymbol{\upsilon}\_{k}] | \boldsymbol{\Psi}\_{k-1} \} \\ &= \sum\_{i=0}^{2L} z\_{i}^{(m)} \Big[ h\{ X\_{k|k-1'}^{x}, \boldsymbol{\mu}\_{k} \} + X\_{k-1,i}^{\upsilon} \Big] \\ &= \sum\_{i=0}^{2L} z\_{i}^{(m)} y\_{k|k-1} \end{array} \tag{24}$$

(4) Calculation of the filter gain matrix:

*Lk* = *Pxy*,*kP*<sup>−</sup><sup>1</sup> *y*,*k* <sup>=</sup> 2*<sup>L</sup> i*=0 *z* (*c*) *i* - *Xx* <sup>k</sup>|*k*−1,*<sup>i</sup>* <sup>−</sup> *<sup>x</sup>*ˆ*k*|*k*−<sup>1</sup> *yk*|*k*−1,*<sup>i</sup>* − *y*ˆ*<sup>k</sup> T* 2*L i*=0 *z* (*c*) *i* - *yk*|*k*−1,*<sup>i</sup>* − *y*ˆ*<sup>k</sup> yk*|*k*−1,*<sup>i</sup>* − *y*ˆ*<sup>k</sup> T* −<sup>1</sup> (25)

(5) Optimal state estimation:

$$\mathfrak{a}\mathfrak{x}\_{k} = \mathfrak{x}\_{k|k-1} + L\_{k}(\mathfrak{y}\_{k} - \mathfrak{y}\_{k}) \tag{26}$$

where *yk* is the actual measured value of the system output.

(6) Estimation of mean square error:

$$P\_{x,k} = P\_{x,k|k-1} - L\_k P\_{y,k} L\_k^T \tag{27}$$

As the process noise and measurement noise are measured in real-time, and to update the covariance of the process noise and measurement noise in real time, the following is needed:

$$\begin{cases} \mu\_k = y\_k - h[\pounds\_k, \mu\_k] \\ F\_k = \mu\_k \mu\_k^T \\ R\_k^v = \left( F\_k + \sum\_{i=0}^{2L} z\_i^c (y\_{k|k-1,i} - y\_k)(y\_{k|k-1,i} - y\_k)^T \right) / 2 \\ R\_k^w = L\_k F\_k L\_k^T \end{cases} \tag{28}$$

where μ*<sup>k</sup>* is the residual error of the system measured output, and *yk*|*k*−1,*<sup>i</sup>* is the residual error of the system measured output estimated by the sigma points. Real-time updates of the process noise and measurement noise can be achieved by Equation (28). Thus, the establishment of the AUKF is completed.

#### **4. Reasonable Optimization of Model Parameters**

In the process of identifying the parameters of the battery model, it was found that the amount of experimental data is always limited for a battery that works continuously for a period of time. When these raw data are substituted into the model for SOC estimation, the parameter values before and after each measured model parameter are prone to abrupt changes.

This paper reasonably optimizes the model parameters when the SOC is 0.1 *X* (where *X* is an integer from 1 to 9). In Figure 9, the ohmic resistance *R*<sup>0</sup> at the time when the experimentally measured SOC is equal to an integer multiple of 0.1 is shown.

When these raw data are substituted into the model for SOC estimation, there will be more abrupt changes or "spikes" in the variation curve of the parameters when the SOC is equal to an integer multiple of 0.1, as shown by the black circle in Figure 9. However, in actual continuous operation conditions, the change trend of each parameter of the battery is continuous and relatively flat, without an obvious inflection point. To make the data curve as flat as possible on the premise of approaching the real value, reasonable methods are considered to optimize the data measured in the test to eliminate the spike. Conventional polynomial fitting results in a very smooth fitting curve, but often has a large deviation from the true value, and the error trend at both ends of the curve is the most serious.

**Figure 9.** Original data curve.

Given the data accuracy and the slowly changing curve characteristics, this paper adopts a reasonable optimization method to weaken data mutation. In time, *SOC* = 0.1*X*(*X* = 1, 2 ...... 9), the model parameters are the identified real values, and the model parameters are optimized before and after 0.1 *X*. That is, point (0.1*X* − 0.01) and point (0.1*X* + 0.01) are selected around point *SOC* = 0.1*X* to optimize the model parameters of the point *SOC* = 0.1*X*. When *SOC* = (0.1*X* − 0.01), the model parameter values at times *SOC* = 0.1*X* and *SOC* = (0.1*X* − 0.1) are connected by a straight line, and the model parameter value at time *SOC* = (0.1*X* − 0.01) is obtained by the functional relationship of the straight line and set as *m*. When *SOC* = (0.1*X* + 0.01), the model parameter values when *SOC* = 0.1*X* and *SOC* = (0.1*X* + 0.1) are connected by a straight line, and the model parameter value at time *SOC* = (0.1*X* + 0.01) is obtained as *n* according to the functional relationship of the straight line. (*m* + *n*)/2 is taken as the optimal value of the model parameter at time *SOC* = 0.1*X*. After optimization, there is a small error in the model parameters at *SOC* = 0.1*X*, but, in return, the model parameters are closer to the curve of continuous changes in the actual working conditions.

The data curve of the battery internal resistance *R*0 measured at 25 ◦C and *SOC* = 0.5 is taken as an example. The *R*0 values of group *SOC* = 0.49 and group *SOC* = 0.51 are selected as model reference data before and after the mutation point, and the average value of the times *SOC* = 0.49 and *SOC* = 0.51 are calculated as the analysis values of the original "peak" point. The obtained optimized data curve is shown in Figure 10.

**Figure 10.** Optimized data curve.

From the optimized data curve, it can be seen that each data mutation and spike is weakened or smoothed, making the change trend of each parameter of the model closer to the actual operating conditions of the EVs. Similarly, the rest of the model parameter data are reasonably optimized (where the *Voc* original data curve is relatively smooth and will not be optimized).

#### **5. Comparison of Online SOC Estimation before and after Model Parameter Optimization**

The verification of the AUKF based on the optimization of the model parameters is divided into three aspects. First, the estimation accuracy of the AUKF is compared to that of the normal UKF. Second, the superiority of the AUKF estimation accuracy before and after the optimization of the model parameters is verified. Finally, the convergence of the AUKF to the initial value error of the SOC before and after the optimization of the model parameters is verified.

In the process of setting the battery charge and discharge state, the UDDS (Urban Dynamometer Driving Schedule) waveform [27] is used as a reference. To ensure correspondence with the experimental object (SANYO 18650 Li-ion battery) in the process of OCV–SOC calibration, the current signal in Figure 11 is adopted to describe the increase or decrease of the current in the discharging or charging process of the power battery. In one period, the average output current is 1.77 A, the maximum discharging current is 5.28 A, and the maximum charging current is 2.42 A. Each period is 1367 s, and the condition lasts for two periods.

**Figure 11.** The input current waveform.

#### *5.1. Comparison of the Estimation Accuracy between the AUKF and the UKF*

The third section explains in detail the process of establishing the AUKF. This section verifies the superiority of the AUKF and the UKF in terms of estimation accuracy in MATLAB/Simulink. The verification process is based on the model parameters before optimization. In this simulation model, the input current *I*(*k*) is integrated using the ampere-hour integral method. As there is no error in the current measurement due to outside disturbances, no accumulative error exists. Thus, the integration of the current in the simulation model can be regarded as the theoretical value of the SOC. The model simulation results and errors are shown in Figures 12 and 13, respectively.

**Figure 12.** The model simulation results.

**Figure 13.** The model simulation errors.

From the simulation results, the AUKF and the UKF can estimate the SOC online well under different charge and discharge rates. The UKF maximum estimation error is 2.98%, and the maximum estimated error of the AUKF is 2.28%; compared with the UKF, its estimation accuracy is significantly improved.

#### *5.2. Comparison of the AUKF Estimation Accuracy before and after Optimization of the Model Parameters*

Substituting the optimized model parameters into the AUKF simulation model, the results are compared with the estimated values before the optimization of the model parameters. The simulation output is shown in Figure 14, and the error is shown in Figure 15.

**Figure 14.** SOC estimation before and after optimization of the model parameters.

**Figure 15.** SOC estimation error before and after optimization of the model parameters.

From the simulation results shown in Figures 14 and 15, the online estimation of the SOC is closer to the theoretical value after the reasonable optimization of the model parameters. The maximum SOC estimation error obtained after optimization of the model parameters is 1.22% and compared to the maximum error of 2.28% before parameter optimization, the error decreases significantly. It is verified that the optimization process can significantly improve the online estimation accuracy of the AUKF based on the original data.

#### *5.3. Convergence Comparison of the Initial Error of the AUKF before and after Parameter Optimization*

The Kalman filter has a strong convergence effect on the initial error of the SOC. This section verifies the effect of the reasonable optimization of the model parameters on the initial error convergence ability of the AUKF. The initial value of the system state SOC in the AUKF is set to the same value, and the output of the AUKF estimation before and after model parameter optimization is compared. The initial SOC value is 0.9, and the simulation lasts for one cycle. A waveform comparison between the SOC value of the simulation model and the theoretical value is shown in Figure 16. Figure 17 shows the convergence of the initial error before and after the model parameter optimization of the first 400 s. The SOC estimation values before and after optimization at 200 s are given.

**Figure 16.** The comparison of the SOC estimation between the simulated value and the theoretical value.

**Figure 17.** The comparison of the estimation values in the first 400 s.

As seen from Figure 16, before and after model parameter optimization, the SOC values estimated by the AUKF converge to the theoretical value quickly. In Figure 17, at 200 s, the estimated value before optimization is 0.967, which is closer to the theoretical value than the optimized estimated value of 0.964, but the difference is only 0.003; that is, the difference is only 0.3%, which can be ignored. This shows that the convergence ability of the AUKF algorithm to the initial error before and after optimization of the model parameters is approximately the same, and optimization of the model parameters has no effect on the convergence of the algorithm.

#### **6. Conclusions**

The SOC provides vital state information for EVs and is strongly nonlinear and time-varying. For the purpose of the accurate online estimation of the SOC, the parameters of a second-order RC equivalent circuit model are identified by combining circuit principles and optimized rationally. The AUKF with optimized model parameters is used to estimate the SOC. The conclusions of this work can be summarized as follows:


In the future, to further improve the accuracy of SOC estimation, the influence of temperature and SOH on SOC estimation must be considered.

**Author Contributions:** Conceptualization, X.G. and X.X.; methodology, X.G.; software, X.G.; validation, J.G., Y.G. and X.X.; formal analysis, Z.L.; investigation, X.G.; resources, X.H.; data curation, X.H.; writing—original draft preparation, X.G.; writing—review and editing, X.X. and J.G.; visualization, X.X.; supervision, X.G.; project administration, X.G. and X.X.; funding acquisition, X.G. and X.X.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 61703145. The Key Scientific Research Projects of Higher Education Institutions in Henan Province, grant number 19A470001, 18A470014. The Key Scientific and Technological Research Projects in Henan Province, grant number 192102210073. The APC was funded by 192102210073.

**Acknowledgments:** I would like to express my deepest gratitude to my supervisor, Longyun Kang, who has provided me with valuable guidance at every stage of writing this paper. I would also like to thank the anonymous reviewers for dedicating the time to review my paper despite their busy schedules.

**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* **Cooperative Optimization of Electric Vehicles and Renewable Energy Resources in a Regional Multi-Microgrid System**

#### **Jin Chen, Changsong Chen \* and Shanxu Duan**

State Key Laboratory of Advanced Electromagnetic Engineering and Technology, School of Electrical and Electronic Engineering, Huazhong University of Science and Technology, Wuhan 430074, China; eufoniuso@163.com (J.C.); duanshanxu@hust.edu.cn (S.D.)

**\*** Correspondence: ccsfm@163.com

Received: 18 April 2019; Accepted: 27 May 2019; Published: 31 May 2019

**Featured Application: The main goal of this work is to provide a cooperative optimization method for regional multi-microgrid system, optimizing dispatching strategy and capacity allocation of electric vehicles and renewable energy sources. Meanwhile, across-time-and-space energy transmission of electric vehicles is considered in the optimization model.**

**Abstract:** By integrating renewable energy sources (RESs) with electric vehicles (EVs) in microgrids, we are able to reduce carbon emissions as well as alleviate the dependence on fossil fuels. In order to improve the economy of an integrated system and fully exploit the potentiality of EVs' mobile energy storage while achieving a reasonable configuration of RESs, a cooperative optimization method is proposed to cooperatively optimize the economic dispatching and capacity allocation of both RESs and EVs in the context of a regional multi-microgrid system. An across-time-and-space energy transmission (ATSET) of the EVs was considered, and the impact of ATSET of EVs on economic dispatching and capacity allocation of multi-microgrid system was analyzed. In order to overcome the difficulty of finding the global optimum of the non-smooth total cost function, an improved particle swarm optimization (IPSO) algorithm was used to solve the cooperative optimization problem. Case studies were performed, and the simulation results show that the proposed cooperative optimization method can significantly decrease the total cost of a multi-microgrid system.

**Keywords:** electric vehicles; renewable energy sources; microgrid; economic dispatching; capacity allocation; cooperative optimization

#### **1. Introduction**

Due to the characteristics of low carbon emission and sustainability, renewable energy sources (RESs) have attracted much attention in recent years [1–3]. In order to cope with the intermittency and uncertainty of RESs, the concept of microgrid (MG) was proposed to increase the utility of RESs by integrating them with energy storage units [4]. However, available storage units are expensive and a massive usage of them would significantly add to the cost of operation. The complexity of a microgrid's controlling is also an issue which challenges the traditional ways a power system operates.

However, increasing of electric vehicles' (EV) penetration brings more uncertainty to the operation of a power system. Uncontrolled charging behaviors of EV owners in load peak hours aggravate the burden of a power grid, posing great threats to the grid's operation.

Since EVs have the ability of storing energy, they can serve as storage units in RES-equipped power systems. Taking EVs as energy storage units and integrating them with RESs can attenuate the burden of disorderly charging while reducing the installation of expensive storage batteries [5]. Many studies focusing on the integration of EVs and RESs have been published in recent years; the main topics of these can be briefly separated into two categories, economic dispatching and system capacity allocation.

Most research has discussed the issue from the aspect of economic dispatching. The majority of the studies elaborated the optimal charging management of EVs [6,7]. However, as a kind of storage units, EVs can also discharge to the main grid or microgrid, providing a vehicle-to-grid (V2G) service to participate into the operation of the power systems [8,9]. Dispatching EVs which provide V2G services was also discussed in many papers. In [10], a MILP (Mixed-integer linear programming) model is proposed to find the optimal charging and discharging strategy of hybrid electric vehicles (HEVs). The proposed method in the paper improved the automation of HEV navigation systems and find a way to a more economic and sustainable transportation system. In [11], an integrated microgrid system considering both renewable energy source and electric vehicles is constructed. EV owners are regarded as a kind of flexible load of the demand response (DR), and the proposed optimization method significantly reduced the cost of the system. Though there is a large number of works that focus on the V2G service of EVs, some of them achieving significant progress, there are still some issues to be addressed. Most models proposed in the works are applied on only single microgrids. Interaction among different microgrids is not considered. On the other hand, the mobility of EVs is not fully exploited, and the across-time-and-space energy transmission of EVs in not considered.

Capacity allocation of the integration system is also an important issue in the optimization of the system, but not frequently mentioned in research. Existing works focus only on the planning of BESSs (Battery Energy Storage Systems) and RESs. Most of the research allocated either charging devices or RESs rather than cooperatively optimizing them at the same time. In [12], planning different types of plug-in electric vehicles (PEV) that charge infrastructures were studied. On the other hand, in [13], capacities of RESs were optimized to minimize the whole system's investment cost while meeting the load demand. In [14], an algorithm for microgrid planning was proposed considering massive connection of EV charging demands, and the system investment cost as well as CO2 emission is reduced. However, on the one hand, the relationship between economic dispatching and capacity configuration of the microgrid system is not considered in the papers. On the other hand, only a few of the studies discuss the allocation for both RESs and EVs, and the concept of cooperative optimization is not found in any research.

To sum up, there has already been much research on the topic of optimal dispatching and system allocation of RES-EV integrated systems [15,16]. Nevertheless, EV's ability of across-time-and-space energy transmission (ATSET) in the context of multi-microgrid systems is not discussed in the literature. The relationship between system economic dispatching and capacity configuration is not considered, and cooperative optimization method considering both allocation and dispatching are not mentioned in the studies. On the one hand, an allocation without system dispatching would cause a redundant installation of RESs and increase the total cost of the system. On the other hand, EV dispatching without optimal allocation of the system would decrease the efficiency of their integration into the microgrid system. Therefore, a cooperative optimization that simultaneously considers allocation and dispatching is necessary for the system's economic operation. For the reference, an optimization method using both the configuration of RESs and EVs is not considered in most of the works.

In order to fully exploit the potential of EV's mobile storage ability as well as reduce the redundant installation of RESs, a cooperative optimization considering EV's across-time-and-space energy transmission is presented in this paper. Both the installation capacity of RESs and the number of EV charging/discharging infrastructures (EVCDIs) are considered as allocation optimization variables.

The main contributions of this paper are listed as follow:


3. To fully exploiting the potential of the EV's mobile energy storage ability, both the installation capacity of RESs and the number of EVCDIs are optimized in the cooperative optimization process.

The following manuscript is organized as follows. In Section 2, basic concept of cooperative optimization is illustrated. Mathematical model of a typical regional multi-microgrid system is constructed and elaborated in Section 3. In Section 4, objective functions and restrictions are described. The cooperative optimization model is demonstrated in Section 5 and case studies are performed in Section 6. Finally, in Section 7, conclusions are drawn.

#### **2. Concept of Cooperative Optimization**

The theoretical structure of the cooperative optimization method proposed in this paper can be illustrated as in Figure 1.

**Figure 1.** Brief process of cooperative optimization of EVs (Electric Vehicles) and RESs (Renewable Energy Sources).

The model can be separated into two sections: one is the forecasting of load demands and RES generations—the accuracy of the forecast can significantly influence the availability of the optimization results—and the other is the main part of the cooperative optimization. In the optimization process, forecasting data are derived according to the mathematical character of EVs and RESs, sent to the dispatching and allocation sections. The dispatching model will dispatch the operation of microgrid system including EVs and RESs, and the allocation section can optimize the installation capacities of EVBCDIs (Electric Vehicle bidirectional charging/discharging infrastructure) and RESs. An allocation module generates the initial system allocation. A dispatching model then optimizes the operation of the system according to the system allocation and feeds the results back to the allocation module. The system allocation is then updated by an allocation module according to the data from the dispatching module. Through several repetitions of this interaction, both allocations and dispatching of the system are finally optimized. Details of the cooperative optimization process will be discussed in following sections.

In this paper, three respects of cooperative optimization are considered.

Firstly, RESs and EVs are cooperatively optimized. By optimally dispatching EVs, energy generated by RESs can be redistributed across different times and spaces, and the utilization efficiency of RESs can be improved. Furthermore, EVs can also profit from interacting with RESs and participating into the operation of the system.

Secondly, microgrids in the regional multi-microgrid system are cooperatively optimized. With the ATSET of EVs, energy can be transmitted through different microgrids, and microgrids with lower electricity prices indirectly sell their energy to microgrids with higher electricity prices. Through such a cooperation, microgrids with lower prices can sell energy to make profits, and microgrids with higher prices can purchase electricity form EVs with lower prices.

Thirdly, economic dispatching and system allocation are cooperatively optimized. An economic dispatching of the system can reduce the cost of a redundant installation of RESs and EVCDIs, and an appropriate system allocation can in turn better exploit the potential of generation and energy storage units to serve the load with less costs.

#### **3. Mathematical Model**

#### *3.1. Structure of the Multi-Microgrid System*

The structure of a classic multi-microgrid system is demonstrated as follows in Figure 2. The system is constructed with several independent microgrids, which can be briefly categorized into residential microgrids (RMGs) and office building microgrids (OBMGs). The load curves of these two kinds of microgrids are slightly different, and OBMGs have a heavier load demand in general.

**Figure 2.** Structure of a classic multi-microgrid system.

In every single microgrid, photovoltaic (PV) generation modules and domestic small-size wind turbines (WTs) are installed according to the microgrid's load demand. To satisfy the requests for charging of both residents and staff, charging plots were installed in all microgrids, and parts of them include EVBCDIs for participating in the power system's operation. For the sake of power balance and the stability of system's operation, all microgrids are connected to the main grid through connection transformers.

#### *3.2. Renewable Energy Sources*

Most load demand of microgrid systems are supplied by distributed energy sources, the majority of which are the RESs mentioned in this paper. Since RESs' characters are highly relevant to the weather and environment, their generations are volatile and intermittent. For the sake of safety and availability, an accurate and reliable prediction of RES generation is necessary. In this paper, artificial neural networks (ANNs) are adopted as the prediction model, and two models from [17,18] are used in

the proposed optimization model. Due to the restriction of the paper length, details of the two models are not demonstrated here.

#### *3.3. Electric Vehicles*

The multi-microgrid system considered in this paper is a living–working combined system. Residents are all staff from office building areas. Thus, all the EVs considered in this paper are commuting cars. They park in the RMG to charge during the knock-off hours from 19:00 to 8:00 and park in the OBMG during working time from 9:00 to 18:00. Considering the random behaviors of the staff, the arrival and departing times of EVs were derived by Monte Carlo experiment and random variables were then added.

Energy variation of EV batteries came from two aspects. One is the energy consumption of EVs for their driving on the commuting route, and the other is charging and discharging through EVBCDIs. State of capacity (SOC) was used to present the energy remaining in the mobile batteries. During the driving distance, the SOC of PEVs could be calculated as (1).

$$\text{SOC}\_{i,j}^{EV}(\mathbf{t}) = \text{SOC}\_{i,j}^{EV}(\mathbf{t} - d\_{i,j}) - D\_{i,j}^{EV} \cdot \mathbb{C}\_d \tag{1}$$

where *SOCEV <sup>i</sup>*,*<sup>j</sup>* (*t*) is the SOC of the *j*th EV in the *i*th microgrid in the *t*th hour, and *di*,*<sup>j</sup>* is the driving time. *Cd* is the energy consumption rate of EV's driving, and *DEV <sup>i</sup>*,*<sup>j</sup>* is the driving distance. Since EV's driving distances are slightly different due to each owner's random behaviors, the distance was also generated from Monte Carlo experiment, and can be demonstrated as Equation (2) [19].

$$D\_{i,j}^{EV} = E(D\_{i,j}^{EV}) + k\_{i,j} \cdot \sigma\_{i,j}^{EV} \tag{2}$$

where *DEV <sup>i</sup>*,*<sup>j</sup>* is the actual value of *DEV <sup>i</sup>*,*<sup>j</sup>* ; *<sup>E</sup>*(*DEV <sup>i</sup>*,*<sup>j</sup>* ) is the mathematical expectation of the driving distance; *ki*,*<sup>j</sup>* is a random parameter and is normally distributed; σ*EV <sup>i</sup>*,*<sup>j</sup>* is the standard deviation of the driving distance. Notice that according to the statistical data of automobile's driving, the driving distance is also normally distributed.

When parking in the OBMG or RMG, the SOCs of EVs were determined only by the charging and discharging behaviors. Considering the self-discharging of EV batteries, the SOC can be calculated according to Equation (3).

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *SOCEV <sup>i</sup>*,*<sup>j</sup>* (*t*) = *SOCEV <sup>i</sup>*,*<sup>j</sup>* (*t* − 1) · (1 − δ*SD*) −*PEV <sup>i</sup>*,*<sup>j</sup>* (*t*) · <sup>Δ</sup>*<sup>t</sup>* · <sup>η</sup>*<sup>C</sup> if PEV <sup>i</sup>*,*<sup>j</sup>* (*t*) < 0 *SOCEV <sup>i</sup>*,*<sup>j</sup>* (*t*) = *SOCEV <sup>i</sup>*,*<sup>j</sup>* (*t* − 1) · (1 − δ*SD*) −*PEV <sup>i</sup>*,*<sup>j</sup>* (*t*) · <sup>Δ</sup>*t*/η*<sup>D</sup> if PEV <sup>i</sup>*,*<sup>j</sup>* (*t*) > 0 (3)

where *PEV <sup>i</sup>*,*<sup>j</sup>* (*t*) is the charging and discharging power of the *j*th EV in the *i*th microgrid in the *t*th hour, while *PEV <sup>i</sup>*,*<sup>j</sup>* (*t*) <sup>&</sup>lt; 0 for charging and *PEV <sup>i</sup>*,*<sup>j</sup>* (*t*) > 0 for discharging. δ*SD* is the self-discharging rate. η*<sup>C</sup>* and η*<sup>D</sup>* are the charging and discharging efficiency rates, respectively.

#### *3.4. Operation Mechanism*

As illustrated in Section 3.1, the system is mainly formulated by EVs, RESs, EVBCDIs, and other assistant devices as connection transformers. In order to optimize the system's operation, data of SOC are collected by the intelligent meters installed in the EVBCDIs and then uploaded to the microgrid management system (MMS) and finally to the general management system (GMS) of the whole multi-microgrid system. The GMS optimizes system operation data and then send the dispatching results to the EVBCDIs to be performed. A brief description of the operation mechanism for the cooperative optimization is demonstrated as Figure 3.

**Figure 3.** Mechanism of the cooperative optimization.

Since there are many EVs in each microgrid, it is not impractical to optimize the charging and discharging powers of every EV, nor is it necessary. In this paper, a hierarchical structure of the management system is utilized. EVs are separated into several EV fleets according to the residential areas they belong to, and each fleet is consisted of all the EVs form the same RMG. Charging and discharging powers of EV fleets rather than single EVs are optimized by GMS, and the data are sent to the management systems of EV fleets. The power of single EVs is then optimized according to the data from the GMS. In this paper, only the optimal dispatching of EV fleets is discussed, and the dispatching method of single EVs can refer to [20] and is not comprehensively discussed in this work.

#### **4. Optimization Problem Formulation**

As a typical mathematical programming problem, objective function and optimization constraints are significant components in the problem construction. The objective function and constraints of the cooperative optimization studied in this paper are formulated as follows.

#### *4.1. Objective Function*

The main optimization object is to improve the economics of the system allocation and decrease the cost of the microgrid operator. The multi-microgrid system discussed in the proposed model, which consists of several microgrids, is assumed to be owned by one single operator. Therefore, minimizing the total operation cost of the whole system is regarded as the objective function. The total cost can be described by Equations (4)–(7).

$$\mathcal{C}\_{MMS} = \mathcal{C}\_{RMG} + \mathcal{C}\_{OBMG} \tag{4}$$

$$\mathbb{C}\_{\text{RMG}} = \sum\_{r=1}^{N\_{\text{RMG}}} \mathbb{C}\_{r}^{\text{RMG}} \tag{5}$$

$$\mathcal{C}\_{\text{OBMG}} = \sum\_{o=1}^{N\_{\text{OBMG}}} \mathcal{C}\_o^{\text{OBMG}} \tag{6}$$

$$N\_{\rm MC} = N\_{\rm RMG} + N\_{\rm OBRG} \tag{7}$$

where *CMMS* is total cost of the whole multi-microgrid system, *CRMG* and *COBMG* are the total costs of RMGs and OBMGs, and *CRMG <sup>r</sup>* and *COBMG <sup>o</sup>* are the costs of the *r*th RMG and the *o*th OBMG. *NMG* is the number of all microgrids in the system. *NRMG* and *NOBMG* are numbers of RMGs and OBMGs.

The cost of each microgrid includes the costs of wind generation, solar generation, EVs and the cost of exchanging power with the main grid. It is calculated by Equation (8).

$$\mathbb{C}\_{i}^{\rm MG} = \mathbb{C}\_{i}^{PV} + \mathbb{C}\_{i}^{\rm WT} + \mathbb{C}\_{i}^{EV} + \mathbb{C}\_{i}^{G} \tag{8}$$

where *CPV <sup>i</sup>* , *<sup>C</sup>WT <sup>i</sup>* , *<sup>C</sup>EV <sup>i</sup>* and *<sup>C</sup>EV <sup>i</sup>* are costs of PV generation, wind power, EV operation, and energy exchanging with the main grid. Costs of PV generation and wind power are formulated by capital costs and costs of maintenance—they are constant values. Therefore, the costs of PV and wind generation can be calculated as Equation (9).

$$\begin{cases} \mathbf{C}\_i^{PV} = \mathbf{C}\_{PV} \cdot \mathbf{N}\_i^{PV} \\ \mathbf{C}\_i^{WT} = \mathbf{C}\_{WT} \cdot \mathbf{N}\_i^{WT} \end{cases} \tag{9}$$

where *CPV* and *CWT* are costs of single PV module and wind turbine, and *NPV <sup>i</sup>* and *<sup>N</sup>WT <sup>i</sup>* are numbers of PV modules and wind turbines in the *i*th microgrid.

Since the battery will depreciate with the increasing of charging/discharging cycles, extra charging and discharging caused by V2G services can significantly raise the cost of EV owners. In order to compensate for the EV owners, the extra cost caused by V2G should be considered. Considering depreciation expense, the cost the EVs can be calculated by Equation (10).

$$\mathbf{C}\_{i}^{EV} = \mathbf{N}\_{i}^{EV} \cdot \mathbf{C}\_{EVBCDl} + \sum\_{t=1}^{T} \sum\_{j=1}^{N\_{i}^{EV}} P\_{i,j}^{EV}(t) \cdot p^{EV}(t) + \sum\_{j=1}^{N\_{i}^{EV}} \mathbf{C}\_{depi,j}^{EV} \tag{10}$$

where *NEV <sup>i</sup>* is the number of EVs in the *i*th microgrid, *T* is the optimization period, *CEVBCDI* is the cost of EVBCDI, *PEV <sup>t</sup>* is the price of EV charging and discharging, and *<sup>C</sup>EV dep* is the depreciation cost which is related to the charging/discharging cycles and is calculated according to [21].

Most of the energy in the main grid is generated from thermal plants, which produce large amount of pollutant as CO2 and SO2. In order to produce as little pollution as possible and to increase the utilization rate of RESs, energy exchanging with the main grid should be restricted. Therefore, in this paper, an environmental penalty is introduced. The costs of exchanging energy with the main grid can be calculated by Equation (11).

$$C\_i^G = \sum\_{t=1}^T p^G(t) \cdot P\_i^G(t) + \sum\_{t=1}^T C\_i^{Gpen}(t) \tag{11}$$

where *pG*(*t*) is the prices of trading energy with the main grid, *PG <sup>i</sup>* (*t*) is the exchanging power of the *<sup>i</sup>*th microgrid with the main grid in the *t*th hour, and *CGpen <sup>i</sup>* (*t*) is the environmental penalty of the *<sup>t</sup>*th hour.

$$C\_i^{\text{Gper}}(t) = \left[ \max(0, P\_i^G(t)) \right] \cdot k\_{\text{env}} \tag{12}$$

where *kenv* is the pollution factor, and it is determined by the types and proportion of pollutants [22].

#### *4.2. Constraints*

Some restrictions exist in the real system due to the characteristic of power devices and system balance. The constraints of the optimization problem in this paper are listed as follows.

#### 4.2.1. Balance of Power Supply and Demand

The power generated and supplied by the microgrid should meet the power demand of end users in every dispatching time interval.

$$P\_i^{EV}(t) + P\_i^{WT}(t) + P\_i^{PV}(t) + P\_i^G(t) = P\_i^L(t) \tag{13}$$

where *PPV <sup>i</sup>* (*t*), *PWT <sup>i</sup>* (*t*), and *PL <sup>i</sup>* (*t*) are power of PV generation, wind generation, and load demand respectively.

#### 4.2.2. Capacity Constraint of EV Batteries

Overcharge and discharge influence the life span of the batteries. However, in order to guarantee the driving demand of EV owners, SOC of EVs should be always higher than a certain minimum value.

$$\text{SOC}\_{\text{min}}^{EV} \le \text{SOC}\_{i,j}^{EV}(t) \le \text{SOC}\_{\text{max}}^{EV} \tag{14}$$

#### 4.2.3. Power Limits

Due to the physical limitations of EVBCDIs and connection transformers, both charging/discharging power and exchanging power are limited to a certain range.

$$\begin{cases} P\_{\text{min}}^{EV} \le P\_{\bar{i},\bar{j}}^{EV} \le P\_{\text{max}}^{EV} \\\ P\_{\text{min}}^{G} \le P\_{\bar{i},\bar{j}}^{G} \le P\_{\text{max}}^{G} \end{cases} \tag{15}$$

#### 4.2.4. Installation Constraints

To satisfy the charging demand of all EV users in the microgrid system, a total number of EVBCDIs in RMGs and the number of that in OBMGs should be equal, as demonstrated in (16) and (17).

$$\sum\_{i=1}^{N\_{MG}} N\_i^{EV} = \sum\_{r=1}^{N\_{RMG}} N\_r^{EVR} = \sum\_{o=1}^{N\_{ORAG}} N\_o^{EVO} \tag{16}$$

$$N\_{\rm MG} = N\_{\rm RMG} + N\_{\rm OBMG} \tag{17}$$

where *NRMG* and *NOBMG* are numbers of RMGs and OBMGs, and *NEVR <sup>r</sup>* and *NEVO <sup>o</sup>* are the numbers of EVBCDIs in the *r*th RMG and the *o*th OBMG.

#### **5. Two-Loop Optimization**

#### *5.1. Improved Particle Swarm Optimization*

Since the optimization problem is a non-linear problem with a complex formulation of objective function, classic mathematical programming is not a proper method to solve it. Instead, an improved particle swarm optimization algorithm (IPSO) is used in this paper. In IPSO, each potential solution is regarded as a particle, and during the optimization, the particles update their positions and velocities in each loop iteration. After several rounds of iteration, the global best value of all the particles is adopted as the optimal solution.

As a kind of intelligent algorithm, IPSO can efficiently search for the optimal solution of a complex nonlinear optimization problem. Details of the algorithm can be seen in [23,24].

#### *5.2. Two-Loop Optimization Process*

The cooperative optimization has a two-loop structure. The process is demonstrated as follows in Figure 4.

**Figure 4.** Process of the cooperative optimization. LBV: local best value, GBV: global best value.

As shown in the figure, the process consists of two loops, which correspond to optimizations of dispatching and system allocation.

The inner-loop is optimization of system dispatching. In this process, the system operation is optimized to realize the lowest total cost while the allocation of the system is fixed.

In the outer-loop, the allocation of EVs and RESs is optimized. By comparing the optimal total cost of different allocation cases, the allocation with the lowest total cost is adopted as an optimal allocation.

By employing the cooperative optimization method, economic benefits can be realized, and the system costs can be decreased.

#### **6. Case Study and Discussion**

#### *6.1. Case Description*

#### 6.1.1. Case 1

In this case, neither numbers of EV charging/discharging infrastructures nor installation capacities of renewable energy sources are optimized. EVs will be discharged within the limits of SOCs in OBMG and be charged as much as possible when they return to RMGs. Economic dispatching of a multi-microgrid system is optimized by general management system (GMS). Case 1 is adopted as a reference in comparison to other cases.

#### 6.1.2. Case 2

In this case, the numbers of EV charging/discharging infrastructures in different microgrids are optimized while the installation capacities of RESs are fixed. Similarly, economic dispatching of multi-microgrid system is optimized by GMS.

#### 6.1.3. Case 3

In this case, the numbers of EV charging/discharging infrastructures are fixed while the installation capacities of RESs are optimized. Similarly, economic dispatching of multi-microgrid system is optimized by GMS.

#### 6.1.4. Case 4

In this case, economic dispatching and capacity allocation of RESs and EVs in multi-microgrid system are optimized synergistically. EVs and RESs can cooperate best in this case, and there is no redundant installation of EVBCDIs or RESs.

#### *6.2. Simulation System Construction*

For simplification of the simulation, a multi-microgrid system with two RMGs and one OBMG is constructed. All EV owners are living in the two RMGs, identified as RMG1 and RMG2, and working in the OBMG. Exchanging power limit of the RMGs is <sup>−</sup>95 kW <sup>≤</sup> *PG <sup>i</sup>* (*t*) ≤ 95 kW, and the limit for OBMG is <sup>−</sup>250 kW <sup>≤</sup> *<sup>P</sup><sup>G</sup> <sup>i</sup>* (*t*) ≤ 250 kW. Since the three microgrids are geographically near, PV generation and wind generation curves are similar and are forecasted by artificial neural networks using historical weather data derived from local weather records. Load demands of the three microgrids are derived from data of real microgrids with similar sizes.

A time-of-use electricity price is adopted as the energy exchanging prices with the main grid. The prices are demonstrated in Figure 5 [23]. In order to guarantee the benefits of the EV owners, charging and discharging prices of EVs should be restricted, as in (18)

$$p\_C^{EV}(t) \le \eta\_C \cdot \eta\_D \cdot p\_D^{EV}(t) \tag{18}$$

where *pEV <sup>C</sup>* (*t*) and *pEV <sup>D</sup>* (*t*) are charging and discharging prices of EVs. They are assumed to be time-invariable in this paper. Additionally, to encourage participation of microgrids, discharging prices of EVs should be lower than the highest exchanging price in OBMG, and charging prices should be higher than the lowest exchanging price in RMG.

$$\begin{cases} \begin{array}{l} p\_C^{EV}(t) \ge \min\left(p\_{\text{BMG}}^G(t)\right) \\ p\_D^{EV}(t) \le \max\left(p\_{\text{OBMGG}}^G(t)\right) \end{array} \tag{19} $$

**Figure 5.** Prices of energy exchanging and EV charging and discharging.

For the EV model, BYD E6 is selected in the simulation. The capacity of the EV battery is 80 kWh and charging and discharging power limits are both 7 kW in normal charging mode. Charging and discharging efficiency are also identical, at 0.9. The daily cost of installing an EVBCDI is 0.45\$. Some other data are referred to in [25].

#### *6.3. Simulation Results*

#### 6.3.1. Case 1

In case one, the numbers of EVBCDIs and installation capacities of RESs are set to a fixed number according to the load demand and neither of them are optimized. The numbers of EVCDIs for RMG1, RMG2 and OBMG are 10, 10, and 20, respectively. WT installation capacities for RMG1, RMG2, and OBMG are 240 kW, 160 kW, and 320 kW. PV installation capacities for RMG1, RMG2, and OBMG are 180 kW, 120 kW, and 240 kW. The two-day operation profiles of RMG1 and OBMG in case one are shown in Figure 6. Profiles of RMG2 are similar to RMG1 and are not demonstrated here.

**Figure 6.** Dispatching results in Case 1: (**a**) Dispatching profiles of RMG1; (**b**) Dispatching profiles of OBMG.

As demonstrated in Figure 6, the EVs tend to charge in RMGs to satisfy their driving demands for the next day, and discharge in the OBMG to participate in the operation of the microgrid and earn some benefit by trading energy with the microgrid operator.

For most systems, the electricity prices of the main grid are always related to the load demand. The electricity prices are high in load peak hours, and the prices are low in load valley hours [26]. However, as shown in Figure 6, the microgrids tend to sell energy to the main grid when the electricity prices are high, and purchase energy from the main grid the prices of electricity are low. Therefore, the energy trading among microgrids and the grid can contribute to peak shaving of the grid.

In Case 1, the total cost of one year is \$195,195.80. The costs for RMG1, RMG2, and OBMG are \$67,872.80, \$38,063.00, and \$89,260.00, respectively.

#### 6.3.2. Case 2

In Case 2, the numbers of EVBCDIs in three microgrids are optimized, while installation capacities of RESs are kept the same as Case 1. The optimal numbers of EVCDIs for RMG1, RMG2, and OBMG are 14, 15, and 29. WTs installation capacities for RMG1, RMG2, and OBMG are 240 kW, 160 kW, and 320 kW. PV installation capacities for RMG1, RMG2, and OBMG are 180 kW, 120 kW, and 240 kW. Since the two-day operation profiles of RMG1 are identical to those in Case 1, the operation profiles of RMG2 are shown here. The two-day operation profiles of RMG2 and OBMG in Case 2 are shown in Figure 7.

Tendency of EV charging and discharging is similar with Case 1. However, the numbers of EVCDIs in both RMGs and OBMG are optimized. Comparing Figure 7 with Figure 6, it is obvious that the total charging and discharging power of EVs increased. The total energy transmitted between RMGs and OBMG increased, and more energy with lower prices are transported to the places and times with higher prices. Therefore, the operation cost of the system decreased, and the total cost also decreased.

In Case 2, total cost of one year is \$193,867.50. The costs for RMG1, RMG2, and OBMG are \$67,491.10, \$37,374.50, and \$89,001.90, respectively.

#### 6.3.3. Case 3

In Case 3, the numbers of EVBCDIs in three microgrids are fixed to be the same with Case 1, and installation capacities of RESs are optimized. The numbers of EVCDIs for RMG1, RMG2, and OBMG are 10, 10, and 20 respectively. Optimal WT installation capacities for RMG1, RMG2, and OBMG are 80 kW, 40 kW, and 360 kW. The installation capacity of PV modules for RMG1, RMG2, and OBMG are 180 kW, 0 kW, and 400 kW. The 2-day operation profiles of RMG1 and OBMG in Case 3 are shown in Figure 8.

**Figure 7.** Dispatching results in Case 2: (**a**) Dispatching profiles of RMG1; (**b**) Dispatching profiles of OBMG.

**Figure 8.** Dispatching results in Case 3: (**a**) Dispatching profiles of RMG1; (**b**) Dispatching profiles of OBMG.

According to the figures, we can see that the profiles of RMG1 are not much different. The installation capacities of RESs in RMGs and OBMG have been optimized. For the RMGs, electricity prices of the distribution grid are low and it is not beneficial to sell energy to the grid. Therefore, the redundant installation of RESs is decreased, and less energy are sold to the grid while the load demands of RMGs are still satisfied. For OBMG, the electricity prices of the distribution grid are high, and it is beneficial to sell extra energy to the grid. Therefore, more RESs are installed to sell more energy to the grid, and the microgrid earns more from selling energy.

In Case 3, the total cost of one year is \$161,090.30. The costs for RMG1, RMG2, and OBMG are \$55,138.10, \$36,348.80, and \$69,603.40, respectively.

#### 6.3.4. Case 4

In Case 4, both the numbers of EVBCDIs of the three microgrids and installation capacities of RESs are optimized. The numbers of EVCDIs for RMG1, RMG2, and OBMG are 19, 8, and 27, respectively. WTs installation capacities for RMG1, RMG2, and OBMG are 80 kW, 80 kW, and 360 kW. PV installation capacities for RMG1, RMG2, and OBMG are 120 kW, 0 kW, and 400 kW. The two-day operation profiles of RMG1 and OBMG in Case 1 are shown in Figure 9.

In this case, both EV charging devices and RESs are optimally allocated. Since the cost of PV is higher than wind generation and since the outputs of PV modules are unstable, the installation of PV modules is decreased. Meanwhile, the number of EVCDIs and the installation capacities of RESs are optimally matched. Neither redundant RES installation nor redundant installation of EVCDIs exists. Therefore, the total cost in this case is the lowest.

In Case 4, the total cost of one year is \$159,264.80. The costs for RMG1, RMG2, and OBMG are \$56,012.30, \$36,306.40, and \$66,946.10, respectively.

#### *6.4. Comparison*

Optimization results of system allocation are demonstrated in Table 1. Results of costs and the net exchanged energy between microgrids and the main grid in four cases are demonstrated in Table 2. The total costs and costs of the three microgrids in the four cases are listed. In addition, the costs of allocation and dispatching are also calculated and listed in the table. According to the figures in Section 6.3 and the data in the tables, conclusions can be drawn as follows:


**Figure 9.** Dispatching results in Case 4: (**a**) Dispatching profiles of RMG1; (**b**) Dispatching profiles of OBMG.


**Table 1.** Comparison of Allocations in all Four Cases.

**Table 2.** Comparison of Costs and Net Exchanged Energy in all Four Cases.


#### **7. Conclusions**

In this paper, a cooperative optimization method for a multi-microgrid system was proposed. Mathematical models of multi-microgrid systems including RESs and EVs were constructed. Energy cooperation of microgrids in the multi-microgrid system was accomplished via the mobile energy storage of EVs. The impact of ATSET of EVs on economic dispatching and capacity allocation of a multi-microgrid system was analyzed. A two-loop optimization methodology that considered the cooperation of economic dispatching and system allocation was proposed, and an IPSO algorithm was used to solve the optimization problem. Both EVs and RESs were succesfully optimized in the optimization. Four cases were presented to verify the cooperative optimization. By comparing the results of the four cases, the following conclusions are drawn:


It is worth noting that the simulation system discussed in this paper was significantly simplified in order to clearly demonstrate the main methodology of the paper. Future works may focus on the application of the model to more complex systems. Real systems are always more complex, with various kinds of generation and load units, and some realistic factors such as the influence of the battery's remnant capacity on the behavior of EV owners can also be considered. A real system may also include storage batteries and micro-turbines. The way of applying the proposed model on real systems is an important topic of the future works. Moreover, since the charging and discharging prices

of EVs considered in this paper are constant, a more flexible pricing system for EVs' charging and discharging is also an issue that worth studying.

**Author Contributions:** Conceptualization, J.C. and C.C.; methodology, J.C. and C.C.; software, J.C.; validation, J.C.; formal analysis, J.C. and C.C.; investigation, J.C. and C.C.; resources, J.C. and C.C.; data curation, J.C. and C.C.; writing—original draft preparation, J.C.; writing—review and editing, J.C. and C.C.; visualization, J.C.; supervision, C.C. and S.D.; project administration, C.C. and S.D.; funding acquisition, C.C.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 51477067 and the Lite-On Power Electronics Technology Research Fund, grant number PRC20161047.

**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*

## **Adaptive Dual Extended Kalman Filter Based on Variational Bayesian Approximation for Joint Estimation of Lithium-Ion Battery State of Charge and Model Parameters**

#### **Jing Hou 1,†, Yan Yang 1,\*, He He <sup>2</sup> and Tian Gao <sup>1</sup>**


Received: 4 April 2019; Accepted: 23 April 2019; Published: 26 April 2019

**Abstract:** An accurate state of charge (SOC) estimation is vital for the safe operation and efficient management of lithium-ion batteries. At present, the extended Kalman filter (EKF) can accurately estimate the SOC under the condition of a precise battery model and deterministic noise statistics. However, in practical applications, the battery characteristics change with different operating conditions and the measurement noise statistics may vary with time, resulting in nonoptimal and even unreliable estimation of SOC by EKF. To improve the SOC estimation accuracy under uncertain measurement noise statistics, a variational Bayesian approximation-based adaptive dual extended Kalman filter (VB-ADEKF) is proposed in this paper. The variational Bayesian inference is integrated with the dual EKF (DEKF) to jointly estimate the lithium-ion battery parameters and SOC. Meanwhile, the measurement noise variances are simultaneously estimated in the SOC estimation process to compensate for the model uncertainties, so that the adaptability of the proposed algorithm to dynamic changes in battery characteristics is greatly improved. A constant current discharge test, a pulse current discharge test, and an urban dynamometer driving schedule (UDDS) test are performed to verify the effectiveness and superiority of the proposed algorithm by comparison with the DEKF algorithm. The experimental results show that the proposed VB-ADEKF algorithm outperforms the traditional DEKF algorithm in terms of SOC estimation accuracy, convergence rate, and robustness.

**Keywords:** state of charge (SOC); joint estimation; lithium-ion battery; variational Bayesian approximation; dual extended Kalman filter (DEKF); measurement statistic uncertainty

#### **1. Introduction**

Electric vehicles (EVs) are believed worldwide to be one of the most important development directions in the vehicle industry because of their advantages in low pollution and energy saving. Lithium-ion batteries, by virtue of their high energy and power density, are the fundamental power source of EVs [1]. Nevertheless, lithium-ion batteries have hidden dangers in safety. Once over-charged or over-discharged, the battery capacity will drop, leading to reduced lifetime or damage, or even explosion. Moreover, the battery characteristics change with different operating conditions. Therefore, to ensure the safety, reliability, and efficiency of the battery system, a battery management system (BMS) is developed and utilised. The BMS is used to control and monitor the battery operating conditions to guarantee the safe operation and longevity of the battery.

The state of charge (SOC) is a significant parameter which indicates the amount of remaining energy in the battery for further service. It needs to be estimated accurately online in order to prevent any over-charge and over-discharge, while providing the information and support for the effective and flexible operation of the vehicle. However, it is difficult to acquire precise SOC estimates since the battery itself is a highly nonlinear system and has a lot of uncertainties.

Extensive studies have been made into SOC estimation. Hannan et al. [2] gave a comprehensive review of these research fruits. Overall, these methods can be mainly classified into four categories: (1) the open-circuit voltage (OCV) method; (2) current integration or coulomb counting (CC) method; (3) data-driven methods; and (4) model-based methods.

The OCV method is based on the one-to-one corresponding relationship of open-circuit voltage and SOC. It is simple but susceptible to temperature, battery aging, and other factors, and the acquisition of accurate measurement of OCV needs lots of rest time, making it almost impossible for moving vehicles. The CC method is widely used in many applications, but it has two weak points. First, its estimation accuracy is strongly dependent on the correctness of the initial SOC value, and second, it can easily diverge due to error accumulation in the current measurement.

The data-driven methods include neural network (NN) algorithms [3–5], fuzzy logic (FL) algorithms [6–8], support vector machines (SVMs) [9–11], and so on. For instance, He et al. [3] employed a multilayer feedforward neural network to estimate the SOC of Lithium-ion batteries by use of the discharge current, terminal voltage, and temperature of the battery as input. Zhao et al. [4] combined a back propagation neural network and adaptive Kalman filter to estimate the SOC, and used a forgetting-factor recursive least-squares (FFRLS) algorithm to identify the time-varying battery model parameter. In [6], a fuzzy logic algorithm was applied to estimate SOC of lithium-ion batteries for the application of a portable defibrillator. AC impedance and voltage recovery measurements were used as the input parameters for the FL model. In [7], a more advanced algorithm named adaptive neuro-fuzzy inference system (ANFIS) was developed to estimate SOC. In addition, Hu et al. [9] proposed an SOC estimation based on an optimized support vector machine for regression with double search optimization process. In [10], an SOC estimation method based on fuzzy least-squares SVM was proposed. These data-driven algorithms do not need to know any battery characteristics and have a good ability of nonlinear mapping and self-learning. However, they require a large amount of experimental data to train the intelligent model beforehand. Different types of training data and training methods exert a great influence on the model error. When the training data cannot cover all the operating conditions, for example, the experiments are incomplete or the battery characteristics have changed, the SOC estimation accuracy will decrease.

The model-based methods are dominated by Kalman filter and its derivatives. Chen et al. [12] used an extended Kalman filter (EKF) along with a nonlinear battery model to estimate the SOC of the lithium-ion battery. In order to further improve the accuracy, adaptive extended Kalman filter (AEKF), unscented Kalman filter (UKF), and adaptive unscented Kalman filter (AUKF) are proposed. He et al. [13] identified the parameters of an improved Thevenin battery model using EKF, and then adopted an AEKF for obtaining correct and robust SOC of the lithium-ion batteries. In [14], a UKF-based method was used to self-adjust the model parameters and provide state of charge estimation of the battery. References [15,16] adopted AUKF to realize SOC estimation. In addition, an adaptive cubature Kalman filter (CKF) was proposed in [17] to improve the convergence rate and SOC estimation accuracy. To overcome the accuracy degradation caused by non-Gaussian noise, particle filters were utilized to estimate the SOC of batteries in [18,19]. To compensate for the time-variability of battery parameters due to variational operating conditions and battery aging, a dual EKF and a dual UKF were employed in [20,21] for simultaneous SOC and parameters estimation, respectively. The model-based methods eliminate the need for correct initial SOC values and accurate measurement, which are requisite for the CC method, and have no demand for a large amount of train data, so it is widely applied and is studied extensively. However, the estimation accuracy is strongly affected by the battery model and parameters. Once there is a model mismatch, the estimation performance will rapidly decline.

In practical applications, the battery model parameters will change with SOC, temperature, and the degree of battery aging. Moreover, the statistic information of the process noise and sensor noise may be unknown or time-varying. Under these situations, the traditional KF-based methods will have low estimation accuracy and poor robustness. To address this issue, some researchers [22–24] resort to a H∞ filter, which takes the time-varying battery parameter into account and has no need to know the statistics of the process noise and the measurement noise. It has strong robustness under uncertain conditions. However, the H∞ filter is a tradeoff between robustness and accuracy, so the SOC estimation accuracy is sacrificed to some extent for the robustness. Liu et al. [25] combined the idea of square root filter and adaptive unscented Kalman filter (ASRUKF) algorithm based on improved Sage-Husa estimation to estimate the SOC of a lithium-ion battery. This method adaptively adjusted the values of the process and measurement covariances in the estimation process to improve the accuracy of SOC estimation. However, it does not consider the uncertainties brought by varying battery model parameters. EI Din et al. [26] proposed a multiple-model EKF (MM-EKF) and an autocovariance least squares (ALS) method for estimating the SOC of lithium-ion battery cells. MM-EKF reduced the dependence of the EKF algorithm on the correct assumptions of the measurement noise statistics by weighted summation of the estimates of multiple hypothesized EKFs. The ALS method extracted the possible correlation in the innovation sequence to estimate the measurement noise covariance. However, both methods leave out consideration of the time-varying battery parameters. Furthermore, the computation load is larger compared with the conventional EKF algorithm.

In fact, the Bayesian approach is the most general approach of solving the problem with uncertain parameters. But it is hard to get the analytical solution for most Bayesian approaches due to complex probability density function or high dimension of integration. Recently, the variational Bayesian (VB) inference method [27–31] has drawn extensive attentions, which utilizes a new simpler, analytically tractable distribution to approximate the true posterior distribution so as to avoid the direct complex calculation of multi-dimensional probability density function. Sarkka et al. [27] adopted the VB method for joint recursive estimation of the dynamic state and the time-varying measurement noise parameters in linear state space models. Li et al. [28] employed VB approximation for the unscented Kalman filter to estimate the time-varying measurement noise covariance so as to improve algorithm adaptability. Sun et al. [29] proposed a VB method to estimate the system states with unknown inputs. Hou et al. [30] combined the VB method with the shifted Rayleigh filter to jointly estimate the target position and the clutter probability so that improving the performance of bearings-only target tracking. In [31], the VB approach was applied to estimate the ARX model parameters along with time delays.

In this paper, we combine the idea of variational Bayesian inference with the dual EKF algorithm (VB-ADEKF) to jointly estimate the battery parameters and SOC of lithium-ion batteries of electric vehicles. Meanwhile, the measurement noise variances are simultaneously estimated with the state estimation to account for the battery model uncertainties and measurement uncertainties, minimizing the impact of model mismatch. The effectiveness of the proposed algorithm has been verified through experiments under different operating conditions. The results show that the proposed VB-ADEKF algorithm outperforms the dual EKF algorithm in terms of SOC estimation accuracy, convergence rate, and robustness. The mean SOC estimation error of VB-ADEKF is under 1% for most of the time.

The remainder of this paper is organized as follows: Section 2 describes the battery model, the definition of SOC and establishes the state space models for SOC estimation and battery parameter estimation. Section 3 illustrates the variational Bayesian approximation-based adaptive Kalman filter and Section 4 presents the proposed variational Bayesian approximation-based adaptive dual extended Kalman filter. The experimental verification and analysis are presented in Section 5. Finally, Section 6 provides a conclusion.

#### **2. Battery Modeling**

#### *2.1. Battery Model*

For the accurate estimation of SOC, a reliable battery model is required. Considering the model accuracy, the structure complexity, and the computation time, the first order resistor–capacitor (RC) model as shown in Figure 1 is adopted to model the lithium-ion battery.

**Figure 1.** A first-order resistor–capacitor (RC) battery model.

The electrical behavior of the model can be written as follows:

$$
\mathcal{U}\mathcal{U}\_l = \mathcal{U}\_{\text{0c}} - \mathcal{U}\_1 - I\_L \mathcal{R}\_0 \tag{1}
$$

$$
\dot{\mathcal{U}}\_1 = \frac{I\_L}{\mathcal{C}\_1} - \frac{\mathcal{U}\_1}{R\_1 \mathcal{C}\_1} \tag{2}
$$

where *Ut* denotes the terminal voltage of the battery, *Uoc* is the open-circuit voltage, *U*<sup>1</sup> is the polarization voltage of the RC network, *IL* is the load current, *R*<sup>0</sup> represents the ohmic internal resistance, and *R*<sup>1</sup> and *C*<sup>1</sup> represent the polarization resistance and polarization capacitance, respectively.

The nonlinear relationship between the open-circuit voltage and the SOC is described using the fifth-order polynomial model as:

$$dL\_{\rm ox}(SOC) = k\_0 + k\_1 SOC + k\_2 SOC^2 + k\_3 SOC^3 + k\_4 SOC^4 + k\_5 SOC^5 \tag{3}$$

where *k*0, *k*1, *k*2, *k*3, *k*4, *k*<sup>5</sup> are the parameters to be identified.

#### *2.2. Definition of State of Charge*

The state of charge (SOC) is defined as the ratio of the remaining capacity in a battery over the rated battery capacity. Using the CC method, the battery SOC can be calculated as follows:

$$SOC\_t = SOC\_{t\_0} - \frac{1}{Q\_{rate}} \int\_{t\_0}^t \eta\_c I\_{L,t} dt\tag{4}$$

where *η<sup>c</sup>* is the coulomb efficiency, *IL*,*<sup>t</sup>* is the load current at time *t*. *Qrate* is the rated capacity of the battery.

#### *2.3. State-Space Model*

#### 2.3.1. State Space Model for SOC Estimation

Taking *X* = [*SOC*, *U*1] *<sup>T</sup>* as the state vector, the load current *IL* as the input and the terminal voltage *Ut* as the output, we can obtain the discrete state-space model as

$$\begin{cases} X\_{k+1} = f(X\_{k'} I\_{L,k'} \theta\_k) + w\_k \\\ y\_k = \mathcal{U}\_{t,k} = h(X\_{k'} I\_{L,k'} \theta\_k) + v\_k \end{cases} \tag{5}$$

where *θ<sup>k</sup>* = [*R*0, *R*1, *τ*1] <sup>T</sup> represents the time-varying battery model parameter vector; *τ*<sup>1</sup> = *R*1*C*<sup>1</sup> is the time constant of the RC network; *wk* ∼ N (0, *<sup>Q</sup><sup>x</sup> <sup>k</sup>* ) is the Gaussian process noise with covariance *Q<sup>x</sup> <sup>k</sup>* ; *vk* ∼ N (0, <sup>Σ</sup>*<sup>x</sup> <sup>k</sup>* ) is the measurement noise with variance <sup>Σ</sup>*<sup>x</sup> <sup>k</sup>* . The initial state has a Gaussian prior distribution *<sup>X</sup>*<sup>0</sup> ∼ N (*X*<sup>ˆ</sup> 0, *<sup>P</sup>*0). The state prior and process noise are assumed to be known, while the measurement noise variance Σ*<sup>x</sup> <sup>k</sup>* is assumed to be unknown. In addition, the process noise, measurement noise, and initial state value are independent of each other.

*f*(·) and *h*(·) represent the nonlinear functions of state vector *Xk*, input *IL*,*k*, and battery model parameter vector *θk*. Their mathematical expressions are

$$f(\cdot) = \begin{bmatrix} 1 & 0 \\ 0 & e^{-\frac{\Delta t}{\eta\_1}} \end{bmatrix} \begin{bmatrix} \text{SOC}\_k \\ \text{U}\_{1,k} \end{bmatrix} + \begin{bmatrix} -\frac{\eta\_s \Delta t}{Q\_{ctr}} & 0 \\ 0 & R\_1(1 - e^{-\frac{\Delta t}{\eta\_1}}) \end{bmatrix} I\_{1,k} \tag{6}$$

$$h(\cdot) = \mathcal{U}\_{\text{OC}}(SOC\_k) - \mathcal{U}\_{1,k} - I\_{L,k}\mathcal{R}\_0 \tag{7}$$

where Δ*t* is the sampling interval of the current.

The state transition matrix and the input control matrix are respectively

$$F\_k = \begin{bmatrix} 1 & 0 \\ 0 & e^{-\frac{\Delta t}{\tau\_1}} \end{bmatrix} \tag{8}$$

$$\mathbf{G}\_k = \begin{bmatrix} -\frac{\eta\_c \Delta t}{Q\_{rate}} & 0\\ 0 & R\_1 (1 - e^{-\frac{\Delta t}{T\_1}}) \end{bmatrix} \tag{9}$$

Linearizing the measurement function, we can get the Jacobian measurement matrix as

$$H\_{\mathbf{x},k} = \frac{\partial h(\cdot)}{\partial X\_k} \bigg|\_{X\_k = \hat{X}\_k^-} = \begin{bmatrix} \frac{\partial L\_{\mathbf{OC}}}{\partial \mathbf{SC}\_k} & -\mathbf{1} \end{bmatrix} \tag{10}$$

#### 2.3.2. State Space Model for Battery Parameter Estimation

Because the battery model parameters vary with changes in the batteries' SOC, degree of aging, and environmental temperature, online recursive battery parameter estimation is needed. So here we establish the state-space equations of the battery parameters as:

$$\begin{cases} \theta\_{k+1} = \theta\_k + r\_k\\ d\_k \quad = h(\mathbf{X}\_{k'} I\_{L,k'} \theta\_k) + \varepsilon\_k \end{cases} \tag{11}$$

where *rk* ∼ N (0, *<sup>Q</sup><sup>θ</sup> <sup>k</sup>* ) is a small white noise with covariance *<sup>Q</sup><sup>θ</sup> <sup>k</sup>* that reflects the time-varying parts of the parameters, *dk* is a measurement function of *<sup>θ</sup>k*, and *ek* ∼ N (0, <sup>Σ</sup>*<sup>θ</sup> <sup>k</sup>* ) is the measurement noise to account for the sensor noise and modeling uncertainties. Σ*<sup>θ</sup> <sup>k</sup>* is assumed unknown here.

The Jacobian measurement matrix of *θ<sup>k</sup>* is calculated as follows [20]:

$$H\_{\theta,k} = \frac{d\boldsymbol{h}(\cdot)}{d\theta\_k}\bigg|\_{\theta\_k = \theta\_k^-} = \frac{\partial \boldsymbol{h}(\cdot)}{\partial \hat{\theta}\_k^-} + \frac{\partial \boldsymbol{h}(\cdot)}{\partial \hat{X}\_k^-} \cdot \frac{d\hat{X}\_k^-}{d\hat{\theta}\_k^-} \tag{12}$$

$$\frac{\partial \mathbf{h}(\cdot)}{\partial \hat{\theta}\_k^-} = \begin{bmatrix} \ -I\_{L,k} & 0 & 0 \ \ \end{bmatrix} \tag{13}$$

$$\frac{\partial h(\cdot)}{\partial \hat{X}\_k^{-}} = H\_{x,k} \tag{14}$$

*Appl. Sci.* **2019**, *9*, 1726

$$\frac{d\hat{X}\_k^-}{d\hat{\theta}\_k^-} = \begin{bmatrix} 0 & 0 & 0\\ 0 & a\_{22} & a\_{23} \end{bmatrix} \tag{15}$$

where

$$\begin{aligned} a\_{22} &= -I\_{L,k} \cdot \left( \exp(\Lambda t/\tau\_1^2) - 1 \right) \\ a\_{23} &= \left( \Delta t/\tau\_1^2 \right) \cdot \left( \pounds\_{2,k}^{-} - R\_1 I\_{L,k} \right) \cdot \exp(\Lambda t/\tau\_1) \end{aligned} \tag{16}$$

#### **3. Variational Bayesian Approximation-Based Adaptive Kalman Filter Algorithm**

As we know, when the measurement noise covariance is unknown or time varying, the classical approach of solving the problem is to use adaptive filters. The Bayesian approach is a typical choice. But it is usually hard to get the analytical form due to complex probability density function or high dimension of integration. In [27], variational Bayesian (VB) approximation was firstly used for the Kalman filter to estimate the joint posterior distribution of the states and noise covariances.

Given a discrete-time linear state space model as:

$$X\_{k+1} = F\_k X\_k + w\_k \tag{17}$$

$$y\_k = H\_k X\_k + v\_k \tag{18}$$

where *wk* ∼ N (0, *Qk*) is a Gaussian distributed process noise, *vk* ∼ N (0, Σ*k*) is a Gaussian measurement noise with diagonal covariance Σ*k*. Note *Qk* is assumed known but Σ*<sup>k</sup>* is unknown.

Assuming the state vector and the measurement noise covariance are independent, the joint posterior distribution of the state and covariance *p*(*Xk*,Σ*k*|*y*1:*k*) can be solved by VB approximation as follows:

$$p(X\_k \: \Sigma\_k | y\_{1:k}) \approx Q\_x(X\_k) Q\_{\Sigma}(\Sigma\_k) \tag{19}$$

The VB approximation can now be formed by minimizing the Kullback–Leibler (KL) divergence between the separable approximation and the true posterior:

$$\begin{split} KL\left[Q\_{\mathbf{x}}(\mathbf{X}\_{k})Q\_{\Sigma}(\Sigma\_{k})||p(\mathbf{X}\_{k},\Sigma\_{k}|y\_{1:k})\right] \\ = \int Q\_{\mathbf{x}}(\mathbf{X}\_{k})Q\_{\Sigma}(\Sigma\_{k}) \times \log\left(\frac{Q\_{\mathbf{x}}(\mathbf{X}\_{k})Q\_{\Sigma}(\Sigma\_{k})}{p(\mathbf{X}\_{k'},\Sigma\_{k}|y\_{1:k})}\right)dX\_{k}d\Sigma\_{k} \end{split} \tag{20}$$

Minimizing the KL divergence with respect to the probability densities *Qx*(*Xk*) and *Q*Σ(Σ*k*) in turn, while keeping the other fixed, we can get the following equations:

$$Q\_x(X\_k) \propto \exp\left(\int \log p(y\_{k'} | X\_{k'} \Sigma\_k | y\_{1:k-1}) Q\_\Sigma(\Sigma\_k) d\Sigma\_k\right) \tag{21}$$

$$Q\_{\Sigma}(\Sigma\_k) \propto \exp\left(\int \log p(y\_{k'}X\_{k'}\Sigma\_k | y\_{1:k-1}) Q\_x(X\_k) dX\_k\right) \tag{22}$$

Computing the above equations, we can get the following densities [27]:

$$Q\_x(X\_k) = \mathcal{N}(X\_k | m\_{k'} P\_k) \tag{23}$$

$$Q\_{\Sigma}(\Sigma\_k) = \prod \text{Inv-Gamma}(\sigma\_{k,i}^2 | a\_{k,i}, \beta\_{k,i}) \tag{24}$$

where the parameters *mk*, *Pk*, *αk*,*i*, *βk*,*<sup>i</sup>* can be calculated as:

$$\begin{aligned} m\_k &= m\_k^- + P\_k^- H\_k^T (H\_k P\_k^- H\_k^T + \Sigma\_k)^{-1} (y\_k - H\_k m\_k^-) \\ P\_k &= P\_k^- - P\_k^- H\_k^T (H\_k P\_k^- H\_k^T + \hat{\Sigma}\_k)^{-1} H\_k P\_k^- \\ \alpha\_{k,i} &= 1/2 + \alpha\_{k-1,i} \\ \beta\_{k,i} &= \beta\_{k-1,i} + \frac{1}{2} [(y\_k - H\_k m\_k)\_i^2 + (H\_k P\_k H\_k^T)\_{ii}] \end{aligned} \tag{25}$$

where *i* = 1, ..., *d* and *d* denote the dimensionality of the measurement. *m*− *<sup>k</sup>* and *P*<sup>−</sup> *<sup>k</sup>* are the predicted state estimate and its covariance, respectively. Here we assume each component of the measurement noise variance is mutually independent, and then the covariance matrix is diagonal, estimated as:

$$\hat{\Sigma}\_k = \text{diag}\left(\beta\_{k,1} / \alpha\_{k,1}, \dots, \beta\_{k,d} / \alpha\_{k,d}\right) \tag{26}$$

Furthermore, in order to describe the dynamics of the measurement noise variance, the inverse Gamma distribution parameters are assumed to change by a scale factor *ρ* ∈ (0, 1]. The formulas are given as:

$$\begin{aligned} \alpha\_{k,i}^- &= \rho\_i \alpha\_{k-1,i} \\ \beta\_{k,i}^- &= \rho\_i \beta\_{k-1,i} \end{aligned} \tag{27}$$

Note that the value *ρ* = 1 corresponds to stationary variances and lower values represent larger time-fluctuations.

The above are the basic equations of adaptive Kalman filter based on VB approximation. Moreover, when the state equation and the measurement equation are nonlinear, the VB method can be rewritten in the EKF framework.

#### **4. Variational Bayesian Approximation-Based Adaptive Dual Extended Kalman Filter**

In order to handle the joint estimation of the SOC and the battery model parameters as well as the unknown measurement noise covariances, we propose a variational Bayesian approximation-based dual extended Kalman filter (VB-ADEKF) in this paper.

First, let us rewrite the state-space equations for SOC estimation as

$$\begin{cases} X\_{k+1} = F\_k X\_k + G\_k I\_{L,k} + w\_k \\\ y\_k = h(X\_{k'} I\_{L,k'} \theta\_k) + v\_k \\\ w\_k \sim \mathcal{N}(0, Q\_k^x) \\\ v\_k \sim \mathcal{N}(0, \Sigma\_k^x) \end{cases} \tag{28}$$

and the state-space equations for battery model parameters as

$$\begin{cases} \theta\_{k+1} = \theta\_k + r\_k \\ d\_k = h(\mathcal{X}\_{k'} I\_{k,k'} \theta\_k) + c\_k \\ r\_k \sim \mathcal{N}(0, Q\_k^{\theta}) \\ c\_k \sim \mathcal{N}(0, \Sigma\_k^{\theta}) \end{cases} \tag{29}$$

where the process noise covariances *Q<sup>x</sup> <sup>k</sup>* , *<sup>Q</sup><sup>θ</sup> <sup>k</sup>* are assumed known, and the measurement noise variances Σ*x <sup>k</sup>* , <sup>Σ</sup>*<sup>θ</sup> <sup>k</sup>* are unknown, being assumed as stochastic variables.

Then, VB-ADEKF is to alternatively solve the joint posterior distribution *p*(*Xk*, Σ*<sup>x</sup> <sup>k</sup>* |*y*1:*k*) of the SOC and its measurement noise variance and the joint posterior distribution *p*(*θk*, Σ*<sup>θ</sup> <sup>k</sup>* |*y*1:*k*) of the battery parameter and its measurement noise variance by VB approximation as follows:

$$p(X\_k \mid \underline{\boldsymbol{X}}\_k^{\boldsymbol{x}} | \underline{\boldsymbol{y}}\_{1:k}) \approx Q\_{\boldsymbol{x}}(X\_k) Q\_{\boldsymbol{\Sigma}}(\boldsymbol{\Sigma}\_k^{\boldsymbol{x}}) \tag{30}$$

$$p(\theta\_k, \Sigma\_k^{\theta} | y\_{1:k}) \approx Q\_{\theta}(\theta\_k) Q\_{\Sigma}(\Sigma\_k^{\theta}) \tag{31}$$

where *Qx*(*Xk*), *Q*Σ(Σ*<sup>x</sup> <sup>k</sup>* ), *<sup>Q</sup><sup>θ</sup>* (*θk*), and *<sup>Q</sup>*Σ(Σ*<sup>θ</sup> <sup>k</sup>* ) can be solved by Equations (23) and (24). It should be noted that the parameters in these approximating densities are actually obtained using the extended Kalman filters since the measurement function *h*(·) is nonlinear. The filtering procedure of the VB-ADEKF algorithm is summarized in Algorithm 1.

**Algorithm 1** : VB-ADEKF.

**(1) Initialization**: *X*ˆ 0, ˆ *θ*0, *Px*,0, *Pθ*,0, *Q<sup>x</sup>* <sup>0</sup>, *<sup>Q</sup><sup>θ</sup>* <sup>0</sup>, *<sup>α</sup>*<sup>ˆ</sup> *<sup>x</sup>*,0, *<sup>β</sup>*<sup>ˆ</sup> *<sup>x</sup>*,0, *<sup>α</sup>*<sup>ˆ</sup> *<sup>θ</sup>*,0, *<sup>β</sup>*<sup>ˆ</sup> *θ*,0 **(2) Prediction:** *X*ˆ <sup>−</sup> *<sup>k</sup>* <sup>=</sup> *Fk*−1*X*<sup>ˆ</sup> *<sup>k</sup>*−<sup>1</sup> <sup>+</sup> *Gk*−<sup>1</sup> *IL*,*k*−1, <sup>ˆ</sup> *θ*− *<sup>k</sup>* <sup>=</sup> <sup>ˆ</sup> *θk*−<sup>1</sup> *P*− *<sup>x</sup>*,*<sup>k</sup>* <sup>=</sup> *Fk*−<sup>1</sup>*Px*,*k*−1*F<sup>T</sup> <sup>k</sup>*−<sup>1</sup> <sup>+</sup> *<sup>Q</sup><sup>x</sup> <sup>k</sup>* , *P*<sup>−</sup> *<sup>θ</sup>*,*<sup>k</sup>* <sup>=</sup> *<sup>P</sup>θ*,*k*−<sup>1</sup> <sup>+</sup> *<sup>Q</sup><sup>θ</sup> k α*ˆ − *<sup>x</sup>*,*<sup>k</sup>* <sup>=</sup> *<sup>ρ</sup>xα*<sup>ˆ</sup> *<sup>x</sup>*,*k*−1, *<sup>β</sup>*ˆ<sup>−</sup> *<sup>x</sup>*,*<sup>k</sup>* <sup>=</sup> *<sup>ρ</sup>xβ*<sup>ˆ</sup> *x*,*k*−1 *α*ˆ − *<sup>θ</sup>*,*<sup>k</sup>* <sup>=</sup> *ρθα*<sup>ˆ</sup> *<sup>θ</sup>*,*k*−1, *<sup>β</sup>*ˆ<sup>−</sup> *<sup>θ</sup>*,*<sup>k</sup>* <sup>=</sup> *ρθβ*<sup>ˆ</sup> *θ*,*k*−1

where *αx*, *β<sup>x</sup>* and *αθ*, *βθ* are the inverse gamma distribution parameters of the measurement noise covariance, *ρ<sup>x</sup>* and *ρθ* are the scale factors.

**(3) Update:** the update of VB-ADEKF utilizes iterate filtering framework.

**First** set: *<sup>X</sup>*<sup>ˆ</sup> (0) *<sup>k</sup>* <sup>=</sup> *<sup>X</sup>*<sup>ˆ</sup> <sup>−</sup> *<sup>k</sup>* , *<sup>P</sup>*(0) *<sup>x</sup>*,*<sup>k</sup>* = *P*<sup>−</sup> *<sup>x</sup>*,*k*, <sup>ˆ</sup> *θ* (0) *<sup>k</sup>* <sup>=</sup> <sup>ˆ</sup> *θ*− *<sup>k</sup>* , *<sup>P</sup>*(0) *<sup>θ</sup>*,*<sup>k</sup>* = *P*<sup>−</sup> *θ*,*k α*ˆ *<sup>x</sup>*,*<sup>k</sup>* = 1/2 + *α*ˆ <sup>−</sup> *<sup>x</sup>*,*k*, *<sup>β</sup>*ˆ(0) *<sup>x</sup>*,*<sup>k</sup>* <sup>=</sup> *<sup>β</sup>*ˆ<sup>−</sup> *<sup>x</sup>*,*<sup>k</sup> α*ˆ *<sup>θ</sup>*,*<sup>k</sup>* = 1/2 + *α*ˆ <sup>−</sup> *<sup>θ</sup>*,*k*, *<sup>β</sup>*ˆ(0) *<sup>θ</sup>*,*<sup>k</sup>* <sup>=</sup> *<sup>β</sup>*ˆ<sup>−</sup> *θ*,*k* **For** *n* = 1 : *N*, iterate the following *N* (*N* denotes iterated times) steps:

• **Measurement variances**:

Σˆ (*n*) *<sup>x</sup>*,*<sup>k</sup>* <sup>=</sup> *<sup>β</sup>*ˆ(*n*) *<sup>x</sup>*,*<sup>k</sup>* /*α*ˆ (*n*) *<sup>x</sup>*,*<sup>k</sup>* , <sup>Σ</sup><sup>ˆ</sup> (*n*) *<sup>θ</sup>*,*<sup>k</sup>* <sup>=</sup> *<sup>β</sup>*ˆ(*n*) *<sup>θ</sup>*,*<sup>k</sup>* /*α*ˆ (*n*) *θ*,*k* • **State estimate and its covariance**: *X*ˆ (*n*+1) *<sup>k</sup>* <sup>=</sup> *<sup>X</sup>*<sup>ˆ</sup> <sup>−</sup> *<sup>k</sup>* + *P*<sup>−</sup> *x*,*kH*<sup>T</sup> *x*,*k*(*Hx*,*kP*<sup>−</sup> *x*,*kH*<sup>T</sup> *<sup>x</sup>*,*<sup>k</sup>* <sup>+</sup> <sup>Σ</sup><sup>ˆ</sup> (*n*) *<sup>x</sup>*,*<sup>k</sup>* )−1(*yk* <sup>−</sup> *<sup>h</sup>*(*X*<sup>ˆ</sup> <sup>−</sup> *<sup>k</sup>* , *IL*,*k*, <sup>ˆ</sup> *θ*− *<sup>k</sup>* )) *P*(*n*+1) *<sup>x</sup>*,*<sup>k</sup>* = *P*<sup>−</sup> *<sup>x</sup>*,*<sup>k</sup>* − *P*<sup>−</sup> *x*,*kH*<sup>T</sup> *x*,*k*(*Hx*,*kP*<sup>−</sup> *x*,*kH*<sup>T</sup> *<sup>x</sup>*,*<sup>k</sup>* <sup>+</sup> <sup>Σ</sup><sup>ˆ</sup> (*n*) *<sup>x</sup>*,*<sup>k</sup>* )−<sup>1</sup>*Hx*,*kP*<sup>−</sup> *x*,*k* • **Battery parameters estimate and its covariance**: ˆ *θ* (*n*+1) *<sup>k</sup>* <sup>=</sup> <sup>ˆ</sup> *θ*− *<sup>k</sup>* + *P*<sup>−</sup> *θ*,*kH*<sup>T</sup> *θ*,*k*(*Hθ*,*kP*<sup>−</sup> *θ*,*kH*<sup>T</sup> *<sup>θ</sup>*,*<sup>k</sup>* <sup>+</sup> <sup>Σ</sup><sup>ˆ</sup> (*n*) *<sup>θ</sup>*,*<sup>k</sup>* )−1(*yk* <sup>−</sup> *<sup>h</sup>*(*X*<sup>ˆ</sup> <sup>−</sup> *<sup>k</sup>* , *IL*,*k*, <sup>ˆ</sup> *θ*− *<sup>k</sup>* )) *P*(*n*+1) *<sup>θ</sup>*,*<sup>k</sup>* = *P*<sup>−</sup> *<sup>θ</sup>*,*<sup>k</sup>* − *P*<sup>−</sup> *θ*,*kH*<sup>T</sup> *θ*,*k*(*Hθ*,*kP*<sup>−</sup> *θ*,*kH*<sup>T</sup> *<sup>θ</sup>*,*<sup>k</sup>* <sup>+</sup> <sup>Σ</sup><sup>ˆ</sup> (*n*) *<sup>θ</sup>*,*<sup>k</sup>* )−1*Hθ*,*kP*<sup>−</sup> *θ*,*k* • **Parameters for the measurement noise variances estimation**: *β*ˆ(*n*+1) *<sup>x</sup>*,*<sup>k</sup>* <sup>=</sup> *<sup>β</sup>*ˆ<sup>−</sup> *<sup>x</sup>*,*<sup>k</sup>* <sup>+</sup> <sup>1</sup> 2 *yk* <sup>−</sup> *<sup>h</sup>*(*X*<sup>ˆ</sup> (*n*+1) *<sup>k</sup>* , *IL*,*k*, <sup>ˆ</sup> *θ* (*n*+1) *<sup>k</sup>* ) 2 + <sup>1</sup> <sup>2</sup> *Hx*,*kP*(*n*+1) *<sup>x</sup>*,*<sup>k</sup> <sup>H</sup>*<sup>T</sup> *x*,*k β*ˆ(*n*+1) *<sup>θ</sup>*,*<sup>k</sup>* <sup>=</sup> *<sup>β</sup>*ˆ<sup>−</sup> *<sup>θ</sup>*,*<sup>k</sup>* <sup>+</sup> <sup>1</sup> 2 *yk* <sup>−</sup> *<sup>h</sup>*(*X*<sup>ˆ</sup> (*n*+1) *<sup>k</sup>* , *IL*,*k*, <sup>ˆ</sup> *θ* (*n*+1) *<sup>k</sup>* ) 2 + <sup>1</sup> <sup>2</sup> *<sup>H</sup>θ*,*kP*(*n*+1) *<sup>θ</sup>*,*<sup>k</sup> <sup>H</sup>*<sup>T</sup> *θ*,*k* **End for**. **And** set *β*ˆ *x*,*k*=*β*ˆ(*N*) *<sup>x</sup>*,*<sup>k</sup>* , *<sup>β</sup>*<sup>ˆ</sup> *θ*,*k*=*β*ˆ(*N*) *<sup>θ</sup>*,*<sup>k</sup>* , *<sup>X</sup>*<sup>ˆ</sup> *<sup>k</sup>* <sup>=</sup> *<sup>X</sup>*<sup>ˆ</sup> (*N*) *<sup>k</sup>* , *Px*,*<sup>k</sup>* <sup>=</sup> *<sup>P</sup>*(*N*) *<sup>x</sup>*,*<sup>k</sup>* , <sup>ˆ</sup> *θ<sup>k</sup>* = ˆ *θ* (*N*) *<sup>k</sup>* , *<sup>P</sup>θ*,*<sup>k</sup>* <sup>=</sup> *<sup>P</sup>*(*N*) *θ*,*k*

By alternatively using two VB-based extended Kalman filters for online estimation of the battery SOC and model parameters, while compensating for the uncertainties in the model parameters by simultaneous estimation of the measurement noise variances, the adaptability of the proposed algorithm to dynamic changes in battery characteristics is greatly improved. Hence, it is very promising to further increase the SOC estimation accuracy and robustness.

#### **5. Experimental Verification and Analysis**

#### *5.1. Experimental Settings*

The proposed method is experimentally evaluated in this section. The experimental setup for the tests is shown in Figure 2. The setup consists of three lithium-ion battery cells connected in series, an electronic load, and a host computer. The tested lithium-ion battery cells are type 18,650, whose

nominal capacity is 2200 mAh, nominal voltage is 3.7 V, charging and discharging cutoff voltages are 4.2 V and 3 V, respectively. The type of the electronic load is IT8516S produced by ITECH, whose current measurement accuracy is ±(0.1% + 0.1% full scale), and voltage measurement accuracy is ±(0.02% + 0.02% full scale). The load current, terminal voltage, and SOC can be recorded via the host computer during the discharge test.

First, a sequence of pulsed discharging experiments were implemented to determine the relationship of open circuit voltage (OCV) and SOC. By using the curve-fitting toolbox in MATLAB, *k*0, *k*1, *k*2, *k*3, *k*4, *k*<sup>5</sup> were identified, which are shown in Table 1. The measured data and fitted curve are presented in Figure 3. The R-square is used to represent the goodness of fit. The normal value range of the R-square is 0–1 and a value closer to 1 indicates a better fitting curve [32]. It can be seen that the curve fits well with the measurement data, indicating that the selected fifth-order polynomial model can describe the relationship of OCV and SOC very well.

**Figure 2.** The experimental setup.

**Table 1.** The identification results of *k*0, *k*1, *k*2, *k*3, *k*4, *k*5.


Then, a constant current discharge test, a pulse current discharge test, and an urban dynamometer driving schedule (UDDS) test were performed. They are commonly used to verify the performances of SOC estimation methods in EVs. In the constant current discharge test, the current keeps invariant but the terminal voltage declines continuously. In the pulse current discharge test, the current stays at 1 A for 10 min and then decreases to 0 and lasts for 30 min. The process is repeated until the battery reaches the lower cut-off voltage. UDDS, also known as FTP72, was used as a test procedure to certify vehicle emissions by the US Environmental Protection Agency. It can simulate the actual driving conditions of vehicles on the road. The battery current and voltage are both sampled at 1 s. In each test, the true SOC was obtained using CC method. The estimation accuracy, convergence rate, and robustness of the proposed VB-ADEKF are evaluated by comparison with DEKF under different tests.

**Figure 3.** The relationship curve of open circuit voltage (OCV) versus state of charge (SOC).

#### *5.2. Constant Current Discharge Test*

The experiment was performed with a constant discharge current of 1 A. The initial SOC value is set to be 0.8, rather than the real SOC of 1.0. The process noise covariances are set as *Q<sup>x</sup> <sup>k</sup>* = <sup>1</sup> × <sup>10</sup>−<sup>6</sup> *<sup>I</sup>*<sup>2</sup> and *Q<sup>θ</sup> <sup>k</sup>* = <sup>1</sup> ×10−<sup>6</sup> *<sup>I</sup>*3. The measurement noise variances used for DEKF are <sup>Σ</sup>*<sup>x</sup> <sup>k</sup>* <sup>=</sup> 0.001 and <sup>Σ</sup>*<sup>θ</sup> <sup>k</sup>* = 0.0005. The scale factors for VB-ADEKF are set to 1 × <sup>10</sup>−4. The initial parameters *<sup>α</sup>*<sup>0</sup> and *<sup>β</sup>*<sup>0</sup> for battery parameters and SOC are both set as 10 and 0.001, respectively.

The estimated values of the battery model parameters by VB-ADEKF are presented in Figure 4. The ohmic resistance is stable at the beginning of the discharge and increases at the end of the discharge. The polarization resistance decreases at first and increases later with the depth of the discharge and again decreases in the last 1000 s. The polarization capacitance shows a declining trend overall in the process of discharge.

**Figure 4.** Results of parameter identification using variational Bayesian approximation-based adaptive dual extended Kalman filter (VB-ADEKF) in the constant current discharge test.

Furthermore, to verify the accuracy of online identification of the battery parameters by the VB-ADEKF algorithm, the measured terminal voltage and the estimated terminal voltage by VB-ADEKF were compared, as shown in Figure 5. The maximum absolute error was 0.023 V, except for the first big error caused by an incorrect initial SOC value. The mean absolute error was 0.0050 V and the relative mean absolute error was 0.052%. It is clear that the estimated terminal voltage agrees well with the measured voltage. This illustrates the effectiveness of the battery model whose parameters are identified by the proposed VB-ADEKF method.

**Figure 5.** Experimental terminal voltage results using VB-ADEKF in the constant current discharge test.

The SOC estimation results using VB-ADEKF and DEKF are shown in Figure 6. Clearly, VB-ADEKF has a much more accurate SOC estimation than DEKF. The SOC estimation error of VB-ADEKF is bounded within ±1% for most of the time, but DEKF goes outside of this interval. The detailed error values are shown in Table 2. The maximum absolute estimation error of VB-ADEKF is 1.28% and the mean absolute error is 0.64% after discharge for 12 min. Both are significantly smaller than the errors of DEKF, which are correspondingly 2.76% and 1.39%. Meanwhile, from the figure it can be seen that VB-ADEKF converges much faster than DEKF under initial SOC errors. This is also verified in Table 2. The convergence time, which is defined as the first time instant at which the SOC estimation error decreases to ±5%, are 10 s and 335 s for VB-ADEKF and DEKF, respectively.


Maximum Absolute Error 2.76% 1.28% 4.93% 5.00% 4.72% 4.10% Mean Absolute Error 1.39% 0.64% 1.01% 0.68% 1.26% 0.89% Convergence Time 335 s 10 s 698 s 690 s 675 s 603 s

**Table 2.** Statistical analysis of state of charge (SOC) estimation error (after 12 min) in three tests.

#### *5.3. Pulse Current Discharge Test*

In the pulse current discharge test, the initial SOC value is set to be 0.8 while the real SOC is 1.0. The process noise covariances are set as *Q<sup>x</sup> <sup>k</sup>* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>6</sup> *<sup>I</sup>*<sup>2</sup> and *<sup>Q</sup><sup>θ</sup> <sup>k</sup>* = <sup>1</sup> × <sup>10</sup>−<sup>6</sup> *<sup>I</sup>*3. The measurement noise variances used for DEKF are Σ*<sup>x</sup> <sup>k</sup>* <sup>=</sup> 0.01 and <sup>Σ</sup>*<sup>θ</sup> <sup>k</sup>* = 0.0005. The initial parameter values of VB-ADEKF are set as *αθ*,0 = 100, *βθ*,0 = 0.0005, *αx*,0 = 10, *βx*,0 = 0.001.

**Figure 6.** SOC estimation results using VB-ADEKF and dual extended Kalman filter (DEKF) in the constant current discharge test.

Figure 7 presents the estimated values of the battery model parameters. Clearly, there is a stepped increase in the ohmic resistance from 34 mΩ to 49 mΩ. The polarization resistance retains stable during the entire discharge process. The polarization capacitance first decreases rapidly then stays stable until the end of the discharge. Figure 8 shows the the measured terminal voltage and the estimated terminal voltage by VB-ADEKF to verify the effectiveness of the battery parameters identification. The maximum and mean absolute estimation errors are 0.086 V and 0.0027 V, respectively. It implies that the estimated terminal voltage has a good agreement with the measured voltage. It further shows the battery model parameters are well identified.

**Figure 7.** Results of parameters identification using VB-ADEKF in the pulse current discharge test.

The SOC estimation results of VB-ADEKF and DEKF are plotted in Figure 9. It shows that VB-ADEKF and DEKF have comparable performance in convergence rate. However VB-ADEKF has more accurate SOC estimation than that of DEKF. This is also verified in Table 2, in which the mean absolute error of VB-ADEKF is 0.68%, while it is 1.01% for DEKF after discharge for 12 min.

**Figure 8.** Experimental terminal voltage results using VB-ADEKF in the pulse current discharge test.

**Figure 9.** SOC estimation results using VB-ADEKF and DEKF in the pulse current discharge test.

#### *5.4. UDDS Test*

To evaluate the SOC estimation performance under dynamic loading profiles, an UDDS cycle was performed on the battery cells. According to the actual tolerable currents of the lithium-ion battery cells, the loading currents are scaled down, as shown in Figure 10. The initial SOC value is set to 0.8. The process noise covariances are set as *Q<sup>x</sup> <sup>k</sup>* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>6</sup> *<sup>I</sup>*<sup>2</sup> and *<sup>Q</sup><sup>θ</sup> <sup>k</sup>* = <sup>1</sup> × <sup>10</sup>−<sup>6</sup> *<sup>I</sup>*3. The measurement noise variances used for DEKF are Σ*<sup>x</sup> <sup>k</sup>* <sup>=</sup> 0.01 and <sup>Σ</sup>*<sup>θ</sup> <sup>k</sup>* = 0.0005. The initial parameter values of VB-ADEKF are set as *αθ*,0 = 100, *βθ*,0 = 0.0005, *αx*,0 = 10, *βx*,0 = 0.001.

The estimated values of the battery model parameters are presented in Figure 11. It can be seen that the range of parameter values of *R*0, *R*<sup>1</sup> and *C*<sup>1</sup> are consistent with these in constant current discharge test and pulse current discharge test, but the changing rules are slightly different. Figure 12 presents the the measured terminal voltage and the estimated terminal voltage by VB-ADEKF to verify the effectiveness of the battery parameters identification. The maximum and mean absolute estimation errors are 0.062 V and 0.0011 V, respectively. It is clear that the estimated terminal voltage has a good consistency with the measured voltage. It further shows the battery model parameters are well identified.

**Figure 10.** Current profiles in the urban dynamometer driving schedule (UDDS) test.

**Figure 11.** Results of parameters identification using VB-ADEKF in the UDDS test.

**Figure 12.** Experimental terminal voltage results using VB-ADEKF in the UDDS test.

Figure 13 presents the SOC estimation results of VB-ADEKF and DEKF. It shows that both VB-ADEKF and DEKF have good estimation accuracy when SOC is between 20% and 90%. But VB-ADEKF still outperforms the traditional DEKF in SOC estimation accuracy and convergence rate. A comparative summary of the two methods is given in Table 2. When SOC is decreased to 20%, the estimation error of both methods begins to increase, but is no larger than 5%. This may be because the polarization effect of the battery is further aggravated at lower SOC level. As a result, the terminal voltage measurement error has increased, thereby reducing the SOC estimation accuracy.

**Figure 13.** SOC estimation results using VB-ADEKF and DEKF in the UDDS test.

#### *5.5. Convergence Ability with Initial SOC Error*

Because it is difficult to determine the initial SOC value precisely in practical applications, it is important and indispensable for the algorithms to have the ability to correct the uncertainty brought by the initial SOC error. Therefore, the convergence rate with initial SOC error is adopted as another indicator for evaluating the SOC estimation algorithms. The true initial SOC value is 1.0. Figures 14–16 present the estimation results of VB-ADEKF and DEKF algorithms for different initial SOC values from 30% to 90% under the above three tests, respectively. Overall, the convergence time increases as the initial SOC error rises for both VB-ADEKF and DEKF. But the growth rates are very different for the two filters under different tests.

**Figure 14.** SOC estimation results of the two filters with different initial SOC values in the constant current discharge test.

**Figure 15.** SOC estimation results of the two filters with different initial SOC values in the pulse current discharge test.

Specifically, from the results of the convergence time shown in Table 3, it can be seen that the convergence time of VB-ADEKF is stabilized about 10 s into the constant current discharge test when the initial SOC value is larger than 50%, and increased to more than 100 s when the initial SOC values are 40% and 30%. But the convergence time of DEKF, which is between 300 s and 500 s, is much longer than VB-ADEKF. It shows that VB-ADEKF can quickly converge to the true SOC values without resulting in the accumulation of errors caused by the initial error of the SOC. In the pulse current discharge test, the convergence rates are comparable for VB-ADEKF and DEKF in the case of small initial SOC errors, for example, 10% or 20%. However, when the initial SOC error is relatively large, the convergence time of VB-ADEKF becomes much smaller than that of DEKF. In the UDDS test, VB-ADEKF only exhibits a slightly increasing trend with the increase of the initial SOC error. However, DEKF converges slower than VB-ADEKF, in the meantime, with a larger SOC estimation error. Its convergence time goes up quickly as the initial SOC error increases. This implies that the initial SOC error has a noticeable impact on the performance of DEKF. But, from the overall perspective, VB-ADEKF is not very sensitive to the initial SOC error. It shows that the proposed VB-ADEKF has better robustness for initial SOC errors than the traditional DEKF.

**Figure 16.** SOC estimation results of the two filters with different initial SOC values in the UDDS test.


**Table 3.** Convergence time (s) of variational Bayesian approximation-based adaptive dual extended Kalman filter (VB-ADEKF) and dual extended Kalman filter (DEKF) with different initial SOC values.

#### *5.6. Effect of Mistuning*

If the working condition of the battery changes abruptly, the SOC measurement error would probably be varied largely with before, so the prior tuning of the measurement variance of the DEKF will not give an optimal estimate of the SOC. But there is no such issue in the proposed VB-ADEKF since the measurement variances are estimated online. This effects of mistuning brought about by inappropriate measurement variance of the traditional DEKF and VB-ADEKF are mainly reflected in the SOC estimation error, as shown in Figures 17–19. Here, the measurement noise variances of DEKF are mistuned to Σ*<sup>x</sup> <sup>k</sup>* <sup>=</sup> 0.01 and <sup>Σ</sup>*<sup>θ</sup> <sup>k</sup>* = 0.001 in constant current discharge test, and <sup>Σ</sup>*<sup>x</sup> <sup>k</sup>* = 0.1 and Σ*θ <sup>k</sup>* = 0.005 in the pulse current discharge test and UDDS test. The initial estimates of the measurement variances of VB-ADEKF are also correspondingly mistuned. From the results, we can see that the convergence rate of DEKF slows down and the SOC estimation accuracy also declines in the case of mistuning, while the SOC estimation performance of the proposed VB-ADEKF remains almost unchanged. This suggests that the proposed VB-ADEKF is more robust than the traditional DEKF.

**Figure 17.** SOC estimation results using VB-ADEKF and DEKF in the case of mistuning in the constant current discharge test.

**Figure 18.** SOC estimation results using VB-ADEKF and DEKF in the case of mistuning in the pulse current discharge test.

**Figure 19.** SOC estimation results using VB-ADEKF and DEKF in the case of mistuning in the UDDS test.

#### **6. Conclusions**

To deal with the measurement statistical uncertainties and inaccurate battery model, a variational Bayesian approximation-based adaptive dual extended Kalman filter (VB-ADEKF) is proposed in this paper for SOC estimation of lithium-ion batteries. First, the variational Bayesian inference is integrated with the extended Kalman filter to jointly estimate the states and the measurement noise covariances. Then, two VB-based extended Kalman filters are alternatively used for online estimation of the battery SOC and model parameters, while simultaneously estimating the measurement noise variances to compensate for the uncertainties in the measurement and battery parameters. Therefore, the adaptability of the proposed algorithm to dynamic changes in battery characteristics is greatly improved. The effectiveness and superiority of the proposed algorithm have been verified by comparing with the dual EKF (DEKF) algorithm through experiments under the constant current discharge test, pulse current discharge test, and UDDS test. The results show that the proposed VB-ADEKF algorithm outperforms the traditional DEKF approach in terms of SOC estimation accuracy and convergence rate. Especially, when the quality of measurements changes with the operating conditions, the proposed VB-ADEKF exhibits better robustness than DEKF.

**Author Contributions:** J.H. conceived this paper, designed the experiments, and analyzed the data; T.G. and H.H. performed the experiments; Y.Y. revised the paper and provided some valuable suggestions.

**Funding:** This research was funded by Shaanxi Provincial Key Research and Development Programs (No. 2017ZDXM-GY-06, No. 2017GY-057), and Xi'an Science and Technology Planning Project–Scientific and Technological Innovation Guidance Project ( No. 201805042YD20CG26 (8)).

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


c 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/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18