A novel composite material with a NF@NCO/NCO NFAs-like structure was fabricated on an NF substrate using a facile hydrothermal process. The excellent supercapacitance performance of the NF@NCO/NCO NFAs was found to depend mainly on easy pathways between the transportation of electrons and ions, a large electrochemical surface area, more accessible electroactive sites, and great adhesion to the NF substrate to enhance energy storage and exhibit higher electrochemical performance than the NF@NCO NLAs electrode. The nanofile-like NF@NCO/NCO NFAs covered the NF@NCO NLAs surface area fully. As a result, the synthesized hierarchical nanofile-like NF@NCO/NCO NFAs composite electrode exhibited superior energy storage performance, including a high specific capacitance of 2312 F g−1 at a current density of 2 mA cm−2, as well as a higher cycling stability at a current density of 5 mA cm−2 (98.7% retention after 5000 cycles) than the NF@NCO NLAs (96.3% retention after 5000 cycles). These results demonstrate the outstanding cycling stability, excellent electrochemical performance, and higher electrical conductivity of the hierarchical nanofile-like NF@NCO/NCO NFAs composite compared to the NF@NCO NLAs composite, highlighting its potential for flexible energy storage systems.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1996-1073/12/7/1308/ s1, Figure S1: Schematic illustration of NF@NCO/NCO NFAs.

**Author Contributions:** Conceptualization, A.K.Y.; Formal analysis, H.-J.K. and E.R.A.; Investigation, A.K.Y. and H.-J.K.; Methodology, A.K.Y.; Supervision, A.K.Y.; Writing—original draft, A.K.Y.; Writing—review & editing, A.K.Y.

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

**Acknowledgments:** This work was financially supported by BK 21 PLUS, the Creative Human Resource Development Program for IT Convergence (NRF-2015R1A4A1041584), Pusan National University, Busan, South Korea. We would like to thank KBSI, Busan for the SEM, TEM, XRD, XPS, and EDX analyses.

**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* **Hybrid Model Predictive Control Strategy of Supercapacitor Energy Storage System Based on Double Active Bridge**

### **Lujun Wang 1,\*, Jiong Guo 1, Chen Xu 1, Tiezhou Wu 1 and Huipin Lin 2**

Received: 19 April 2019; Accepted: 30 May 2019; Published: 4 June 2019

**Abstract:** In order to solve the problem of which the dynamic response of a supercapacitor (SC) is limited due to the mismatch dynamic characteristics between the DC/DC converter and supercapacitor in an energy storage system, this paper proposes a hybrid model predictive control strategy based on a dual active bridge (DAB). The hybrid model predictive control model considers the supercapacitor and DAB in a unified way, including the equivalent series resistance and capacitance parameters of the SC. The method can obtain a large charging and discharging current of the SC, thereby not only improving the overall response speed of the system, but also expanding the actual capacity utilization range of the SC. The simulation results show that compared with the model prediction method of the dual active bridge converter, the proposed control method can effectively improve the overall response speed of the system, which can be improved by at least 0.4 ms. In addition, the proposed method increases the actual upper limit of the SC voltage, reduces the actual lower limit of the SC voltage, and then expands the actual capacity utilization range of the SC by 18.63%. The proposed method has good application prospects in improving the dynamic response performance of energy storage systems.

**Keywords:** dual active bridge; SC; energy storage; fast response; capacity utilization

### **1. Introduction**

Supercapacitors (SCs) have developed rapidly in recent years due to their large capacity, high power density [1], high charge-discharge rate, high conversion efficiency, wide operating temperature range, no pollution, and convenient control. However, due to the wide range of terminal voltages during the storage and release of energy in an SC, the systems used in actual energy storage systems are usually equipped with bidirectional DC/DC converters to implement power conversion [2].

In addition, since an SC voltage unit is low [3], many units should be connected in series in actual use, causing the problem of equalization. For some applications, where the DC bus voltage is relatively large, such as a three-phase grid-connected inverter, a high-gain bidirectional DC converter suitable for a wide voltage range is needed to adjust the output voltage to meet the energy requirements of the DC bus [4]. In recent years, dual active bridge bidirectional DC-DC converters are commonly used in medium power applications [5], which enable a soft switching operation over a wide voltage range. Additionally, the transformer turns ratio can be used to obtain a high gain output [6]. In addition, the DAB converter has the advantages of a high power-density, bidirectional transmission power, high efficiency, and structural symmetry [7,8]. Due to the large number of switches in a DAB converter, the control process is generally complicated, so the response time of the converter is longer than with an

SC [9], which does not reflect the fast response characteristics of SCs. In order to improve the dynamic response speed of the converter when the load is abrupt or the input voltage fluctuates, some scholars have proposed related control methods.

In [10], an inductor current boundary control method was proposed based on tradition inductor current peak control and single-phase shift control. The boundary control method significantly speeds up the dynamic response of the DC-DC converter when the load suddenly changes, but the method is too complicated and requires up to five Hall sensors. In [11], a load current feedforward control method was proposed based on direct power control, which introduces the load current into the control system, significantly improving the load response speed of the converter. However, it did not analyze the response characteristics of the converter when the input voltage was changed. In [12], based on the idea of direct power control, a virtual direct power control strategy was proposed, which introduces input voltage and load current information at the same time, so that the converter achieves a fast dynamic response. In [13], based on the working principle of a dual-active bridge DC-DC converter, the principle of single-phase shift control, and the steady-state characteristics of the output current, an output current control method based on the input voltage feedforward was proposed to improve the dynamic characteristics of the converter. This method is based on the research of the dual active bridge DC-DC converter as the current source. In [14], a model predictive control strategy was proposed based on the establishment of the state space average model of a DAB by analyzing the characteristics of a DAB and using the traditional phase shift control as an example. The proposed method significantly improved the dynamic characteristics of the converter for the input voltage's abrupt change and the load's abrupt change compared with the traditional voltage closed-loop control method.

In order to improve the overall response speed of a super-capacitor energy storage system based on dual active bridges, a hybrid model predictive control strategy with a fast dynamic response time is proposed. First, the overall structure model of the energy storage system based on a combined appropriate SC equivalent model and a DAB converter model is established in this paper. Then, the working principle of a DAB based on traditional phase shift control is analyzed, and the state space average model of the energy storage system composed of a DAB converter and SC is established. By selecting the appropriate evaluation target and establishing the cost function, a model predictive control strategy incorporating the parameters of the SC equivalent model is proposed to manage the energy of the system. Finally, the effectiveness of the proposed control strategy is verified by simulations in Matlab (<IP_ADDRESS>8062 (R2017a), MathWorks, Natick, Massachusetts, America) software and experiments.

### **2. Working Principle of the DAB-SC System**

### *2.1. SC Model*

The structure of the SC unit is mainly composed of electrodes, electrolytes, diaphragms, and packaging, and each component has an important impact on its performance. Thus, it is necessary to build a precise and applicable SC model in circuit analysis. Scholars have proposed many SC equivalent models based on the physical structure and charging and discharging conditions, such as the first-order classical resistor-capacitor (RC) model, RC transmission line model, multi-branch RC parallel model, frequency domain model [15], and artificial neural network (ANN) model. By combining various SC models to obtain a general model, the general SC equivalent circuit model can be simplified based on characteristics, such as: Large-scale modularization of SCs, high charging and discharging power, and high charging and discharging frequency. The simplified model is shown in Figure 1.

**Figure 1.** Simplified capacitor circuit model of the SC unit.

The use of SCs for large-scale energy storage requires identification of the SC parameters. In the simplified equivalent circuit model of the SC unit, RES is the equivalent series resistance, which not only characterizes the heat loss inside an SC, but also causes different voltage drops as the discharge current changes. Thus, RES affects the maximum value of the discharge current to a certain extent. The magnitude of RES generally does not change with external factors [16], and is usually measured by a classical model.

When the charge and discharge conditions of the SC are determined, the main factor that affects the value of the variable capacitor in the equivalent circuit is the actual operating voltage [17]. In order to identify the parameters of the SC unit, several sets of tests are needed, by charging the SC unit to different voltage levels with a constant frequency.

According to Equation (1), the actual capacitance of an SC unit in different sets of tests can be calculated:

$$\mathbf{C\_{SC}} = \frac{\mathbf{Q\_{SC}}}{\Delta l l\_{\rm SC}} = \frac{\int l\_{\rm SC} \, \mathrm{d}t}{\Delta l l\_{\rm SC}} = \frac{\mathbf{I\_{SC}} \Delta T}{\Delta l l\_{\rm C} - \Delta l l\_{\rm ESR}}.\tag{1}$$

In Equation (1), *Q*SC is the amount of charge of an SC; Δ*U*SC is the actual voltage change; ISC is the SC current; ΔT is the charging and discharging time; Δ*U*C is the total voltage change; Δ*U*ESR is the voltage drop caused by ESR.

After processing the (*U*charge, *Cac*tual) data-points by the linear fitting method, the relationship between the capacitance and voltage of a SC can be obtained as shown in Equation (2):

$$\mathbf{C\_{actual}} = \Delta \mathbf{C} + \mathbf{C}\_0 = \mathbf{K\_V} \mathbf{U\_{charge}} + \mathbf{C\_0} \tag{2}$$

In Equation (2), *C*actual is the actual capacitance of an SC at different voltages; Δ*C* is the variable capacitance; C0 is the initial capacitance of an SC; Kv is the coefficient of variation of the SC capacitance with the voltage; *U*charge is the SC's actual voltage.

In this paper, through the linear fitting of graphic data in [18], the voltage and capacitance data of Maxwell's 1 F SC were obtained, and then the data of the 125 V, 63 F SC, which was used in this paper, was obtained according to the string and parallel relationship of the capacitor. The 125 V, 63 F SC is made up of 46 modules in series, and the module was obtained by paralleling 2900 Maxwell's 1 F, 2.7 V SC units. Additionally, the obtained data of the 125 V, 63 F SC were linearly fitted to obtain a relation as shown in Equation (3):

$$C\_{\text{actual}} = 0.1127 U\_{\text{charge}} + 81.384.\tag{3}$$

### *2.2. Principle Analysis of DAB*

The DAB converter acts as an interface between the SC and DC bus. Generally, there are two basic control methods for the DAB converter: The first is pulse width modulation control; the second is phase shift control. Phase shift control is simpler and easier to implement than the first. Therefore, most DAB controls adopt the phase shift control method [19]. Taking single-phase shift as an example, when a single DAB is in normal operation, the switching frequencies of the H-bridges on both sides of the isolation transformer are the same, and the two H-bridges output square-wave voltages, *u*ab and *u*cd, respectively, with a duty ratio of 50%. When the phase of the *u*ab leads the *u*cd, the energy is forwardly transmitted; that is, from *u*dc to *u*o. Otherwise, the energy is reversely transmitted. Therefore, by controlling the degree of the phase shift angle between the two square wave voltages, the transferred energy can be controlled.

Figure 2 shows the topology of the DAB-SC system. Figure 3 shows the main working waveform of the converter with the phase shift control method. The two switches on each bridge arm are complementarily turned on, and the diagonal switches are turned on or off at the same time. The power is controlled by controlling the phase shift angle between the midpoint voltages, *u*ab and *u*cd, of the primary and secondary side arms. Ts is half of the switching frequency, Ts = 1/(2fs), and *D* represents the phase shift between *u*ab and *u*cd, where *D* = *t*on/Ts, and −1 < *D* < 1. The working principles of the converter when transmitting power in the forward and reverse directions are similar. This paper utilized the forward direction as an example to analyze the operating characteristics of the DAB converter. Before the analysis, the following assumptions were made: (1) The converter is in steady state operation; (2) all switches are regarded as ideal switches with a reverse parallel diode, and the parasitic capacitance is connected in parallel; (3) the size of L is equal to the sum of the leakage inductance of the transformer and the applied inductance; (4) the transformer is ideal, the excitation current is zero; and (5) *u*dc > n*u*o, where n is the turns ratio of the transformer.

**Figure 2.** Topology of the DAB-SC system.

**Figure 3.** The main working waveform of the DAB converter with phase shift control: (**a**) Forward power transfer; (**b**) reverse power transfer.

When the converter is working in the forward direction, the driving signals of S1, S4 and S2, S3 are ahead of S5, S8 and S6, S7 respectively, as shown in Figure 3a. In the figure, *u*'cd is the voltage of *u*cd converted to the primary side; *u*L is the voltage across the inductor, L; and *i*L is the inductor current.

Under steady-state conditions, the relationship between the transmission power of the DAB and the amount of phase shift is as shown in Equation (4):

$$P\_o = \frac{\text{n}\mu\_{\text{dc}}\mu\_o}{2\text{f}\_s\text{L}}D(1-|D|). \tag{4}$$

It can be known from Equation (4) that the transmission power reaches maximum when *D* = 0.5 in the forward transmission, and when *D* = −0.5 in the reverse transmission. Typically, the DAB converter operates in the range of *D*∈[−0.5, 0.5].

The output power of the DAB can also be expressed as *P*o = *u*o*i*o, and the relationship between the output current of the DAB converter and the phase shift amount can be obtained as shown in Equation (5):

$$i\_0 = \frac{\mathbf{n}u\_{\rm dc}}{2\mathbf{f}\_{\rm s}\mathbf{L}}D(1-|D|). \tag{5}$$

In the DAB converter, there are two key parameters, which are the transformer ratio, n, and the buffer inductance, L. When the input voltage of the DAB converter is equal to the value of the output voltage converted to the primary side, the current stress of the device is minimized. By considering the actual fluctuation of the input and output voltages, the transformer ratio can be selected to make the rated input voltage equal to the value of the rated output voltage converted to the primary side. For example, when the rated voltages of the input and output voltages are 700 V and 75 V, respectively, the transformer ratio is determined to be n = 700/75 ≈ 9.

The zero voltage switch (ZVS) range of the device and the suppression of the reactive power are mutually contradictory. A larger ZVS range of a device requires larger inductor values; however, for the suppression of the reactive power, the smaller the inductor value, the better [20], which is a contradiction. The focus of this paper was to study the response speed of the system, and thus the influence of ZVS was not considered. Therefore, a small inductance value as selected for the research.

### *2.3. State Space Average Model of DAB-SC*

In order to establish the state space average model of the DAB, the output capacitor voltage and the inductor current were selected for modeling. It can be seen from Figure 3a that when the voltage, *u*ab, leads the voltage, *u*cd, the energy is transmitted forward, 0 < *D* < 0.5. The overall DAB-SC system has six working states in one switching cycle, and a system of differential equations can be established in each working state. According to Kirchhoff's law, the differential equations, Equations (6)–(11), can be obtained:

