Fractional-Order Modeling and Parameter Identification for Ultracapacitors with a New Hybrid SOA Method

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.


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][7][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: • The model structure is connected to UC dynamic characteristics; • The using of fractional order reduces the number of parameters; • The resistance voltage dependence has been taken into account; • The hybrid optimization algorithm improves recognition accuracy.

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: p r is the usual notation for the binomial coefficient, n = [(t − a)/h]. For p = m(m > 0), a D p t = d m t dt m represents the derivative of order m. On the contrary, a D p 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: • Caputo definition: 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.

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. On each of the following days, they need to record the lowest place they have ever been and the corresponding altitude. 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 kth day, ith seeker p i (k) has a tendency to move towards the personal best position p ibest (k) and the best position of population p zbest (k). Meanwhile, he tends to inherit the moving direction d i (k − 1) of the yesterday if f (p i (k)) < f (p i (k − 1)) and to reverse it if not. Another reference direction might also exist, which is represented by d iother (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 ith seeker determines his moving direction d i (k) by generating a random angle. If the angle is as shown in Figure 1, the ith seeker moves positively in horizontal.    The search step length l i (k) is calculated by Gaussian membership function in Equation (8). µ i is a decreasing function of ith seeker's fitness ranking I i . I i means that f (p i (k)) is the I i 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.

138
The search step length ( ) is calculated by Gaussian membership function in Equation (8).

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 X w along the line segment joining X c and X w as seen in Figure 4. If the first step fails, NM no longer explores that line segment but starts to draw closer to X b by shrinking. Figure 5 shows the evolution of simplex when NM is used to find the optimal point of a 2-dimensional objective function.

Basis Candidate
Best vertex substitute for along the line segment joining and as seen in Figure 4. If the first step fails,

155
NM no longer explores that line segment but starts to draw closer to by shrinking. Figure 5 shows 156 the evolution of simplex when NM is used to find the optimal point of a 2-dimensional objective 157 function.  Energies 2019, 12, x FOR PEER REVIEW 6 of 13

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.  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    Rastrigin function has several minima. The global one is at (0, 0). This minimum gives a function value f = 0. The nearest other minima are, respectively, at (0, 1), (0, −1), (1, 0), and (−1, 0), and they all give a value f = 1. 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.

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.

229
The impedance of CPEs can be expressed as follows: With the method proposed in [29] by Zhao C.N. and Xue D.Y., the solution to the fractional- can be written as The impedance of CPEs can be expressed as follows: 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: where R s represents the series resistance and R s is particularly small for UC. Part 2 which consists of CPE 1 and R c 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: 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 γ 1 y(t) + a 2 D γ 2 y(t) + · · · + a n D γ n y(t) = u(t) can be written as For part 1: U 1 (t) = R s ·I(t); For part 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 R t = U t /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, R t increases. If parameters are assigned values: R s = 1.537 mΩ, R c = 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 R s = 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.

257
The resistance of different voltages at 200 A is calculated by using the same method and its 258 variation with voltage is shown in Figure 11. Charging keeps steady in the charging process. In 259 discharging experiment, increases sharply. In the later parameter identification, is 260 determined by the look-up table method. 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 R s keeps steady in the charging process. In discharging experiment, R s increases sharply. In the later parameter identification, R s is determined by the look-up table method.

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

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 R s 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Û(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.

Results and Discussion
After performing the optimization procedure, the result is R c = 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 R s , 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.

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 R s 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].

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