$$\begin{cases} \begin{array}{c} \frac{d\dot{\mathbf{u}}\_{\star}}{dt} = \frac{\mathbf{n}u\_{0}}{\mathbf{L}} + \frac{\mathbf{u}\_{\mathrm{dc}}}{\mathbf{L}}\\ \frac{d\mathbf{u}\_{\mathrm{c}}}{dt} = -\frac{\mathbf{n}u\_{\mathrm{c}}}{\mathbf{C}} - \frac{u\_{\mathrm{c}} - \mathbf{u}\_{\mathrm{c}}}{\mathbf{R}\_{\mathrm{ES}}\mathbf{C}} \end{array} \begin{array}{c} t \in [0, \frac{D}{2}\mathbf{T}\_{\mathrm{s}}], \end{array} \tag{6}$$

$$\begin{cases} \frac{\text{d}\underline{\text{u}}\_{\text{}}}{\text{d}\underline{\text{t}}} = \frac{\text{nu}\_{\text{o}}}{\text{L}} + \frac{\text{u}\_{\text{dc}}}{\text{L}}\\ \frac{\text{d}\underline{\text{u}}\_{\text{o}}}{\text{d}\underline{\text{t}}} = -\left(\frac{\text{nu}\_{\text{o}}}{\text{C}} - \frac{\text{u}\_{\text{SC}} - \text{u}\_{\text{o}}}{\text{R}\_{\text{ES}}\underline{\text{C}}}\right) \end{cases} \text{t} \in \left[\frac{D}{2}\mathrm{T}\_{\text{s}}, D\Gamma\_{\text{s}}\right],\tag{7}$$

$$\begin{cases} \frac{d\bar{\mathbf{u}}\_{\rm L}}{d\mathbf{f}} = -\frac{\mathbf{n}u\_{\rm o}}{\mathbf{L}} + \frac{u\_{\rm dc}}{\mathbf{L}}\\ \frac{d\mathbf{u}\_{\rm o}}{d\mathbf{f}} = \frac{\mathbf{n}\bar{\mathbf{u}}\_{\rm s}}{\mathbf{C}} - \frac{u\_{\rm s} - u\_{\rm SC}}{\mathbf{R}\_{\rm ES}\mathbf{C}} \end{cases} t \in [\mathbf{DT}\_{\rm s}, \mathbf{T}\_{\rm s}]\_{\rm \prime} \tag{8}$$

$$t \left\{ \begin{array}{c} \frac{\text{d} \mathbf{i}\_{\text{L}}}{\text{d} \mathbf{t}} = -\frac{\text{nu}\_{0}}{\text{L}} - \frac{\text{u}\_{\text{dc}}}{\text{L}}\\ \frac{\text{d} \mathbf{u}\_{\text{0}}}{\text{d} \mathbf{t}} = \frac{\text{n} \mathbf{i}\_{\text{L}}}{\text{C}} - \frac{\text{u}\_{\text{p}} - \text{u}\_{\text{SC}}}{\text{R}\_{\text{FS}} \text{C}} \end{array} \mid t \in \left[\mathbf{T}\_{\text{s} \prime} \left(\mathbf{1} + \frac{\text{D}}{2}\right) \mathbf{T}\_{\text{s}}\right] \tag{9}$$

$$\begin{cases} \frac{d\dot{\mathbf{l}}\_{\rm l}}{dt} = -\frac{\mathbf{n}u\_{\rm o}}{\mathbf{L}} - \frac{u\_{\rm dc}}{\mathbf{L}}\\ \frac{d\dot{u}\_{\rm o}}{dt} = \frac{\mathbf{n}\dot{\mathbf{l}}\_{\rm s}}{\mathbf{C}} + \frac{u\_{\rm SC} - u\_{\rm o}}{\mathbf{R}\_{\rm KS}\mathbf{C}} \end{cases} t \in [(1 + \frac{D}{2})\mathbf{T}\_{\rm s}, (1 + D)\mathbf{T}\_{\rm s}]. \tag{10}$$

$$t \left\{ \begin{array}{c} \frac{d\dot{\mathbf{u}}\_{\rm t}}{d\mathbf{t}} = \frac{\mathbf{n}u\_{\rm o}}{\mathbf{L}} - \frac{u\_{\rm cl}}{\mathbf{L}}\\ \frac{d\mathbf{u}\_{\rm o}}{d\mathbf{t}} = -\frac{\mathbf{n}\dot{\mathbf{u}}\_{\rm s}}{\mathbf{C}} - \frac{u\_{\rm o} - u\_{\rm cl}}{\mathbf{R}\_{\rm IS}\mathbf{C}} \end{array} \; t \in [(1+D)\mathbf{T}\_{\rm s}\mathcal{D}\mathbf{T}\_{\rm s}].\tag{11}$$

For each differential equation, it can only represent the voltage–current relationship in this operating state. Therefore, it is necessary to establish a differential equation that can describe the characteristics of the converter over the entire switching period. It is worth noting that in the DAB, the inductor current is a pure alternative current (AC) component whose average value is zero during one switching cycle, so the average value of the inductor current is meaningless. Therefore, a simplified reduction model of the DAB module can be obtained by eliminating the inductor current as shown in Equation (12):

$$\frac{\mathrm{d}\langle\mu\_{\mathrm{o}}\rangle}{\mathrm{d}t} = \frac{\mathrm{n}D(1-D)\langle\mu\_{\mathrm{dc}}\rangle}{2\mathrm{L}f\_{\mathrm{s}}\mathrm{C}} - \frac{\langle\mu\_{\mathrm{o}}\rangle - \langle\mu\_{\mathrm{SC}}\rangle}{R\_{\mathrm{ES}}\mathrm{C}}.\tag{12}$$

Similarly, when the voltage, *u*cd, leads *u*ab, the energy is reversely transmitted; that is, −0.5 < *D* < 0, and the simplified reduced-order model of the DAB module is shown in Equation (13):

$$\frac{\mathrm{d}\langle\mu\_{\mathrm{o}}\rangle}{\mathrm{d}t} = \frac{\mathrm{n}D(1-D)\langle\mu\_{\mathrm{dc}}\rangle}{2\mathrm{L}f\_{\mathrm{s}}\mathrm{C}} - \frac{\langle\mu\_{\mathrm{o}}\rangle - \langle\mu\_{\mathrm{SC}}\rangle}{R\_{\mathrm{ES}}\mathrm{C}}.\tag{13}$$

In conclusion, when −0.5 < *D* < 0.5, the simplified reduced-order model of the DAB module is as shown in Equation (14):

$$\frac{\mathrm{d}\langle\mu\_{\mathrm{o}}\rangle}{\mathrm{d}t} = \frac{nD(1-|\mathrm{D}|)\langle\mu\_{\mathrm{dc}}\rangle}{2\mathrm{L}f\_{\mathrm{s}}\mathrm{C}} - \frac{\langle\mu\_{\mathrm{o}}\rangle - \langle\mu\_{\mathrm{SC}}\rangle}{R\_{\mathrm{ES}}\mathrm{C}}.\tag{14}$$

### **3. Management Strategy**

The state of charge (SOC) reflects the remaining capacity of the SC. In a photovoltaic system equipped with energy storage, by accurately estimating the SOC, how long the SC can work under a current load can be determined. There is a definite relationship between the energy, *E*SC, stored by the SC and the terminal voltage, as shown in Equation (15):

$$E\_{\rm SC} = \frac{\mathcal{C}\_{\rm SC} \mathcal{U}\_{\rm o}^2}{2}. \tag{15}$$

Therefore, it is only necessary to detect the SC voltage to accurately determine its SOC. The relationship between the voltage of the SC and the SOC is as shown in Equation (16):

$$\text{SOC} = \frac{\frac{\text{C}\_{\text{SC}}\text{U}\_{\text{o}}^{2}}{2}}{\frac{\text{C}\_{\text{SC}}\text{U}\_{\text{N}}^{2}}{2}} = \left(\frac{\text{U}\_{\text{o}}}{\text{U}\_{\text{N}}}\right)^{2}. \tag{16}$$

In Equation (16), CSC is the capacity of the SC, *U*o is the SC's terminal voltage, and UN is the rated voltage of the SC.

Based on the relationship between the terminal voltage and the SOC, the system can be managed by controlling the terminal voltage. Six voltage limits were set: The DC bus voltage upper limit, Udc\_max; lower limit, Udc\_min; SC upper limit voltage, Usc\_max; lower limit voltage, Usc\_min; upper limit warning voltage, Usc\_h; and lower limit warning voltage, Usc\_l. The working principle is as follows: When the DC bus voltage exceeds Udc\_max and the SC voltage is close to Usc\_min, the SC can be charged with a large current. As the SC voltage increases, the phase shift duty ratio of the DAB converter is changed by predictive control in order to prevent overcharging and to increase the upper limit of the actual SC charging voltage; when the DC bus voltage exceeds Udc\_min, and the SC voltage is close to Usc\_max, the SC performs a large current discharge. As the SC voltage decreases, the phase shift duty ratio of the DAB converter is also changed by predictive control in order to prevent over discharging and to decrease the lower limit of the actual SC charging voltage.

According to the State Grid Corporation "DC Power System Technical Supervision Regulations" and other relevant standards and specifications, when the DC bus power supply is continuously supplied, normal fluctuation of the grid voltage should be no more than 10% of the rated voltage [21]. The rated DC bus voltage of the microgrid in this paper was set at 700 V with a fluctuation of ±5%; that is, Udc\_min = 665 V, and Udc\_max = 735 V. The SC had a rated voltage of 125 V and the upper limit warning voltage and the lower limit warning voltage were 40%–80%.

In summary, the coordinated control strategy of the energy storage systems is as follows: a) In the shutdown state of an energy storage system, the DAB converter is in an idle state and is not charged or discharged; b) in the SC charging area, when the DC bus voltage, *U*dc, satisfied that *U*dc > 735 V, and *SOC*sc ≤ 80%, the SC is charged to absorb excess power in the DC bus; and c) in the SC discharge area, when the DC bus voltage satisfied that *U*dc < 665 V, and *SOC*sc > 40%, the SC discharges, supplementing the power shortage in the DC bus.

### **4. Model Predictive Control**

Model predictive control (MPC) is an active control strategy based on the mathematical model of the converter and is adjusted by predicting the change trend of the control variables. It predicts the optimality of the next moment based on the circuit parameters and the sampling information at the current time to minimize the deviation between the output and the reference. There are many forms of the MPC algorithm, and all MPC algorithms have three basic characteristics: Prediction model, rolling optimization, and feedback correction. A schematic diagram of the MPC strategy of the DAB-SC system studied in this paper is shown in Figure 4, where *u*o is the DAB output voltage, *u*dc is the DC bus voltage, and *i*SC is the SC current.

**Figure 4.** Schematic diagram of the MPC strategy.

### *4.1. Establishment of the Prediction Model for the DAB-SC System*

The derivative of the output voltage in Equation (14) reflects the trend of the output voltage. Additionally, the Euler forward method can be used to discretize Equation (14), and then Equation (17) can be obtained:

$$\frac{\mathbf{d}u\_{\rm o}}{\mathbf{d}t} = \frac{u\_{\rm o}(t\_{\rm k+1}) - u\_{\rm o}(t\_{\rm k})}{t\_{\rm k+1} - t\_{\rm k}} = \frac{u\_{\rm o}(t\_{\rm k+1}) - u\_{\rm o}(t\_{\rm k})}{2\mathbf{T}\_{\rm s}}.\tag{17}$$

Equation (17) can be simplified with Equation (14) to obtain Equation (18):

$$u\_{\rm o}(t\_{\rm k+1}) = \frac{\rm nD(1-|D|)u\_{\rm dc}(t\_{\rm k})}{2\rm LCf\_{\rm s}^{2}} - \frac{u\_{\rm o}(t\_{\rm k}) - u\_{\rm SC}(t\_{\rm k})}{R\_{\rm ES}\rm Cf\_{\rm s}} + u\_{\rm o}(t\_{\rm k}).\tag{18}$$

The SC has a high-power density characteristic and can provide a large instantaneous current, thus the current flowing through SC at time *t*k+1 was set to be *i*SC(tk+1). Then, Equation (19) can be established: *u*o(*t*k)−*u*SC(*t*k)

$$\frac{\mu\_{\rm o}(t\_{\rm k}) - \mu\_{\rm SC}(t\_{\rm k})}{R\_{\rm ES}} = i\_{\rm SC}(t\_{\rm k+1}). \tag{19}$$

By substituting Equation (19) into Equation (18), Equation (20) can be obtained:

$$u\_{\rm o}(t\_{\rm k+1}) = \frac{\rm nD(1-|D|)u\_{\rm dc}(t\_{\rm k})}{2\rm LC f\_{\rm s}^{2}} - \frac{i\_{\rm SC}(t\_{\rm k+1})}{\rm Cf\_{\rm s}} + u\_{\rm o}(t\_{\rm k}).\tag{20}$$

In Equation (20), *u*o(tk) and *u*dc(tk) represent the DAB output voltage and DC bus voltage, respectively, sampled at time *t*k; *u*o(tk+1) represents the predicted value of the SC voltage at time *t*k+1.

According to the general characteristics of the capacitor, the equation of the equivalent capacitor in the SC equivalent model can be obtained, as shown in Equation (21):

$$\mathbf{C}\_{\mathbf{SC}} \frac{\mathbf{d} \mathbf{u}\_{\mathbf{SC}}}{\mathbf{d}t} = \mathbf{i}\_{\mathbf{SC}} \tag{21}$$

where dt = 2Ts in Equation (21). In order to reflect the fast response characteristic of the SC, we chose *i*SC in Equation (21) to be *i*SC(tk+1). Then, Equation (22) can be obtained:

$$[\mathfrak{u}\_{\rm SC}(t\_{\mathbf{k}+1}) - \mathfrak{u}\_{\rm SC}(t\_{\mathbf{k}})] \mathbb{C}\_{\rm SC} = 2\mathcal{T}\_{\rm s} i\_{\rm SC}(t\_{\mathbf{k}+1}) = \frac{i\_{\rm SC}(t\_{\mathbf{k}+1})}{\mathfrak{f}\_{\rm s}}.\tag{22}$$

According to Equation (2), the relationship between the capacitance of the SC and the actual voltage across the SC is as shown in Equation (23):

$$\mathbf{C}\_{\mathbf{SC}} = \mathbf{K}\_{\mathbf{V}} \boldsymbol{\mu}\_{\mathbf{0}} (\mathbf{t}\_{\mathbf{k}}) + \mathbf{C}\_{\mathbf{0}}.\tag{23}$$

According to Equation (22) and Equation (23), taking the average value of each state in the case of the charging and discharging of the SC, Equation (24) can be obtained:

$$\left[\boldsymbol{u}\_{\rm o}(\boldsymbol{t}\_{\rm k+1}) - \boldsymbol{i}\_{\rm SC}(\boldsymbol{t}\_{\rm k+1})\mathbf{R}\_{\rm ES} - \boldsymbol{u}\_{\rm o}(\boldsymbol{t}\_{\rm k}) + \boldsymbol{i}\_{\rm SC}(\boldsymbol{t}\_{\rm k})\mathbf{R}\_{\rm ES}\right] \cdot \left[\mathbf{K}\_{\rm V}\boldsymbol{u}\_{\rm o}(\boldsymbol{t}\_{\rm k}) + \mathbf{C}\_{0}\right] = \frac{\boldsymbol{i}\_{\rm SC}(\boldsymbol{t}\_{\rm k+1})}{\boldsymbol{\mathbf{f}}\_{\rm s}}.\tag{24}$$

By joining Equation (20) and Equation (24) together to eliminate *i*SC(tk), Equation (25) can be obtained:

$$u\_0(t\_{k+1}) = \frac{nD(1-|D|)u\_{\rm dc}(t\_k)}{2\text{L}\text{Cf}\_{\text{s}}^2} \left(1 - \frac{1}{A}\right) - \frac{1}{A}i\_{\rm SC}(t\_k)\mathbf{R}\_{\text{ES}} + u\_0(t\_k),\tag{25}$$

where the value of A is as shown in Equation (26):

$$A = \frac{\mathbb{C}}{\mathcal{K}\_{\rm v} u\_0(t\_k) + \mathcal{C}\_0} + \mathbb{C}f\_{\rm s} \mathcal{R}\_{\rm ES} + 1. \tag{26}$$

### *4.2. Phase Shift Optimization Calculation and Analysis*

For energy storage systems with a DAB, the primary control objective is to stabilize the DC bus voltage. When the bus voltage is too high, the SC needs to be charged. When the bus voltage is low, the SC is required to discharge.

Therefore, different SC reference voltages were selected in this paper. When the bus voltage is high, the SC reference value takes the upper limit warning voltage, Usc\_h; when the bus voltage is low, the SC reference value takes the lower limit warning voltage, Usc\_l. The value of the SC reference value is shown in Equation (27):

$$\mathbf{U}\_{\rm o}{}^{\*} = \begin{cases} \mathbf{U}\_{\rm SC\\_h\nu} \boldsymbol{\mu}\_{\rm dc}(k) > \mathbf{U}\_{\rm dc\\_max} \mathbf{y} \\\ \mathbf{U}\_{\rm SC\\_l\nu} \boldsymbol{\mu}\_{\rm dc}(k) < \mathbf{U}\_{\rm dc\\_min} \mathbf{y} \end{cases} \tag{27}$$

The deformation of Equation (27) is given by Equation (28):

$$\mathbf{U}\_{\rm o}^{\*} = \frac{\mathbf{U}\_{\rm SC\\_h} + \mathbf{U}\_{\rm SC\\_l}}{2} + \frac{\mathbf{U}\_{\rm SC\\_h} - \mathbf{U}\_{\rm SC\\_l}}{2} \cdot B. \tag{28}$$

The value of *B* in Equation (28) is as shown in Equation (29):

$$B = \frac{\frac{\mu\_{\rm dc}(k) - \mathcal{U}\_{\rm dc,max}}{\left| \mu\_{\rm dc}(k) - \mathcal{U}\_{\rm dc,max} \right|} + \frac{\mu\_{\rm dc}(k) - \mathcal{U}\_{\rm dc,min}}{\left| \mu\_{\rm dc}(k) - \mathcal{U}\_{\rm dc,min} \right|}}{2} = \begin{cases} 1, \mu\_{\rm dc}(k) > \mathcal{U}\_{\rm dc,max} \lambda \\\ -1, \mu\_{\rm dc}(k) < \mathcal{U}\_{\rm dc,min} \end{cases} \tag{29}$$

Therefore, a cost function is established based on the prediction model of Equation (26) of the output voltage as shown in Equation (30):

$$J(k) = (\mu\_{\rm o}(t\_{k+1}) - \mathbb{U}\_{\rm o}")^2. \tag{30}$$

In order to minimize the deviation between the output voltage and the reference voltage in the next moment, the cost function needs to be minimized; that is, the selected optimal phase shift amount should minimize the cost function value. By combining Equations (25), (28), and (29) and deriving Equation (30), the value of *D* can be obtained as in Equation (31):

$$D = \mathbf{B} \cdot \frac{1 - \sqrt{1 - 4 \mathbf{B} \mathbf{M}}}{2}. \tag{31}$$

The value of *M* in Equation (31) is as shown in Equation (32):

$$M = \frac{A}{A - 1} \cdot \frac{2 \text{L} \text{Cf}\_{\text{s}}^{2}}{\text{nu}\_{\text{dc}}(k)} [\text{U}\_{\text{o}}^{\*} - \text{u}\_{\text{o}}(t\_{\text{k}}) + \frac{1}{A} \cdot i\_{\text{SC}}(t\_{\text{k}}) \text{R}\_{\text{ES}}]. \tag{32}$$

A block diagram of the MPC structure constructed according to Equation (32) is shown in Figure 5.

**Figure 5.** MPC structure block diagram.

It can be seen from the above analysis that the optimal phase shift of the DAB converter in the proposed method is obtained based on the sampling voltage, sampling current, and circuit parameters of the DAB-SC system. Additionally, the parameters of the SC can reflect the fast response characteristics of the SC. Compared with predictive control based on the simple construction of a DAB converter, the proposed control strategy makes the overall response of the DAB-SC system faster.

### **5. Simulation Analysis**

A simulation model of the DAB-SC system, as shown in Figure 6, was built in MATLAB/simulink according to the topology diagram of Figure 2 and the control block diagram of Figure 5. The parameters of the simulation are shown in Table 1. Additionally, the experimental simulation pulse signal generator had a frequency of 20 kHz and a duty cycle of 50%.

**Figure 6.** Simulation model of the DAB-SC system based on Matlab/Simulink: (**a**) The overall model; (**b**) model of the DAB converter.

**Table 1.** Simulation Parameters.

### *5.1. SC Discharge in the Normal Working Range*

When the DC bus voltage, *U*dc, is lower than the lower limit voltage of 665 V, which is 651.7 V, and the initial value of the voltage, *U*o, on the SC side is 75 V, the dynamic response time of the system is 1.9 ms when the bus voltage is stabilized to 665 V with only a prediction model of the DAB converter. The simulation waveform is shown in Figure 7a. Under the same conditions, using the

overall prediction model composed of the DAB converter and the SC to control the SC discharge to stabilize the bus voltage to 665 V, the dynamic response time of the system is 1.5 ms, and the simulation waveform is shown in Figure 7b. Under these simulation conditions, a detailed waveform of the drive signal is shown in Figure A1 (see Appendix A).

**Figure 7.** Simulation waveform diagram when the initial voltage, *U*dc, is 651.7 V and the initial voltage, *U*o, is 75 V: (**a**) MPC only for the DAB; (**b**) MPC for the DAB-SC system.

When the DC bus voltage, *U*dc, is lower than the lower limit voltage of 665 V, which is 651.7 V, and the initial value of the voltage, *U*o, on the SC side is 58.6 V, the dynamic response time of the system is 24.5 ms when the bus voltage is stabilized to 665 V by using the DAB converter only to establish a prediction model to control the SC discharge. The simulation waveform is shown in Figure 8a. Under the same conditions, using the overall prediction model composed of the DAB converter and the SC to control the SC discharge to stabilize the bus voltage to 665 V, the dynamic response time of the system is 7.5 ms, and the simulation waveform is shown in Figure 8b.

It can be seen from the simulation results in Figures 7 and 8 that when the SC is discharged from different SOCs, the method of using the overall model of the DAB converter and the SC for predictive control to stabilize the bus voltage always obtains faster response times than the method of using the MPC only for the DAB converter to stabilize the bus voltage.

**Figure 8.** Simulation waveform diagram when the initial voltage, *U*dc, is 651.7 V and the initial voltage, *U*o, is 58.6 V: (**a**) MPC only for the DAB; (**b**) MPC for the DAB-SC system.

### *5.2. SC Charge in the Normal Working Range*

When the DC bus voltage, *U*dc, is higher than the higher limit voltage of 735 V, which is 749.7 V, and the initial value of the voltage, *U*o, on the SC side is 75 V, the dynamic response time of the system is 1.7 ms when the bus voltage is stabilized to 735 V by using the DAB converter only to establish a prediction model to control the SC charge. The simulation waveform is shown in Figure 9a. Under the same conditions, using the overall prediction model composed of the DAB converter and the SC to control the SC charge to stabilize the bus voltage to 735 V, the dynamic response time of the system is 1.1 ms, and the simulation waveform is shown in Figure 9b. Under the simulation conditions, a detailed waveform of the drive signal is shown in Figure A2 (see Appendix B).

When the output power of the photovoltaic is high, the voltage, *U*dc, on the DC bus side is higher than the higher limit voltage of 735 V, which is 749.7 V, and the initial value of the voltage, *U*o, on the SC side is 93.6 V, the dynamic response time of the system is 14 ms when the bus voltage is stabilized to 735 V by using the DAB converter only to establish a prediction model to control the SC charge. The simulation waveform is shown in Figure 10a. Under the same conditions, using the overall prediction model composed of the DAB converter and the SC to control the SC charge to stabilize the bus voltage to 735 V, the dynamic response time of the system is 5.9 ms, and the simulation waveform is shown in Figure 10b.

It can be seen from the simulation results in Figures 9 and 10 that when the SC is charged from different SOCs, the method of using the overall model of the DAB converter and the SC for predictive control to stabilize the bus voltage always obtains faster response times than the method of using the MPC only for the DAB converter to stabilize the bus voltage.

**Figure 9.** Simulation waveform diagram when the initial voltage, *U*dc, is 749.7 V and the initial voltage, *U*o, is 75 V: (**a**) MPC only for the DAB; (**b**) MPC for the DAB-SC system.

**Figure 10.** Simulation waveform diagram when the initial voltage, *U*dc, is 749.7 V and the initial voltage, *U*o, is 93.6 V: (**a**) MPC only for the DAB; (**b**) MPC for the DAB-SC system.

### *5.3. SC Discharge Closer to the Lower Limit*

Through the observations of Figures 7–10, it can be found that in these four cases, compared with the method of using the MPC only for the DAB converter to stabilize the bus voltage, the method of using the overall prediction model composed of the DAB converter and the SC to stabilize the bus voltage makes the SC voltage have a larger amount of change. Figure 11 contains the simulation waveforms obtained by using two different control methods when the voltage, *U*dc, on the DC bus side is 651.7 V and the initial value of the voltage, *U*o, on the SC side is 58 V. As can be seen from Figure 11a, in this case, a method of establishing a prediction model only for the DAB converter to control the discharge of the SC cannot stabilize the bus voltage to 665 V. It can be seen from Figure 11b that the dynamic response time of the system is 12 ms by using the overall prediction model composed of the DAB converter and the SC to control the SC discharge to stabilize the bus voltage to 665 V.

**Figure 11.** Simulation waveform diagram when the initial voltage, *U*dc, is 651.7 V and the initial voltage, *U*o, is 58 V: (**a**) MPC only for the DAB; (**b**) MPC for the DAB-SC system.

When the voltage, *U*dc, on the DC bus side is 651.7 V, and the initial value of the voltage, *U*o, on the SC side is 56.58 V, using the overall prediction model composed of the DAB converter and the SC to control the SC charge to stabilize the bus voltage to 665 V, the dynamic response time of the system is 15 ms, and the simulation waveform is shown in Figure 12.

**Figure 12.** Simulation waveform diagram of the MPC for the DAB-SC system when the initial voltage, *U*dc, is 651.7 V and the initial voltage, *U*o, is 56.58 V.

### *5.4. SC Charge Closer to the Upper Limit*

Figure 13 contains the simulation waveforms obtained by using two different control methods when the voltage, *U*dc, on the DC bus side is 749.7 V and the initial value of the voltage, *U*o, on the SC side is 94 V. As can be seen from Figure 13a, in this case, the method of establishing a prediction model only for the DAB converter to control the charge of the SC cannot stabilize the bus voltage to 735 V. It can be seen from Figure 13b that the dynamic response time of the system is 6 ms by using the overall prediction model composed of the DAB converter and the SC to control the SC charge to stabilize the bus voltage to 735 V.

When the voltage, *U*dc, on the DC bus side is 749.7 V, and the initial value of the voltage, *U*o, on the SC side is 98.1 V, using the overall prediction model composed of the DAB converter and the SC to control the SC charge to stabilize the bus voltage to 735 V, the dynamic response time of the system is 22.3 ms, and the simulation waveform is shown in Figure 14.

It can be seen from the simulation results of Figures 11–14 that compared with the method of establishing a prediction model only for the DAB converter to stabilize the bus voltage, the method of using the overall prediction model composed of the DAB converter and the SC to stabilize the bus voltage increases the actual upper limit of the SC voltage, reduces the actual lower limit of the SC voltage, and expands the actual capacity utilization range of the SC by 18.63%.

*Energies* **2019**, *12*, 2134

**Figure 13.** Experimental waveform diagram when the initial voltage, *U*dc, is 749.7 V and the initial voltage, *U*o, is 94 V: (**a**) MPC only for the DAB; (**b**) MPC for the DAB-SC system.

**Figure 14.** Simulation waveform diagram of the MPC for the DAB-SC system when the initial voltage, *U*dc, is 749.7 V and the initial voltage, uo, is 98.1 V.

### **6. Conclusions**

In this paper, a hybrid MPC strategy composed of a DAB and SC was proposed in order to improve the dynamic characteristics of an energy storage system. The structural characteristics of the SCs were analyzed, and a suitable equivalent model of the SC was chosen for its fast response characteristics. The working principle of the DAB converter in each stage under single-phase shift control was analyzed and the prediction model of the DAB-SC system was established. By combining the SC model and the DAB model, an appropriate evaluation target was selected to establish the objective function for model prediction. A simulation model of the DAB-SC system was built in MATLAB/Simulink for verification. The simulation results show that the proposed method with the overall prediction model composed of a DAB converter and SC can effectively improve the overall response speed of the system and expand the actual capacity utilization range of the SC. The validity of the proposed method was verified and the proposed method has good application prospects in improving the dynamic response performance of energy storage systems.

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

**Funding:** This research was funded by National Natural Science Foundation Youth Project of China, grant number 51607060.

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

### **Abbreviation**

### **Appendix A**

It can be seen from Figure A1 that the secondary side drive signal of the converter leads the primary side drive signal in SC discharge state. And comparing Figures A1a and A1b, it can be seen that compared with the method of establishing a prediction model only for the DAB converter to control SC discharge, the method of controlling SC discharge by using the overall prediction model composed of the DAB converter and SC has a larger amplitude of inductor current, and due to the internal resistance RES of SC, the voltage on the low voltage side declines instantaneously when SC is discharged, making *U*dc > n*U*o, therefore, the inductor current in the power transfer mode of the converter decreases linearly when SC is discharged.

**Figure A1.** Detailed simulation waveform diagram of the drive signal when *U*dc initial voltage is 651.7 V and *U*o initial voltage is 75 V: (**a**) MPC only for DAB; (**b**) MPC for DAB-SC system.

### **Appendix B**

It can be seen from Figure A2 that the primary side drive signal of the converter leads the secondary side drive signal in SC discharge state. And comparing Figures A2a and A2b, it can be seen that compared with the method of establishing a prediction model only for the DAB converter to control SC discharge, the method of controlling SC discharge by using the overall prediction model composed of the DAB converter and SC has a larger amplitude of the inductor current, and due to the internal resistance RES of SC, the voltage on the low voltage side rises instantaneously when SC is charged, making *U*dc < n*U*o, therefore, the inductor current in the power transfer mode of the converter decreases linearly when SC is charged.

**Figure A2.** Detailed Simulation waveform diagram of the drive signal when *U*dc initial voltage is 749.7 V and *U*o initial voltage is 75 V: (**a**) MPC only for DAB; (**b**) MPC for DAB-SC system.

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

### **Design and Validation of Ultra-Fast Charging Infrastructures Based on Supercapacitors for Urban Public Transportation Applications**

**Fernando Ortenzi 1,\* , Manlio Pasquali 1, Pier Paolo Prosini 1, Alessandro Lidozzi 2 and Marco Di Benedetto 2**

Received: 20 May 2019; Accepted: 17 June 2019; Published: 19 June 2019

**Abstract:** The last few decades have seen a significant increase in the number of electric vehicles (EVs) for private and public transportation around the world. This has resulted in high power demands on the electrical grid, especially when fast and ultra-fast or flash (at the bus-stop) charging are required. Consequently, a ground storage should be used in order to mitigate the peak power request period. This paper deals with an innovative and simple fast charging infrastructure based on supercapacitors, used to charge the energy storage system on board electric buses. According to the charging level of the electric bus, the proposed fast charging system is able to provide the maximum power of 180 kW without exceeding 30 s and without using DC–DC converters. In order to limit the maximum charging current, the electric bus is charged in three steps through three different connectors placed between the supercapacitors on board the bus and the fast charging system. The fast charging system has been carefully designed, taking into account several system parameters, such as charging time, maximum current, and voltage. Experimental tests have been performed on a fast charging station prototype to validate the theoretical analysis and functionality of the proposed architecture.

**Keywords:** supercapacitors; fast and flash charging; urban public transportation

### **1. Introduction**

The electric vehicles (EVs) in public transportation could be the keystone to limit a city's traffic problem, noise level, and pollution within urban complexes with high population density [1,2]. In fact, having electric public vehicles in the city center will have a large positive environmental effect on both noise and pollution. According to the many advantages provided by EVs, more and more cities have been replacing their public vehicles with electric vehicles [3–8]. In [7], a cost-benefit method was employed for the scheduling of an electric city bus fleet on a single route. Three different charging infrastructure scenarios were considered. In the first scenario, only one charging station was used. The second scenario considered two charging stations that were located at the same terminus. In the third scenario, two charging stations were located at opposite terminuses. The first scenario is the cheaper solution compared to the others. The second solution is more suitable in terms of energy demand and traffic conditions. However, the electric vehicles, such as electric buses, are limited by the onboard storage system and the long charging times, which depends on the charging infrastructure. There are several different methods to charge the electric buses [8,9].

1. Overnight charging method. In this case, the electric buses are equipped with a huge energy storage system in order to satisfy the bus routes during the day. The large storage system (i.e.,

for a 12 m bus there are about 350 kWh of lithium batteries) must be recharged during the night (in about 6 h).

Compared to the slow charging stations, the main advantages of the fast charging stations are an extended range autonomy and a smaller storage system. Unfortunately, the connection of the unexpected loads, such as an energy storage system, to the electrical grid may cause unwanted effects on the grid due to the great power demands for short times. For instance, the Hess-ABB Tosa (Trolleybus Optimisation System Alimentation), an 18 m bus fitted out by ABB with an optimized complete traction chain, is composed by 40 kWh energy-storage units and an automatic energy transfer system at bus-stops. This electric bus is charged two times on the 18 km route. The first charging point, called Palexpo Flash Station, is equipped with a 400 kW power supply unit and the charging time is only 15 s, whereas the second charging point is composed by a 200 kW power supply unit, and the charging time is 3–4 min. These storage systems inevitably require expensive components with a large effect on investment costs. Performance and cost optimization of charging stations and electric storage systems pass through technology solutions able to integrate devices with different functions.

In [1], the evaluating methodology of power-quality field measurements of an electric bus charging station in Taiwan was proposed. In this case, the charging station under study had 11 charging poles, including two fast charging and nine slow charging poles. Experimental measurements of voltage flickers, harmonic distortion, including electrical power measurements, and three phase voltages were performed. In [2], a control structure for the DC charging station was proposed. The charging station makes use of the frontend AC–DC converter, including the power grid filters, and backend DC–DC converter connected via DC link capacitor. The DC–DC converter is controlled through the P-controller, whereas the AC–DC converter is controlled using the PI (Proportional-Integral) controller, where the input error is provided by the difference between the measured grid voltage and reference voltages. The charging station proposed in [5] is equipped with a modular multilevel converter (MMC) and the control structure is based on finite control set-model predictive control.

Such charging processes can be realized using the supercapacitors. The advantages of supercapacitors can be highlighted in three different characteristics [10]: (1) high specific power—supercapacitors can be totally discharged and charged in the order of seconds; (2) long life, to the order of 106 cycles; (3) low specific energy.

The use of supercapacitors in both hybrid vehicles and hybrid storage configurations, including microcars or forklifts, has been discussed in different papers [11–13]. In an electric microcar tested on a dynamometer chassis running the urban part of the previous type approval procedure driving cycle (ECE15) [11], the addition of supercapacitors to the battery storage brought a reduction of the maximum current from lead-acid batteries and then advantages in terms of vehicle range (more than 50%) and battery life (more than double). The same results have been obtained also on a forklift [12]. Such advantages (in terms of range) have not been obtained if the vehicle is equipped with lithium-ion batteries due to their lower internal resistance [13].

In this paper, a flash charging infrastructure, a three-year project called "Ricerca di Sistema Elettrico", has been designed, evaluated in terms of investment costs [14,15], and tested using a reduced scale model prototype [16]. Such charging infrastructure needs to be equipped with a storage system based on supercapacitors due to the high-power transfer. The realized electric bus has a hybrid storage system composed of lead acid batteries for the auxiliaries and backup, as well as the supercapacitors for the traction. Furthermore, the proposed system presents an innovative concept, where the charging is directly performed between two supercapacitors without use of the converter: the first one placed in the charging station and the second one placed on board the electric bus. To limit the charging current, an inductor is connected between the supercapacitors and the charge is made in three different configurations by changing the connection of the supercapacitors in the charging station. The advantages of the present configuration are that the charging station is simpler and does not have any expensive devices like a converter; the only components used are contactors and an inductance, cheaper than a DC–DC converter designed for such a purpose.

The paper is organized in six parts. The background of the charging station in the electric buses is discussed in Section 1. The electric bus details (Tecnobus Gulliver U500), equipped with the pantograph and the hybrid storage system, adapted to be compliant with the present charging system, is explained in Section 2. The charging architecture and the charging control strategy is discussed in Section 3, while experimental results are presented in Section 4. Section 5 investigates the possible structure and the management of a complex charging network for public transport. Finally, Section 6 gives the conclusions.

### **2. The Pantograph Equipped Vehicle**

The vehicle used to test the flash charge infrastructure was a minibus, Tecnobus Gulliver U500, also used for previous fast charge experimentation [17,18]. The specifications of the minibus are reported in Table 1.

**Table 1.** Tecnobus Gulliver U500 Specification.

This vehicle was originally equipped with lead acid batteries. The previous experimental results, carried out at the ENEA (Italian Agency for New Technologies, Energy and Sustainable Economic Development) Casaccia Research Center, were also obtained using the LFP (lithium iron phosphate ) batteries [18].

The energy storage system on board the bus was specifically designed in order to use the flash charging method. The energy storage was composed by a hybrid system, where 8.6 kWh lead acid batteries (72 V/120 Ah) were used together with three Maxwell supercapacitor modules (125 V, 63 F) connected in series, as shown in Figure 1a. The lead acid batteries were used for both the auxiliaries and as a backup battery, when the supercapacitors were totally discharged during the route of the bus. Furthermore, the batteries helped to partially discharge the supercapacitors, when their voltage was too high to be charged. The nominal energy of the supercapacitors was 410 Wh; nevertheless, since the

supercapacitors worked at half their nominal voltage, only 75% of the supercapacitors' nominal energy was used. The description of the architecture on board the bus is explained in more detail in [19]. Figure 1b shows the pantograph installed on the bus. The pantograph was a fast charging device manufactured by Schunk with the following characteristics:

**Figure 1.** (**a**) Hybrid Storage system on board the bus; (**b**) and the pantograph for ultrafast charge.

### **3. Charging Architecture**

The main characteristics of the flash charging method proposed in this paper are the following:

In [14–16] the charging station was designed and tested in a reduced scale prototype. In that case, the power levels were not taken into account and the electric scheme of the charging station was created to optimize the simulation model. Such a system was composed by six ultra-capacitors in the ground charging station and three onboard the vehicle, each one of 16 V and 250 F. Starting from this analysis, in this paper a real scale protype was built. In this paper, the proposed charging station was built with the architecture explained below.

The block scheme of the proposed charging station is depicted in Figure 2. Five contactors were used to operate the charging process. The charging was made in three steps: (1) four modules, connected two in series and two in parallel (2s2p), were used to charge the vehicle; (2) an extra module was added in series to the others; and (3) a sixth module was connected in series to the others.

**Figure 2.** Block scheme of the charging station.

In the first step, the contactors RL5, RL11, and RL7 were closed to allow the four-supercapacitor module to transfer energy from the charging station to the onboard supercapacitors. When the charging transfer was almost complete (or the current fell below 20 A), the contactor RL5 was opened, while the contactor RL3 was closed to start the second step. When the second step was finished, the closure of the RL1 started the third step. Finally, four AC–DC converters with 2 kW rated power were used to charge the supercapacitors after the energy transfer process. The maximum charging current was limited, using the inductor and resistor with the parameters listed in Table 2. The inductor is shown in Figure 3. The supercapacitor configuration together with the inductance is an efficient strategy to manage the flash charge in a cheap way. The reason for the use of three recharge steps instead of one is related to the possibility of reducing the maximum recharge current and optimizing the energy transfer.

**Figure 3.** Inductor used in the charging station.

**Table 2.** Inductor and resistor parameters.

L = 5.4 mH Resistance = 0.155 Ω Maximum Current = 500 A Ku = 0.5 IRMS = 187 A

### **4. Results**

Different tests were carried out in order to validate the design of the proposed charging station.

The simulation results performed in one step (1-step) and three steps (3-step) recharge processes for the proposed system are shown in Figures 4 and 5. The 1-step charging is faster, but the 3-step charging is more efficient and the correspondent current is lower. Table 3 shows the comparison between the two charging methods.

**Figure 4.** Comparison of 1-step and 3-step charge. (**a**) Charging current as a function of time, (**b**) bus supercapacitors voltage.

Δ

**Figure 5.** Voltages and currents as a function of time: (**a**) station side, (**b**) onboard side.

**Table 3.** Comparison of 1 and 3 step recharge.

It has been found that the test results depend on the initial voltage of the supercapacitors. In Table 2 the extreme conditions of recharge are considered. The voltage and current profile during a charging are shown in Figure 5. The six modules connected as in Figure 2 present different voltage values during the charge. It can be seen from Figure 5 that the three steps can be easily identified, observing three different peak current and the three different voltages (supercapacitor voltages on board the vehicle and the voltage on the ground station).

The maximum peak current, more than 400 A, was obtained during the second step and the supercapacitor module connected in the second step reached the lowest value of voltage at the end of the charge process. The charging time, without the delay times due to the manual drive of the charging steps, was within the 30 s. The charging time of the ground supercapacitors was about 2.5 min.

The charging process when the initial voltage of the onboard supercapacitors was changed from 190 V to 250 V is depicted in Figures 6 and 7. With reference to different onboard voltages, the transfer energy was reduced during the first step due to the different voltage between ground and onboard storage. The current changed from 300 A (for voltage = 190 V) to 220 A (voltage = 230 V). Additionally, the first current step is missing when the onboard voltage is 250 V (Figure 7b).

**Figure 6.** Charging process with on board initial voltage of: (**a**) 190 V, (**b**) 210V.

**Figure 7.** Charge process with onboard initial voltage of: (**a**) 230 V, (**b**) 250V.

Thanks to the lighter storage system compared to the original version (powered with 43 kWh lead acid batteries) the unladen weight passed from 4300 kg to about 3500 kg and the gross laden weight from 6300 to 5600 kg. Consequently, the range on a reference circuit, internal to the Enea Casaccia Research Center of about 0.7 km, may range from 320 Wh/km up to 450 Wh/km, depending on the number of passengers transported.

The bus route was performed in different conditions: a full and empty bus. Moreover, the state of charge (SOC) of the supercapacitors was changed from 200 V to 275 V. Figure 8 shows the reference speed profile and the SOC for different bus conditions and different voltages of the supercapacitor tank. As can be seen from Figure 8, the charging station is able to charge the onboard supercapacitors independently of their initial voltage.

**Figure 8.** (**a**) Vehicle speed profile and state of charge (SOC) for different bus weights, (**b**) voltage of the supercapacitor tank and related SOC.

Figure 9 shows the final voltage at the end of the charge process as a function of the initial voltage of the supercapacitors. When the onboard SOC is close to zero, the final voltage is about 360 V, and it increases with an increase in the initial voltage. When the initial voltage reaches the maximum value (250 V), the charging station uses the two charge steps to bring the final voltage at 375 V. This is the rated voltage for the supercapacitors; it is not possible to charge the energy storage system of the bus if the initial voltage has higher values. However, the bus has a hybrid storage system composed by batteries and supercapacitors. Some of the energy from the supercapacitors is transferred to the batteries in order to decrease the supercapacitor voltage. Thus, the battery storage on board the bus can be used to discharge the supercapacitors when their voltage is too high. Furthermore, battery storage is used as a backup system when the supercapacitors are totally discharged during the route.

**Figure 9.** Final onboard supercapacitor voltage as a function of the initial voltage (as recorded before the charging process).

### **5. Structure of a Charging Network and its Management**

In this paper, the structure of a flash charging station in the bus route was considered. Once the feasibility of the idea is verified, it is necessary to consider the charging network building for a complex transport system.

Figure 10 shows how it is possible to connect a supercapacitor bank to different bus stops, and it is also possible to connect a bus stop to different supercapacitor banks.

**Figure 10.** Proposed structure for the flash charging network.

In this way, it is possible to create a structured flash recharging network. To build the network the following must be considered:

For example, if two stops are very close, and a large interval of time is expected for the passage of the buses, only one bank of capacitors can be used; if instead there are many stops nearby with a high frequency of passing buses it is necessary to provide more banks of capacitors, greater than the number of stops if needed.

A rule that has been established is that if two stops are more than a hundred meters apart then they are not connected.

To manage the network, it is necessary to establish the value of the supercapacitors' charge power, Figure 11 shows a possible technique for managing the charging power.

**Figure 11.** Proposed criteria for network management.

From the knowledge of waiting times at the bus stop, a solution for choosing the capacitor recharge power is shown in Figure 11.

The subject of study has been introduced for completeness; it will be explained in more detail in other papers, but to give an example, consider a simplified case.

Consider the following charging network (Figure 12), which can be seen as a system where the status is given by the position of the buses and the SOC of the capacitors on the ground. The input is the direction and speed of the buses, or equivalent, the waiting times at the stops and the output are the SuperCapacitorsconnections with the bus stops and their recharging powers.

**Figure 12.** Example of flash charge network.

In a first example, the waiting time at the stops is considered a multiple of 3 min, the supercapacitors' recharge power is a constant and equal to 8 kW. In another case, it is possible to choose between three thresholds—8, 4, and 2.7 kW.

The following Table 4 shows the transition times at the stops for four buses.

**Table 4.** Bus travel times and relative stops.

Table 5 indicates a possible management of the capacitors by adopting a single recharge threshold; the term Cif indicates that the capacitor is connected to the station f.

**Table 5.** Management for a single recharge threshold, Power = 8 KW.

Table 6 indicates a possible solution with three recharge thresholds.

**Table 6.** Management with three recharge thresholds, P = 2.7,4, and 8 KW.

Figure 13 allows the results to be compared for the same time interval.

The subject of study was introduced for completeness; it will be explained in more detail in the other papers.

15 0 0 0 0 0 0 0 0 0 0

**Figure 13.** Bus position and configuration of the charging network for 6' < T < 9 in the two solutions.

### **6. Conclusions**

The preliminary results of a flash charge system for a public transport vehicle characterized by supercapacitors and batteries storage system onboard the vehicle were presented. The charging process was divided into three steps in order to limit the maximum currents.

It can be noticed that such a charging system could be effective from both a technical and economical point of view when used for reduced distance between two adjacent stops. The necessary grid average power can be definitely fulfilled by a common urban light system due to the quite low average value. In fact, the charging process of the ground station is accomplished with a longer time between two different vehicles. The results have confirmed the theoretical, modeling, and experimental study served for sizing and implementing of the 1:1 scale system.

In the final part of the paper, suggestions have been proposed for future improvements of the charging system and for alternative structures of the flash charging stations.

**Author Contributions:** Funding acquisition, P.P.P.; project administration, F.O.; software, A.L., M.D.B.; writing—original draft preparation, M.P.; writing—review and editing, F.O., M.P and A.L.

**Funding:** This research was funded by Ministero dello Sviluppo Economico, Project C.5, Energy Storage systems for the electric system, PAR 2017.

**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* **Fractional-Order Modeling and Parameter Identification for Ultracapacitors with a New Hybrid SOA Method**

### **Jianhua Guo 1, Weilun Liu 1,\*, Liang Chu 1 and Jingyuan Zhao 2**

Received: 10 October 2019; Accepted: 7 November 2019; Published: 8 November 2019

**Abstract:** This paper deals with an ultracapacitor (UC) model and its identification procedure. To take UC's fractional characteristic into account, two constant phase elements (CPEs) are used to construct a model structure according to impedance spectrum analysis. The different behaviors of UC such as capacitance, resistance, and charge distribution dynamics are simulated by the corresponding part in the model. The resistance under different voltages is calculated through the voltage rebound method to explore its non-linear characteristics and create a look-up table. A nonlinear fractional model around an operation voltage is then deduced by applying the resistance table. This time identification is carried by a proposed hybrid optimization algorithm: Nelder-Mead seeker algorithm (NMSA), which embeds the Nelder–Mead Simplex (NMS) method into the seeker optimization algorithm (SOA). Its time behavior has been compared with the linear fractional model for charging and discharging current profiles at different levels.

**Keywords:** hybrid optimization; ultracapacitor; fractional-order model; time-domain analysis

### **1. Introduction**

Ultracapacitors (UCs) are currently gaining attention rapidly as a new type of energy-storage component with a higher density than ordinary capacitance and longer cycle life than batteries [1,2]. However, its energy storage principle determines its capacity cannot temporarily surpass the battery [3]. Thence, it is always associated with batteries in the energy storage unit when it supplies energy to electric vehicles (EV) or hybrid electric vehicles (HEV). UC can meet vehicle startup and acceleration requirements due to its high power rate. Thanks to the short charging time of UC, the braking energy can be absorbed efficiently during decelerating. With a reasonable power distribution strategy in-vehicle, the application of UC is a good way to reduce battery load and further to improve the battery's cycle life.

In the early stage of UC modeling, the internal properties of UC are analyzed in many electrochemical models. In their model, Gouy and Chapman adopted the principle that ion distribution of the double-layer's outer element conformed to the laws of Bolzmann [4]. The mathematical model proposed in [5] described the effect of pore size on UC performance. On the contrary, the machine learning model describes the relationships between inputs and outputs without considering the internal structure. The subject [6–8] utilized ANN technology to predict the performance of carbon materials in supercapacitors. Among all the modeling approaches, the equivalent circuit model is the most commonly used. Through studying the slow charging behavior of UC, R.L. Spyker and R.M. Nelms used a classical equivalent circuit composed of capacitance in parallel with resistance and an equivalent series resistance [9]. To take into account the fast dynamics, a model consisting of three RC

branches is proposed that reflected the internal charge distribution process within a specific voltage range [10]. Because of the widespread use of activated carbon with high surface area, some groups adopted the transmission line model (TLM) to simulate the frequency response of UC and discussed the contribution of the parameters to the impedance spectrum [11]. The complex determination of different elements and large numbers of parameters leads to inconvenience.

To simplify the complexity of the model, the fractional-order model has been employed. A simple UC model composed of a series resistor and a CPE is used to model the impedance in the frequency range in [12]. The capacitance-voltage dependence is taken into account in [13] by implanting a voltage-dependent capacitance in one of the three branches and a fitting formula with an undetermined coefficient is used. N. Bertrand et al. [14] used a set of linear systems to deduce a global nonlinear model.

In a process of model identification, there are traditional methods such as the least square methods, the gradient-based methods, and the maximum likelihood methods. They can be suitable for nonlinear systems after some efforts made by [15]. However, the recognition process of all the above methods needs to derive the objective function, which is difficult to implement in fractional-order systems. The derivative-free algorithms, such as the simplex search method [16] and Nelder–Mead simplex (NMS) algorithm, [17] are successfully applied to the parameter tuning. The limitations of these techniques are the natural local concern and the sensitivity to the initial solution of the algorithm. Fractional order system identification problems have been carried out using global heuristic techniques such as particle swarm optimization (PSO) [18], genetic algorithm (GA) [19], and artificial bee colony algorithm (ABC) [20]. Compared to the local search algorithm, they have a global performance nature and a greater possibility to converge to the global optimal point. Nevertheless, they often oscillate in the vicinity of the optimal point or miss the optimal area due to strong randomness.

The seeker optimization algorithm (SOA) is also a heuristic algorithm, which was invented by Dai [21]. It explores the experience of human search behavior through the analysis of the human random search process. By exchanging experience, seekers determine their direction and step size, then update their location, and, finally, find the optimal point. Since it is proposed, it has been used by researchers in different engineering applications like system modeling and optimal control [22,23]. To overcome the individual limits of SOA, a hybridization of SOA and NM algorithm is used in this paper, which simultaneously retains the benefits of both techniques while discarding their drawbacks.

The proposed modeling approach leads to a model that has the following characteristics:

### **2. Fractional-Order Calculus**

Fractional-order calculus (FOC) generalizes the integrals and derivatives from integer-order to arbitrary real order. This theory goes back to a note written by Leibniz, in which the meaning of the half-derivative of the function *x* was discussed. Some approaches to generalizations of the notion of calculus can be found in [24]. A great deal of work has been done to propose the models of linear time-invariant system and nonlinear system for analyzing the stability, time, frequency domain response, and other properties of the fractional-order model [25,26]. FOC provides an excellent instrument for the description of stochastic memory phenomena such as electrochemical phenomena, chaos, and viscoelasticity. Consider a function *y* = *f*(*t*), which is integrable and has at least *m*(*m*-*N*) continuous derivatives. The universal expression of integer-order derivatives and integrals is as follows:

$$\,\_{a}D\_{t}^{p}f(t) = \lim\_{h \to 0} h^{-p} \sum\_{r=0}^{n} (-1)^{r} \binom{p}{r} f(t - rh), \tag{1}$$

$$
\binom{p}{r} = \frac{p(p-1)(p-2)\cdots(p-r+1)}{r!} \tag{2}
$$

 *p r* is the usual notation for the binomial coefficient, *n* = [(*t* − *a*)/*h*]. For *p* = *m*(*m* > 0), *aDp t* = *dmt dtm* represents the derivative of order *m*. On the contrary, *aDp t* denotes *m* fold integral if *p* = −*m*(*m* > 0). The unification of integer-order derivatives and integrals is the basis for allowing *p* in (1) to be an arbitrary real number or even a complex number. In this section, three main definitions of fractional order integrals and derivatives are briefly introduced.

• *Grünwald–Letnikov Definition:*

$$\,\_{a}D\_{t}^{a}f(t) = \lim\_{h \to 0} h^{-a} \sum\_{r=0}^{n} (-1)^{r} \binom{a}{r} f(t - rh) \,\_{\ast} (m - 1 < a < m) \tag{3}$$

$$
\sigma \left( \begin{array}{c} \alpha \\ r \end{array} \right) = \frac{1}{(\alpha+1)B(\alpha-r+1, r+1)} \tag{4}
$$

• *Riemann–Liouville Definition:*

$$\,\_aD\_t^a f(t) = \frac{1}{\Gamma(m-a)} \frac{d^m}{dt} \int\_a^t (t-\tau)^{m-a-1} f(\tau) d\tau,\\ (m-1 < a < m) \tag{5}$$

• *Caputo definition:*

$$\,\_aD\_t^a f(t) = \frac{1}{\Gamma(m-\alpha)} \int\_a^t (t-\tau)^{m-\alpha-1} f^{(m)}(\tau) d\tau, \,(m-1 < a < m) \tag{6}$$

Grünwald–Letnikov generalizes integer-order calculus directly and requires that *f*(*t*) can be differentiated *m* = [α] +1 times (α ∈ *R*). However, a limit operation is not convenient. Rieman–Liouville provides an excellent opportunity to weaken the requirement for the differentiability of *f*(*t*) through *m*-α-1 times integral before *m*-order derivative. Another definition was proposed by M.Caputo that exchanges the sequence of derivative and integral. The main advantage of the Caputo approach is permitting the initial conditions for fractional differential equations to take on the same form as for integer-order differential equations, which makes them easy to interpret. Another difference between Caputo's definition and the other two definitions can be shown in the derivative of a constant C. The Caputo derivative of a constant is 0, whereas the other two bring a nonzero value as Equation (7), so the Caputo definition is used in this paper.

$$\, \_aD\_t^{\alpha} \mathbb{C} = \frac{\mathbb{C}t^{-\alpha}}{\Gamma(1-\alpha)}\tag{7}$$

### **3. A New Hybrid Algorithm: NMSA**

### *3.1. Seeker Optimization Algorithm (SOA)*

SOA is a new swarm intelligence optimization algorithm, which has better optimization quality and robustness than GA and PSO [21]. It resolves optimization problems by means of map search. *n* seekers search for the lowest position on the *D*-dimension map if there are *D* optimized variables. The lowest position stands for the set of optimal parameters and the height of a point is calculated by the objective function. Their journey lasts for *g* days, to begin with, they are sent to the initial position randomly on the first day.  Meanwhile, they share their locations and altitudes

with the information network, which can be read by everyone and the best location of the past few days can be selected. Afterward, according to several reference locations, each seeker determines the direction and distance to go tomorrow, following a certain rule, and start their search the next day. They move, communicate, and make plans at the end of the day journey until their whole journey is finished and a position closest to the bottom is found.

On *k*th day, *i*th seeker *pi*(*k*) has a tendency to move towards the personal best position *pibest*(*k*) and the best position of population *pzbest*(*k*). Meanwhile, he tends to inherit the moving direction *di*(*k* − 1) of the yesterday if *f*(*pi*(*k*)) < *f*(*pi*(*k* − 1)) and to reverse it if not. Another reference direction might also exist, which is represented by *diother*(*k*). Those reference directions are marked in Figure 1. Using horizontal as an example, after the number of reference directions toward the positive direction and those toward the negative direction are counted, respectively, draw the probability distribution semicircle. Then *i*th seeker determines his moving direction *di*(*k*) by generating a random angle. If the angle is as shown in Figure 1, the *i*th seeker moves positively in horizontal.

**Figure 1.** Determine the search direction.

The search step length *li*(*k*) is calculated by Gaussian membership function in Equation (8). μ*i* is a decreasing function of *i*th seeker's fitness ranking *Ii*. *Ii* means that *f*(*pi*(*k*)) is the *Ii*th lowest in all *f*(*p*(*k*)). σ*i* is negatively correlated with *k*. μ*i* allows the seeker near the optimal point to search for more details and σ*i* makes seekers search more carefully at the later stage of the search. Figure 2 is an intuitive expression for determining the step length. '→' represents the tendency to lessen step length.

$$\mu\left(l\_i\right) = e^{-\frac{l\_i^2}{2\sigma\_i^2}}\tag{8}$$

**Figure 2.** Determine the search step length.

#### *3.2. Nelder–Mead Simplex Algorithm*

As a derivative-free method, the NM algorithm is used to construct a simplex to solve unconstrained minimization problems. To minimize a mathematical function *f* that has *n* real parameters, the simplex is designed with *n* + 1 vertexes and to replace the worst vertex with a better one each iteration via reflection, expansion, contraction, or reduction.

After the initial simplex is determined, the function *f* is evaluated at each of the vertices, which are later renamed according to its value and then comes the relationship in Equation (9). The basis points and the candidate points that will be used later are listed in Table 1. And they are elected pursuant to the rule shown in Figure 3. The NM algorithm is the first step towards finding an effective substitute for *Xw* along the line segment joining *Xc* and *Xw* as seen in Figure 4. If the first step fails, NM no longer explores that line segment but starts to draw closer to *Xb* by shrinking. Figure 5 shows the evolution of simplex when NM is used to find the optimal point of a 2-dimensional objective function.

$$f(X\_1) \le f(X\_2) \le \cdots \le f(X\_{n+1}) \tag{9}$$

**Table 1.** Basis points and the candidate points of Nelder–Mead (NM).

**Figure 3.** NM evolutionary process.

**Figure 5.** A case of NM algorithm.

Ļėļ:Expansion

### *3.3. New Hybrid SOA: NMSA*

The random factors in the seeker's direction and step length bring SOA good global exploration ability. Nevertheless, it is also the randomness that causes its poor local exploration ability. SOA's later search is always wasted. Besides, when dealing with a multi-extremum problem, SOA may converge to a wrong extreme point because of insufficient exploration to the main extreme point. As a beam search method, the NM simplex method has opposite capability compared with SOA. However, it also produces incorrect outputs due to the quality of the initial simplex. In this paper, a synthetic algorithm of NM and SOA is proposed to make their respective advantages complementary to each other. In every day's searching, the top *n* + 1 seekers of SOA do the NM transformation and iterate once. After that, they begin to share their positions and to make the next day's plan. A diagram of the proposed method NMSA (The hybrid algorithm based on the seeker algorithm and the Nelder–Mead simplex) is given in Figure 6.

**Figure 6.** NMSA algorithm.

To evaluate the performance of the optimization strategy, a set of artificial test functions has been proposed and used as benchmarks. Typically, they can be classified in the term of modality. If there is only one peak in the function landscape, the function is said to be unimodal (e.g., sphere function, booth function, etc.), which is often used for convergence velocity comparison [27]. On the contrary, the function is a multimodal function (e.g., Rastrigin function, Ackley function, Schwefel function, etc.). It turns out to be an extremely important performance test for global optimization because it creates difficulties in finding a global minimum. In such a landscape, the search process can direct the search away from the true optimal solutions and trapped in one of the local minima.

For our proposed offline global optimization algorithm NMSA, its reliability and precision should be tested. Rastrigin function, one of the multimodel functions, is based on De Jong's function with the addition of cosine modulation to produce many local minima. In addition, the small fitness differences in the topology increase the complexity of the problem [28]. Therefore, consider the Rastrigin function as the test function of NMSA. The function is given by the following formula and visualized by Figure 7.

$$f(\mathbf{x}, y) = \mathbf{x}^2 + y^2 - 10\cos(2\pi\mathbf{x}) - 10\cos(2\pi y) + 20\tag{10}$$

**Figure 7.** Visualization of Rastrigin function.

Rastrigin function has several minima. The global one is at (0, 0). This minimum gives a function value *f* = 0.  As shown in Figure 8, based on the respective function value obtained by NMSA and SOA after 10 runs, it can be noted that the NMSA algorithm leads successfully to the global minimum compared with the original SOA. With SOA, nearly half of the optimization processes converge to secondary minima and the remainder need at least 90 iterations to converge to the global minimum. NMSA needs less than 60 iterations to get the correct minimum with no wrong judgments in 10 runs. That is because NM offsets the SOA's weakness in local exploration and improves the directionality in NMSA. In general, NMSA has the following characteristics: (1) Its global and local exploration abilities are well balanced for finding the global minimum; (2) It has combined the ergodicity and directionality effectively for a faster rate of convergence. In Table 2, the search results after 100 iterations are presented in two numerical forms: the coordinates of the searched approximate optimal solution and the corresponding fitness value. The former shows that NMSA has advantages over SOA in finding global optimization. As for the latter, the fitness value is numerically equal to the fitness error because the true optimal fitness value in this example is 0. After comparing this value, we can see that the search accuracy of NMSA is about two orders of magnitude higher than that of SOA.

**Figure 8.** The Fitness curve of the seeker optimization algorithm (SOA) and NMSA; (**a**) is the fitness curve of SOA; (**b**) is the fitness curve of NMSA.

**Table 2.** Results of SOA and NMSA.

#### **4. Modeling and Identification of UC**

#### *4.1. Fractional Order Modeling of UC*

UC, also known as electrochemical double layer capacitor, stores the electrical energy by ion absorption in the double-layer or pseudo-capacitive between electrode and electrolyte. A double layer in the order of 0.3–0.8 nm gradually forms during charging or discharging, which leads to a larger capacitance and higher energy density than ordinary capacitors. In addition, carbon with pores is widely used as the electrode of UC because it provides a larger surface area. Because of that energy storage principle and the use of carbon electrodes, UC performance is affected by distributed surface reactivity, inhomogeneity, and roughness, which can be better characterized by CPE according to the electrochemical impedance spectroscopy (EIS). The impedance curve will deflect downward to a squashed semi-circle if the capacitor parallel with the resistor is replaced by CPE. The CPE can also model the imperfect capacitor behavior of a double-layer capacitor; therefore, this paper proposed a fractional model with 2 CPEs as shown in Figure 9.

**Figure 9.** Ultracapacitor (UC) fractional model.

The impedance of CPEs can be expressed as follows:

$$Z\_{\rm CPE\_1}(j\omega) = \frac{1}{\mathbb{C}\_1(j\omega)^a} \tag{11}$$

$$Z\_{\rm CPE\_2}(j\omega) = \frac{1}{\mathcal{C}\_2(j\omega)^\beta} \tag{12}$$

where *C*1 and *C*2 are the capacitance coefficients, α(0 ≤ α ≤ 1) and β(0 ≤ β ≤ 1) are the fractional orders of CPE. According to the equivalent circuit structure, the transfer function of the fractional model can be delineated by the following:

$$G(s) = \frac{\mathcal{U}(s)}{I(s)} = R\_s + \frac{R\_c}{1 + \mathcal{C}\_1 R\_c s^\alpha} + \frac{1}{\mathcal{C}\_2 s^\beta} \tag{13}$$

where *Rs* represents the series resistance and *Rs* is particularly small for UC. Part 2 which consists of *CPE*1 and *Rc* serves to capture the charge diffusion caused by specifically adsorbed ions. *CPE*2 models the behavior of the double layer. To simplify the solving process, the inverse Laplace transformation in Equation (13) can be divided into three parts due to its linear nature:

$$\mathcal{U}I(t) = \mathcal{U}\_1(t) + \mathcal{U}\_2(t) + \mathcal{U}\_3(t) \tag{14}$$

$$\frac{\mathcal{U}I\_1(t)}{\mathcal{R}\_\delta} = I(t) \tag{15}$$

$$\frac{lL\_2(t)}{R\_c} + C\_1 D\_t^a lL\_2(t) = I(t) \tag{16}$$

$$\,^C\!C\_2D\_t^\beta\!L\_3(t) = I(t) \tag{17}$$

With the method proposed in [29] by Zhao C.N. and Xue D.Y., the solution to the fractional-order differential equation

$$a\_1 D^{\gamma\_1} y(t) + a\_2 D^{\gamma\_2} y(t) + \dots + a\_n D^{\gamma\_n} y(t) = u(t) \tag{18}$$

can be written as

$$y(t) = \frac{u(t) - \sum\_{i=1}^{n} \left(\frac{a\_i}{h^{\gamma\_i}} \sum\_{j=1}^{\left[\left(t-t\_0\right)/h\right]} \alpha\_j^{\left(\gamma\_i\right)} y\_{t-jh}\right)}{\sum\_{i=1}^{n} \frac{a\_i}{h^{\gamma\_i}}}\tag{19}$$

$$
\omega\_0^{(\gamma\_i)} = 1, \omega\_j^{(\gamma\_i)} = \left(1 - \frac{\gamma\_i + 1}{j}\right) \omega\_{j-1'}^{\gamma\_i}, j = 1, \ 2, \dots \tag{20}
$$

For part 1: *U*1(*t*) = *Rs*·*I*(*t*); For part 2: *u*(*t*) = *I*(*t*), *y*(*t*) = *U*2(*t*), *a* = [ *C*1 1 *Rc* ], γ = α 0 ; For part 3: *u*(*t*) = *I*(*t*), *y*(*t*) = *U*3(*t*), *a* = *C*2, γ = β.

### *4.2. Time Domain Response Analysis of UC*

After the current is set to zero, the terminal voltage of UC varies with the time. This phenomenon is called voltage rebound and can be used to calculate the effective resistance of UC with the formula *Rt* = *Ut*/*I* [30]. The response of Nesscap to a single pulse is drawn in Figure 10. When current is suddenly set to zero, *t* = 0 and the resistance at this point *R*0 = *U*0/*I*. As *t* increases, *Rt* increases. If parameters are assigned values: *Rs* = 1.537 mΩ, *Rc* = 5.393 mΩ, *C*1 = 7501 F, *C*1 = 2918 F, α = 0.2699, β = 0.9663, the model response fits well with the real response. Decompose the model response into three parts, as shown in Figure 10, according to the previous section. *U*1 imitates the immediate resistance rebound at *t* = 0 and equals to *U*0, so *Rs* = *R*0. It can be seen both *U*2 and *U*3 contribute to the response for *t* > 0 and *t* < 0. *U*2 mainly represents the obvious inertia properties under the sudden current change while *U*3 mainly represents the capacitance characteristics. The existence of α and β provides a more flexible approach to receive a more accurate model.

**Figure 10.** Single-pulse response of UC.

The resistance of different voltages at 200 A is calculated by using the same method and its variation with voltage is shown in Figure 11. Charging *Rs* keeps steady in the charging process. In discharging experiment, *Rs* increases sharply. In the later parameter identification, *Rs* is determined by the look-up table method.

**Figure 11.** Resistance variation at 200 A.

### *4.3. Parameter Estimation in Time Domain*

A Nesscap UC with the nominal capacitance of 3000 F, steady-state resistance of 0.4 mΩ, power density of 512 W/kg, and energy density of 4.4 Wh/kg is used to carry on charge and discharge experiments at different currents. The rectangular current with an amplitude value of 200 A is the input during charging. When *Rs* can be determined using the look-up table, the remaining five parameters are identified by NMSA. Its objective function is defined as Equation (21), which can well evaluate the ability of the model to describe the real capacitance behavior. *U*(*t*) is the measured voltage and *U*ˆ (*t*) is the derived result from identified model. Sampling time *t* = 0.1 s. The key parameters of NMSA are as follows: population size *S* = 20, number of searching days *g* = 100, dimension of map is 5, reflection factor α = 1, expansion factor β = 2, contraction factor γ = 0.5, shrinkage factor δ = 0.5.

$$fitness = \sqrt{\frac{\sum\_{j=1}^{\left[\left(t\_{end} - t\_0\right)/h\right]} \left(\mathcal{U}(j) - \mathcal{U}(j)\right)}{\left[\left(t\_{end} - t\_0\right)/h\right]}} \tag{21}$$

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

After performing the optimization procedure, the result is *Rc* = 14.31 mΩ, *C*1 = 7021 F, *C*2 = 3799 F, α = 0.7944, β = 0.9994 for charging. The validation test has been done with charging and discharging currents at levels of 300 A and 400 A. The UC voltage is constrained to vary in the voltage range [1.35 V, 2.7 V]. The results are presented in Figure 12. Compared to the fitness results obtained by NMSA without considering the voltage dependence of *Rs*, great improvement was achieved under a discharge condition. The error between Nesscap and the nonlinear model is within 0.02 V for charging and 0.05 V for discharging. The result shows that the model gives an accurate voltage response for 300 A and 400 A current level. Though there is a slight difference between the charging and discharging, which may be caused by changes in other parameters, it can be neglected. The proposed model simulates UC behavior with good accuracy in the time domain.

**Figure 12.** Comparison between model voltage response and UC voltage (The left side is the response curve and the right side it the response error compared with UC voltage) (**a**) Validation test with 300 A charging profile; (**b**) Validation test with 400 A charging profile; (**c**) Validation test with 300 A discharging profile; (**d**) Validation test with 400 A discharging profile.

### **6. Conclusions**

In this paper, a fractional model of UC is applied. The influence of two fractional parameters α(0 ≤ α ≤ 1) and β(0 ≤ β ≤ 1) on the UC model's time domain is analyzed. α represents the inertia properties in UC, β represents the undesirable components in UC. The voltage dependence of *Rs* is discussed and a look-up table is set up to add nonlinear characteristics to the model using the voltage rebound method. The resistance of UC varies greatly during discharging. A new hybrid optimization algorithm NMSA is proposed in this paper and has been proven to have sufficient global search capability and accuracy. After applying 300 A and 400 A level current to the identified model, it can be proved that the current magnitude and status of charging or discharging have almost no influence or slight influence on the UCs behavior in the range [200 A, 400 A].

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

**Funding:** The APC was funded by the Key Science and Technology Development Project: Research on the driving cycle development and the control strategy for PHEV based on the Intelligent Transportation System (Project No.<PHONE_NUMBER>5GX).

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

### **Memory E**ff**ect and Fractional Di**ff**erential Dynamics in Planar Microsupercapacitors Based on Multiwalled Carbon Nanotube Arrays**

### **Evgeny P. Kitsyuk 1, Renat T. Sibatov 2,\* and Vyacheslav V. Svetukhin 1**

Received: 6 December 2019; Accepted: 27 December 2019; Published: 2 January 2020

**Abstract:** The development of portable electronic devices has greatly stimulated the need for miniaturized power sources. Planar supercapacitors are micro-scale electrochemical energy storage devices that can be integrated with other microelectronic devices on a chip. In this paper, we study the behavior of microsupercapacitors with in-plane interdigital electrodes of carbon nanotube array under sinusoidal excitation, step voltage input and sawlike voltage input. Considering the anomalous diffusion of ions in the array and interelectrode space, we propose a fractional-order equivalent circuit model that successfully describes the measured impedance spectra. We demonstrate that the response of the investigated micro-supercapacitors is linear and the system is time-invariant. The numerical inversion of the Laplace transforms for electric current response in an equivalent circuit with a given impedance leads to results consistent with potentiostatic measurements and cyclic voltammograms. The use of electrodes based on an ordered array of nanotubes reduces the role of nonlinear effects in the behavior of a supercapacitor. The effect of the disordering of nanotubes with increasing array height on supercapacitor impedance is considered in the framework of a distributed-order subdiffusion model.

**Keywords:** supercapacitor; anomalous diffusion; fractional derivative; equivalent circuit; distributed-order subdiffusion

### **1. Introduction**

Planar micro-supercapacitors (PMSCs) are micro-scale electrochemical energy storage devices that can be embedded in a chip and integrated with other microelectronic devices. Different materials were suggested for electrodes of PMSC including graphene, nanowires, multiwalled carbon nanotubes (MWCNT), and conducting polymers (see review [1] and references therein). With an in-plane interdigital structure (Figure 1a,b), the PMSCs have a better performance than the conventional sandwich-like analogs of similar size [1]. Due to the spatial separation of the electrodes there is no need for a separator and organic binders. This simplifies the design and simulation of the device [2]. Interdigitized PMSCs possess excellent mechanical and electrical characteristics [3]. Due to design, the decreased ionic diffusion path between electrodes leads to small ion transport resistance, that results in higher power capability [3]. The high electrical conductivity of electrodes and their porosity can mitigate the ion diffusion resistance caused by the electrolyte viscosity [4].

An important problem is the development of reliable models for engineering calculations and diagnostics of PMSC. The simplest way is to extend the effective models of standard supercapacitors taking the geometry and specifics of PMSCs into account. One of the successful directions in phenomenological modeling of conventional supercapacitors is based on the use of fractional calculus [5–12]. Usually, these models are formulated in terms of equivalent circuits containing a constant phase element (CPE) sometimes called a capacitor of fractional order or a fractal element [13–15]. This element has an impedance with a frequency dependence proportional to (*j*ω) −α. Bertrand, Sabatier et al., and Allagui et al. [8,16–18] have shown that fractional order models can effectively describe the dynamical behavior of electrochemical power sources and can be used to build SOC (state of charge) and/or SOH (state of health) estimators. Fractional differential models of supercapacitors were used to describe their operation in pulse power systems [19], to describe the impedance in a frequency range [5,7,20], to study temporal response of supercapacitors [11,21–23], to identify parameters from transient characteristics [24,25], to control supercapacitor systems using fractional state models [7], or fractional nonlinear systems [8,26].

**Figure 1.** Planar supercapacitors studied in this work: photographs of several samples (**a**), planar geometry of PMSC (**b**), SEM images of multiwalled carbon nanotube (MWCNT) array in electrodes for different resolutions (**c**,**d**).

In this paper, for the first time, a fractional differential model for PMSC is proposed. We show that the response of the studied supercapacitors with carbon nanotube-based electrodes is linear and the system is time-invariant. The proposed model well describes the behavior of PMSC that is confirmed by various experimental methods. We assume that the use of electrodes with an array of rod-like nanostructures should reduce the role of nonlinear effects in PMSC behavior. The effect of weak disordering of CNT array for increased height is considered in the framework of a distributed-order subdiffusion model.

### **2. Carbon Nanotubes in Planar Microsupercapacitors**

The studied PMSC samples were fabricated on silicon or borosilicate glass substrates by means of standard microelectronics technologies. In the case of silicon substrates, they were pre-oxidized to form a layer of high-quality dielectric. The oxide thickness was ~0.3 μm. Next, a layer of metal (titanium) was applied to the plates. The measurements were carried out under normal atmospheric conditions. During measurement, the samples were placed into the shielded volume to protect against electromagnetic interference. Contacting the electrodes on the chip was carried out using silver-plated clamps. Multiwalled carbon nanotubes were synthesized by catalytic PECVD process (see details in [27]) (Figure 1). The capacities of a supercapacitor with a solid electrolyte based on polyvinyl alcohol (PVA) were investigated. Phosphoric acid H3PO4 or lithium chloride LiCl were used as the ion-conducting agent.

The model developed here is applied to description of PMSC with electrodes composed of MWCNT forest. Two samples of PMSC are analyzed. They are labeled as D33 and A53, and differ in the width of stripes with MWCNTs: 30 μm for D33 and 50 μm for A53. Otherwise, they are the same, including the area. The measured voltammograms (Figure 2) for fabricated PMSCs demonstrate shapes indicating the absence of electrochemical reactions on electrode-electrolyte interface. These shapes are typical for supercapacitors, where electric double layer plays the main role in charging kinetics. As shown below, these samples are characterized by linear response uniform in time that make these devices very convenient for description and prediction of their behavior under various conditions.

**Figure 2.** Cyclic voltammograms of planar micro-supercapacitor with solid electrolyte PVA + H3PO4 for scanning rates 50 mV/s (**a**) and 25 mV/s (**b**) as well as for the plasma treated sample (**c**) and curves measured at two different temperatures (**d**).

### **3. Linear Response and Anomalous Di**ff**usion Model**

In this section, we consider the linear response in terms of the integral equation describing a renewal generally non-Markovian process. By using the Mellin tranform of its kernel, it can be shown that a time-invariant dynamical process obeys a fractional differential equation of a complex order without additional assumptions [28]. The process is characterized by the order (spectral) distribution, so this representation incorporates both the standard Debye process and the non-Debye process family depending on choice of the spectral distribution.

Generally speaking, the object under investigation in the linear response theory (LRT) is an abstract system with two points: input and output. The input signal *v*(*t*) and output signal *u*(*t*) are linked via some functional relation

$$u(t) = \,\_aF\_t[v(\cdot)].$$

Here, *F* is a response functional and *a* denotes an initial moment of input signal appearance: *v*(*t*) = 0 for *t* < *a*. The lower limit *a* of the functional can be taken as −∞, which allows to include in consideration signals appeared in remote past or composed of several segments separated by silence intervals. For short, we will write *Ft* instead of −∞F*t* .

In the frame of LRT, the response functional *F* should satisfy the following two requirements.

This functional can be represented in the following integral form,

$$u(t) = \int\_{-\infty}^{t} K(t, t')v(t')dt'.\tag{1}$$

• If the kernel *K* is invariant with respect to shift, *K*(*t*, *t* ) = *G*(*t* − *t* ), and *G*(*t*) = 0 for *t* < 0, we deal with the Standard Linear Response Theory (SRLT).

The response formula (1) takes the form of a convolution integral

$$u(t) = \int\_0^t G(t - t')v(t')dt'.\tag{2}$$

An interesting behavior of discharging current in Panasonic supercapacitors (0.22, 0.47 and 1 F) was revealed in [22,23]. This behavior conforms neither of known models based on equivalent circuit because the response is nonuniform in time. Here, we will show that the response of PMSC-MWCNT is linear and the system is time-invariant.

The invariance of response kernel with respect to shift in time opens a wide field for using such power analytical tools as integral transforms. The Mellin transform casts additional light on the sense of fractional operators. Representing function *G*(*t*) via its Mellin transform

$$G(t) = \frac{1}{2\pi i} \int\_{c-i\infty}^{c+i\infty} t^{-\upsilon} \overline{G}(\sigma) d\sigma,$$

and inserting the result into integral term of Equation (2):

$$\int\_0^t G(t - t')\upsilon(t')dt' = \int\_{c - i\infty}^{c + \infty} w(\sigma) \, \_0\mathbf{I}\_t^{1 - \sigma} \upsilon(t) d\sigma.$$

Here,

$$w(\sigma) = \frac{\Gamma(1-\sigma)}{2\pi i}\overline{G}(\sigma)\_{\prime\prime}$$

and

$$\, \_0\mathrm{I}\_t^{1-\sigma}u(t) = \frac{1}{\Gamma(1-\sigma)} \int\_0^t \frac{v(t')dt'}{(t-t')^\sigma}$$

is an integral operator of complex order μ = 1 − σ. Passage from the integral to fractional derivatives can be realized say by means of regularization procedure (derivation of a finite part under μ < 0). After multiplying on the weight factor and integrating over σ, this term yields the fractional derivative of distributed order according to the spectral function *w*(*s*) [29,30]. As a result, Equation (2) takes the form

$$\_0D\_t^{\langle w(\cdot)\rangle}v(t) = u(t).$$

The spectral function *w*(σ) defines the LRT kinetics and links it with dynamics of an open system [28].

The application of fractional equations to the analysis of supercapacitors is usually justified by heterogeneity and complexity of porous electrodes [16,23,31–33]. These equations are often interpreted in terms of anomalous diffusion processes. Anomalous diffusion is usually associated with power law expansion of the diffusion packet, Δ(*t*) ∝ *t* α/2. The case 0 <α< 1 corresponds to subdiffusion, and the case α > 1 to superdiffusion.

Different specific mechanisms can lead to anomalous diffusion processes [33,34]. The features of anomalous diffusion are often described within the continuous time random walk (CTRW) model. The decoupled version of CTRW model implies that different jump lengths and waiting (trapping) times between two successive jumps are independent random variables. Fractional-order diffusion equation can be obtained as the asymptotic equation of CTRW with power law distributions of random trapping times τ: *P*{τ > *t*} ∝ *t* −ν (0 <ν< 1). Considering CTRW with space dependent jump probabilities, Barkai et al. [35] derived a time-fractional Fokker–Planck equation (FFPE) for the case, when the mean waiting time diverges. FFPE describes anomalous diffusion in an external force field, it is derived in the following form [35,36],

$$\frac{\partial \mathbf{c}(\mathbf{x},t)}{\partial t} = \,\_0\mathbf{D}\_t^{1-\nu} \mathbf{K} \frac{\partial}{\partial \mathbf{x}} \left[ \frac{\partial \mathbf{c}(\mathbf{x},t)}{\partial \mathbf{x}} + \frac{\mathbf{c}(\mathbf{x},t)}{kT} \frac{\partial V}{\partial \mathbf{x}} \right],\tag{3}$$

where *c*(*x*, *t*) is particle concentration, *V*(*x*) is an external potential, *K* is a generalized diffusion coefficient, *kT* is the thermal energy.

Here,

$$\,\_0\mathrm{D}\_t^{1-\nu}c(\mathbf{x},t) = \frac{1}{\Gamma(\nu)}\frac{\partial}{\partial t}\int\_0^t \frac{c(\mathbf{x},\tau)}{(t-\tau)^{1-\nu}}d\tau, 0 < \nu \le 1,$$

is the Riemann–Liouville derivative of fractional order 1 − ν [37]. Fundamental solutions to Equation (3) can be found in [34,35,38].

The fractional diffusion equation for the simplest case *V*(*x*) = 0 has the following form [39],

$$\frac{\partial \mathbf{c}(\mathbf{x},t)}{\partial t} = K\_0 \mathbf{D}\_t^{1-\nu} \frac{\partial^2 \mathbf{c}(\mathbf{x},t)}{\partial \mathbf{x}^2}. \tag{4}$$

The Fourier transform of Equation (4) leads to the expression

$$(i\omega)^{\vee}\vec{\varepsilon}(\mathbf{x},\omega) = \mathsf{K}\Delta\vec{\varepsilon}(\mathbf{x},\omega). \tag{5}$$

The flux density for diffusion described by (4) takes the form

$$J(\mathbf{x}, t) = K \, \_0\mathrm{D}\_t^{1-\nu} \frac{\partial \mathbf{c}(\mathbf{x}, t)}{\partial \mathbf{x}}.\tag{6}$$

The subdiffusive generalization of Warburg's impedance for a semi-infinite medium can be written in the form [40]

$$Z(j\omega) = B(i\omega)^{-(1-\nu/2)}.\tag{7}$$

Here, *B* is a frequency-independent parameter.

Based on fractional Equations (3) and (4), different generalizations of diffusion impedances can be obtained for different geometries and boundary conditions (see [32,40–42] and references therein). In our recent paper [43], effects produced by ion subdiffusion in lithium-ion cell are described within the fractional-order generalization of the single-particle model. Here, we use relation of fractional-order impedance to anomalous ion diffusion to interpret the behavior of PMSC.

### **4. Equivalent Circuit Model**

Considering the anomalous diffusion of ions in planar electrodes of PMSC and electrolyte in the interelectrode space, we obtained the following simple equivalent circuit of PMSC (Figure 3). It contains CPEs corresponding to ion subdiffusion in the MWCNT array and, in the interelectrode space, a capacitor *C*DL and a resistor *R*DL related to the electric double layer capacity and electrolyte resistance. The corresponding impedance is

$$Z(s) = R + \frac{1}{\mathbb{C}\_{\beta}s^{\beta}} + \frac{1}{\mathbb{C}\_{a}s^{a} + \left[R\_{\rm DL} + 1/\left(\mathbb{C}\_{\rm DL}s\right)\right]^{-1}}\tag{8}$$

Fitting measured impedance spectra, we determined parameters of equivalent circuit for samples D33 and A53 at room temperature (Table 1). To fit the data, we used the program "EIS Spectrum Analyzer" [44]. An example of such fitting for D33 sample is shown in Figure 4. The value of parameter β can be explained by the weak influence of ion diffusion in the interelectrode space on the impedance spectrum; therefore, the interelectrode region can be described using a series-connected resistor and an almost "ideal" capacitance. Parameter α is related to ion subdiffusion in MWCNT array. According to Equation (7), value α = 0.7 corresponds to subdiffusion parameter ν = 2(1 − α) = 0.6.

**Figure 3.** Simple equivalent circuit model of PMSC.

**Table 1.** Parameters of the equivalent circuit model.

**Figure 4.** Nyquist plot (**a**) and frequency dependence of phase shift (**b**) for PMSC (D33) with PVA + LiCl electrolyte. Circles are experimental data, dots are the result of fitting by impedance Equation (8).

To find electric current as a response of PMSC, we use the numerical inverse Laplace transformation of the following relation: ˜*I*(*s*) = *V*˜(*s*)[*Z*(*s*)]−1 . Note that usually this relation leads to disagreement of numerical results with measured data, because *Z*(*s*) depends on operating voltage and it can not be considered as a constant characteristic of the device for all voltages. Moreover, nonlinear effects can take place, whereas *Z*(*s*) is obtained for the linearized process.

We use the following input voltage signals and their Laplace transforms.

For potentiostatic charging *V*(*t*) and its Laplace transform are

$$V(t) = V\_0 H(t), \quad \mathcal{V}(s) = \frac{V\_0}{s},$$

where *H* is the Heaviside function. For the discharging process, we have

$$V(t) = V\_0[H(t) - H(t - \theta)], \quad \mathcal{V}(s) = V\_0 \frac{1 - e^{-\theta s}}{s},$$

where θ is a charging time. To plot discharging curves in log-log scale, we shift obtained *I*(*t*) in time to charging duration θ.

For cyclic voltammetry, the following voltage signal is used

$$V(t) = 2V\_0 \left(\frac{t}{\Theta} + 2\sum\_{k=1}^{\infty} (-1)^k \left(\frac{t}{\Theta} - k\right) \mathbb{H}\left(\frac{t}{\Theta} - k\right) - \frac{1}{2}\right).$$

Its Laplace transform has the form

$$\mathcal{V}(\mathbf{s}) = 2V\_0 \left( \frac{1}{\Theta \mathbf{s}^2} + \frac{2}{\Theta \mathbf{s}^2} \sum\_{k=1}^{\infty} (-1)^k \mathbf{e}^{-k\theta \mathbf{s}} - \frac{1}{2} \right) = V\_0 \left( \frac{2}{\Theta \mathbf{s}^2} \frac{1 - \mathbf{e}^{-\theta \mathbf{s}}}{1 + \mathbf{e}^{-\theta \mathbf{s}}} - 1 \right) = V\_0 \left( \frac{2}{\Theta \mathbf{s}^2} \text{ch} \left( \frac{\theta \mathbf{s}}{2} \right) - 1 \right).$$

We use the Numerical Laplace Inversion package (Wolfram Mathematica) [45] to invert the Laplace transform ˜*I*(*s*) = *V*˜(*s*)[*Z*(*s*)]−1 . Calculated charging and discharging curves, and voltammograms are shown in Figures 5 and 6 and compared with experimental data. There were no fitting parameters, all parameters have been determined from impedance spectra. Thus, we observe quite satisfactory agreement of model responses with experimental data. This means that the PMSC response in potentiostatic and cyclic voltammetry measurements are linear and uniform in time and parameters of equivalent circuit model do not depend on operating voltages.

**Figure 5.** Comparison of measured charging and discharging (θ = 100 and 200 s) curves (lines) for PMSC D33 (**a**) and A53 (**b**) with linear responses predicted by the equivalent circuit model with parameters indicated in Table 1 (points).

**Figure 6.** Comparison of measured voltammograms (lines) for PMSC D33 (**a**) and A53 (**b**) with the linear responses predicted by the equivalent circuit model with parameters indicated in Table 1 (points).

Effective capacities of the studied PMSC are estimated to be 102 μF for D33 and 54 μF for A53. D33 is characterized by a smaller electrode stripe width (30 μm) in comparison with A53 (50 μm). The difference in capacities of these samples can be explained by the difference in surface area of contact between the MWCNT array and electrolyte. Despite the fact that the MWCNT array is porous and saturated with electrolyte, it is quite dense and the surface layers of the MWCNT array-electrolyte contact play the main role.

The deviation of the tails of theoretical discharge curves from the experimental data can be explained by variations in the parameter of the constant phase element considered in the next section. These variations can be related to disordering of the array with the height of MWCNT nanotubes.

### **5. Distributed-Order Fractional Model**

To determine the dependence of PMSC capacitance on nanotube height, the MWCNT arrays with similar morphology, but of different heights, were obtained by varying the synthesis time. The dependence of array height on the synthesis time is shown in Figure 7 (left panel). Next, the capacitances of PMSC with PVA+LiCl electrolyte were measured. Figure 7 (right panel) shows the dependence of specific capacity on the height of the array. The capacitive efficiency of using the MWCNT array height *h*, or ratio *C*/*h*, decreases with *h*. The use of hydrophilic treatments can increase the PMSC capacity by activating the unused area of MWCNT array, however, the power of the device does not increase. Such behavior can be explained by the small variation of the morphology of MWCNT forest for increased heights. These changes can modify the type of ion diffusion inside electrode at different heights (e.g., decrease the subdiffusion exponent). Based on this assumption, we propose distributed-order fractional circuit model and derive corresponding impedances for different distribution of subdiffusion parameter on height.

**Figure 7.** Dependence of MWCNT height on duration of array synthesis (**a**) and the efficiency of using the MWCNT array height for the device capacitance (**b**). Circles correspond to experimental data, lines are approximation by a cubic polynomial function (**a**) and a linear function (**b**).

Consider two cases of distributed-order fractional models for PMSC. The first model implies subdiffusive behavior of fractional order homogeneously distributed over entire volume of electrodes. The second case assumes changing of subdiffusion exponent (fractional order) with height of the array. We show that these models lead to different PMSC impedances and the second case is more preferable to describe our measurements.

### *5.1. Subdi*ff*usion Exponent Distributed Over Electrodes*

For distributed-order fractional model, the subdiffusion equation can be written in the following form [33],

$$\frac{\partial \mathbf{c}(\mathbf{x},t)}{\partial t} = \mathbf{K} \, \frac{\partial^2}{\partial \mathbf{x}^2} \int\_0^1 d\nu \, \rho(\nu) \, \_0 \mathbf{D}\_t^{1-\nu} \mathbf{c}(\mathbf{x},t),\tag{9}$$

where ρ(ν) is a probability density function for subdiffusion exponent 0 < ν ≤ 1, *K* is a generalized diffusion coefficient.

Equation (9) can be derived from the continuity equation with the following relation between current density and concentration:

$$J(\mathbf{x}, t) = -\frac{\partial}{\partial \mathbf{x}} \iint d\boldsymbol{\nu} \,\rho(\boldsymbol{\nu}) \, \mathcal{K}\_{\boldsymbol{\nu}} \, \_0 \mathcal{D}\_t^{1-\boldsymbol{\nu}} \, \mathcal{E}(\mathbf{x}, t)$$

Using the Fourier transformation of this distributed-order fractional Fick's law,

$$\tilde{f}(\mathbf{x},\omega) = -\frac{\partial}{\partial \mathbf{x}} \int\_0^1 d\boldsymbol{\nu} \,\rho(\boldsymbol{\nu}) \, \mathcal{K}\_{\boldsymbol{\nu}} \left(i\boldsymbol{\omega}\right)^{1-\nu} \tilde{c}(\mathbf{x},\boldsymbol{\omega}),$$

and solution of the diffusion equation in the Fourier domain,

$$\vec{\pi}(\mathbf{x},\omega) = \mathbb{C}\_0 \exp \left[ -\mathbf{x} \left( i\omega \right)^{1.2} \left( \int\_0^1 d\nu \ \rho(\nu) \ K\_{\nu} \left( i\omega \right)^{1-\nu} \right)^{-1.2} \right].$$

we obtain the following impedance,

$$Z(i\omega) \propto \left[ i\omega \sqrt{\int\_0^1 d\nu \,\rho(\nu) \, K\_{\nu}(i\nu)}^{-\nu} \right]^{-1}.$$

For the case of weak dependence of *K*ν on ν and uniform distribution of order between νmin and νmax, the impedance takes the form

$$Z(i\omega) = B\sqrt{\frac{\ln(i\omega)}{\left(i\omega\right)^{2-\nu\_{\rm min}} - \left(i\omega\right)^{2-\nu\_{\rm max}}}},\tag{10}$$

where *B* is a constant parameter.

### *5.2. Variation of Subdi*ff*usion Order with MWCNT Height*

Let us assume the presence of variation in morphology of CNT forest with height of array. Subdiffusion perpendicular to MWCNT inside the electrodes is described by Equation (4), but the order ν depends on height. A separate layer of thickness *dh* can be characterized by impedance Equation (7). These layers are assumed to operate in parallel, and the equivalent circuit presented in Figure 8 should be used instead of the previous model shown in Figure 3.

**Figure 8.** (**a**) Distributed order fractional circuit model of PMSC. (**b**) Nyquist diagrams for four cases of fractional PMSC model.

As a result, we obtain the following impedance of anomalous diffusion in electrodes,

$$Z\_{d \gets PE}(i\omega) = \left[ \bigcap\_{\nu=\text{min}}^{\nu\_{\text{max}}} b\_{\nu}(i\omega)^{1-\nu/2} \rho(\nu) \, d\nu \right]^{-1}. \tag{11}$$

If we assume uniform distribution of ν between νmin and νmax and constant *b*ν, we have

$$Z\_{d \text{CPE}}(i\omega) = \frac{B \text{ } \ln(i\omega)}{\left(i\omega\right)^{1-\nu\_{\text{min}}/2} - \left(i\omega\right)^{1-\nu\_{\text{max}}/2}}.\tag{12}$$

In Figure 8b, Nyquist plots are presented for four cases of subdiffusion index values. Points denoted as νmax and νmin corresponds to impedance of the circuit presented in Figure 8 with single α-CPE, characterized by νmin = 0.5 or νmax = 0.9. Points denoted by V (volume) and H (height) correspond to the similar circuit including impedances described by Equation (10) or Equation (12), respectively. The shapes of these Nyquist plots look similar, but slope strongly depend on the type of distribution of ν.

Figure 9 demonstrates charging and discharging curves and voltammograms for the circuit presented in Figure 8 with νmin = 0.5, νmax = 0.9. Other parameters are *R* = 1275 Ohm, *C*β = 0.101 × 10−3 sβ/Ohm, β = 0.91, *R*DL = 0, *C*DL = 12.5 μF. Discharging curves indicate the presence of the memory effect, i.e., dependence on the charging prehistory. Cyclic voltammograms of the distributed-order model have the form similar to the single order model. The difference can be seen by variation of νmin. Figure 10 shows two IV-curves for νmin = 0.5 and νmin = 0.7. Parameters

are the same as for Figure 9 except νmin. Due to effect of disordering, smaller values of νmin can be related to higher MWCNT array. For increased height of array, the capacity of PMSC grows, but the dependence of *C* on *h* is not direct proportionality. The shape of cyclic voltammogram for the same scanning rate is changed. This shape is closer to rectangle for higher values of νmin. This is qualitative explanation of the behavior observed in Figure 7.

**Figure 9.** Theoretical charging-discharging curves (**a**) and cyclic voltammograms (**b**) predicted by the distributed-order fractional model.

**Figure 10.** Calculated cyclic voltammograms for the distributed-order fractional model of PMSC for two values of νmin.

### **6. Conclusions**

We studied the behavior of PMSC with electrodes of MWCNT array under sinusoidal excitation, step voltage input, and sawlike voltage input. Studied PMSC samples were fabricated on glass substrates by means of standard microelectronics technologies. MWCNTs were synthesized by a catalytic PECVD process. For the first time, a fractional differential model for planar supercapacitor is proposed. The corresponding equivalent circuit model, taking the anomalous diffusion of ions in the array and interelectrode space into account, successfully describes the measured impedance spectra. We demonstrate that the response of the investigated micro-supercapacitors is linear and uniform in time. The numerical inversion of the Laplace transforms for electric current response in an equivalent circuit with a given impedance leads to results consistent with potentiostatic measurements and cyclic voltammograms. Thus, the use of electrodes based on an ordered array of nanotubes reduces the role of nonlinear effects in the behavior of a supercapacitor. In addition, the effect of the disordering of nanotubes with increasing array length on the impedance of a supercapacitor is discussed in the framework of a distributed-order subdiffusion model. The dependence of PMSC capacitance on the height of MWCNT arrays in electrodes is measured. The capacity grows with *h*, while the ratio *C*/*h* decreases. The explanation of such behavior is given within the fractional differential model

assuming modification of the array morphology for the increased height that leads to a variation of the subdiffusion exponent.

**Author Contributions:** Conceptualization, E.P.K., V.V.S. and R.T.S.; methodology, R.T.S.; software, R.T.S.; validation, E.P.K. and R.T.S.; formal analysis, R.T.S. and E.P.K.; resources, V.V.S.; data curation, V.V.S.; writing–original draft preparation, R.T.S.; writing–review and editing, R.T.S. and E.P.K.; supervision, V.V.S.; project administration, V.V.S. and E.P.K.; funding acquisition, E.P.K. and R.T.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported by the Ministry of Science and Higher Education of the Russian Federation through the assignment of 2020 year for Scientific-Manufacturing Complex "Technological Centre". Renat Sibatov acknowledges the financial support from the Russian Science Foundation (project 19-71-10063).

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

### **References**

Electronics, Control, and Instrumentation (Cat. No.97CH36066), New Orleans, LA, USA, 14 November 1997; Volume 3, pp. 1080–1085.

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

*Article*

### **Development and Experimental Validation of a Supercapacitor Frequency Domain Model for Industrial Energy Applications Considering Dynamic Behaviour at High Frequencies**

### **Gustavo Navarro \*, Jorge Nájera, Jorge Torres , Marcos Blanco, Miguel Santos and Marcos Lafoz**

CIEMAT, Government of Spain, 28040 Madrid, Spain<EMAIL_ADDRESS>(J.N<EMAIL_ADDRESS>(J.T<EMAIL_ADDRESS>(M.B<EMAIL_ADDRESS>(M.S<EMAIL_ADDRESS>(M.L.)

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +34-91-335-71-94

Received: 31 December 2019; Accepted: 27 February 2020; Published: 4 March 2020

**Abstract:** Supercapacitors, one of the most promising energy storage technologies for high power density industrial applications, exchange the energy mostly through power electronic converters, operating under high frequency components due to the commutation. The high frequency produces important effects on the performance of the supercapacitors in relation to their capacitance, inductance and internal resistance (ESR). These parameters are fundamental to evaluate the efficiency of this energy storage system. The aim of the paper is to obtain an accurate model of two supercapacitors connected in series (functional and extrapolated unit) to represent the frequency effects for a wide range of frequencies. The methodology is based on the definition of an appropriate equivalent electric circuit with voltage dependance, obtaining their parameters from experimental tests, carried out by means of electrochemical impedance spectroscopy (EIS) and the use of specific software tools such as EC-Lab® and Statgraphics Centurion®. The paper concludes with a model which reproduces with accuracy both the frequency response of the model BCAP3000 supercapacitors, provided by the manufacturer, and the experimental results obtained by the authors.

**Keywords:** supercapacitor; energy storage; modelling; frequency model; capacitance; ESR; experimental validation

### **1. Introduction**

Among the variety of energy storage systems available for industrial applications, supercapacitors (SCs) have emerged as one of the most selected solutions when high power density, high efficiency, long lifespan, and extended range of operating temperatures are required. Specifically, SCs are widely employed in industry fields where there is a need for storing or releasing a relatively low amount of energy in a very short time and large number of charge-discharge cycles [1,2].

Although SCs are used in many applications by themselves for peak-power pulses reduction or compensation of power oscillations, they can also be hybridized with other energy storage systems, thus increasing the performance of the combined system. Transportation, uninterruptible power supply (UPS), and energy harvesting are the most relevant applications. Firstly, SCs are used in automotive industry to increase the efficiency of hybrid electric vehicles, either for an efficient restart of the engine after a stop, either for limiting the battery current peaks and improving the regenerative breaking in batteries-SCs hybrid vehicles, or for enhancing the dynamic response and load levelling in case of fuel cell-SCs hybrid vehicles [3–7]. Secondly, the UPS industry employs hybrid batteries-SCs systems for