Estimation of Supercapacitor Energy Storage Based on Fractional Differential Equations

In this paper, new results on using only voltage measurements on supercapacitor terminals for estimation of accumulated energy are presented. For this purpose, a study based on application of fractional–order models of supercapacitor charging/discharging circuits is undertaken. Parameter estimates of the models are then used to assess the amount of the energy accumulated in supercapacitor. The obtained results are compared with energy determined experimentally by measuring voltage and current on supercapacitor terminals. All the tests are repeated for various input signal shapes and parameters. Very high consistency between estimated and experimental results fully confirm suitability of the proposed approach and thus applicability of the fractional calculus to modelling of supercapacitor energy storage.


Background
As of today, supercapacitors are the main components of many devices and systems, e.g., backup power and electricity recovery systems as well as automotive applications, hybrid vehicles and many others. The ability to accumulate charge without any chemical reactions makes such elements to have hundreds of times higher number of charge/discharge cycles in comparison to typical batteries [1]. Additionally, high charge/discharge rates make them effective for applications in energy recovery systems used for example in transportation or renewable energy sources [2,3]. In all these applications, the key parameter is the information on the amount of accumulated energy in the supercapacitor [4,5]. Unfortunately, the well known relationship for typical capacitors that allows to determine the information, that is (1/2)CU 2 , cannot be used [6]. The amount of accumulated energy cannot be determined on the basis of the voltage on capacitor terminals only. The main reason for this is the diffusion process associated with the charge redistribution [1,7]. This is why many researchers have been trying to determine a Correspondence: r.kopka@po.opole.pl Opole University of Technology, ul. Prószkowska 76, 45-758 Opole, Poland supercapacitor model that would allow estimating the behavior of a real system. Currently, researchers mainly adopt the combinations of typical electronic elements, e.g., RC quadripole or series and parallel combinations of such elements. However, all of these models assume a relationship between supercapacitor current and voltage on its terminal in form of a typical, integer order differential equation [3][4][5]7].
But it turns out that some completely new possibilities for energy estimation in such systems can be obtained by the application of the fractional calculus [8,9]. The noninteger-order differo-integral calculus was proposed over 300 years ago, but important implementation issues are related with the advent of computers and their use in modeling of discrete-time dynamical systems [10][11][12][13][14]. Application of fractional calculus to the problem of supercapacitor parameter estimation is not a new issue. There are many publications in this field [15][16][17][18][19][20][21][22][23][24][25]. The authors perform the task of estimating parameters in both frequency and time domains [26].
This paper is an extended version of the author's conference presentation [27], in which a fractional-order approach has been briefly introduced to estimate energy accumulated in the supercapacitor.
Accurate estimation of parameters of the supercapacitors is also of utmost importance in assessing their reliability [28][29][30][31]. Permanent degradation processes inside the supercapacitor can change the equivalent series resistance and capacitance. Thus, accurate determination of these parameters, based on the proposed method, also allows to accurately assess the performance of the capacitor.
This paper starts with some preliminaries related to fractional-order integration and differentiation. Next, it presents the parameter estimation method used during tests and proposes new energy calculation method based on the fractional calculus. The Results and Discussion section present the calculated energy for various scenarios and compare them with reference (measured) values. Conclusions and contributions are summarized in Conclusions section.

Methods
The use of porous materials in supercapacitors and specific manner of charge accumulation cause that the traditional approaches based on integer order derivative models are not accurate enough. Many researchers have proposed various solutions in form of combination of typical RC elements with constant or variable values [4,7]. But it turns out that definitely better precision can be obtained using noninteger-order differential calculus for defining the relationships between supercapacitor's current and voltage [17,19]. Additionally, such a solution may result in a very simple model structure, while providing very high accuracy [18].

Fractional Order Differo-Integral Calculus
Fractional order differential calculus has been known for over 300 years. However, only recent several years have brought its popularity in modelling of physical phenomena and processes. It is believed that description of dynamics with a derivative or integral of nonintegerorder can be one of the most effective methods for modelling of real properties of many complex phenomena and industrial processes, especially based on novel materials and technologies [10,12,13,[32][33][34]. Noninteger-order differential or integral calculus is a generalization of classical calculus to order α that belongs to the set of real numbers R. The differo-integral operator of order α ∈ R of function f (t) on the range [ a, t] can be written as follows assuming that the function f (t) is multiple times differentiable and integrable. As for the operator (1), there are many definitions of its realization. Such definitions differ in properties and areas of application. The most popular are the Riemann-Liouville, Caputo and Grünwald-Letnikov (GL) definitions [34]. The latter will be used in this paper in the form where the binomial α j is defined as follows In order to obtain a fractional model at discrete moments of time, the GL definition in a discrete form is simplified as There are several discretization schemes for the GL Eq. (4). The most popular ones include backward differences (Euler), trapezoidal (Tustin), and Al Alaoui operators. Using Euler's method, the fractional derivative at discrete time moments k can be presented as The infinite sum of previous samples must be in real systems limited to a finite value due to the limited memory and limited calculation time. Now, the truncated or finite-length discrete-time approximation of GL is where f (l) = 0 for l < 0 and L is the length of the model (6) [23]. Reducing the number of samples results in decreased calculation accuracy. This is important for systems operating in a continuous time. Some other kinds of solution are algorithms approximating fractional differointegrals with integer order models. An example can be Oustaloup recursive filters [35]. Another effective finitelength model is the FFLD, being a combination of the truncated model (6) and a Laguarre-based difference [24,36,37]. All results of identification as well as energy measurements are obtained based on all samples in the (long) observation window L, i.e., with maximum accuracy. Figure 1 presents the step responses of integration and differentiation obtained based on (6), for k = 0, 1, . . . , L and for various values of integration/differentiation order α.
Assuming different values of order α, one can more accurately model different physical processes, especially diffusion ones. a b Fig. 1 Step responses for integrating (a) and differentiating (b) models with various orders α

Parameter Estimation for Fractional Model
The results of all energy measurements and identification procedures presented in this paper were obtained for a supercapacitor charged from a controlled voltage source. In such a system, the supercapacitor current i C (t) must be limited by a resistor R connected in series with the supercapacitor C (Fig. 2). Estimation of all supercapacitor parameters is performed based on quadripole response u C (t) to voltage step u(t) at its input. Choosing the appropriate value of derivative order α allows to account for a supercapacitor model of the physical phenomena related to diffusion processes associated with the charge redistribution during the charging and discharging processes. The parallel resistor r P additionally enables modelling of the leakage current. Using the fractional differential calculus for modeling supercapacitors, the model structure can be of low complexity. For supercapacitor charged from the voltage source, a model consists of only two elements, i.e. a simple RC quadripole (Fig. 2a). For low capacities, the a b c Fig. 2 Supercapacitor RC models, base model (a), expanded with a series resistance (b), and with additional parallel resistance (c) series resistance r S is of importance (Fig. 2b), while the leakage current I L may be additionally represented by the parallel resistance r P (Fig. 2c). Using the fractional order calculus to model the supercapacitor, the relation between the voltage on capacitor terminals and its current can be expressed as follows where the operator d α /dt α means a differentiation operator of order α and the SI unit of C α is [ F/sec 1−α ]. The basic supercapacitor configuration presented in Fig. 2a can be treated as a first-order inertial system and can be represented by the fractional transfer function where T = RC α . Taking into account the series resistance r S (Fig. 2b), the circuit is treated as phase-delaying correction system with the transfer function (compare [24]) where T 1 = C α (R + r S ) and T 2 = r S C α . Additionally, allowing for the parallel resistor r P representing the leakage current I L (Fig. 2c), the system transfer function can be expressed as where K = R/r P + 1, T 1 = C(Rr s /r P + R + r S ) and T 2 = r S C. In the time domain, Eq. (10) can be presented as The time response of the model defined by (11) was obtained by transforming it into the form presented graphically in Fig. 3, where integration and differentiation operations are of fractional order α. This model was used during the process of estimation of supercapacitor parameters. The tested supercapacitor was identified using the system presented in Fig. 4a. The control procedure of the entire system was developed using the Matlab/Simulink software with xPC Toolbox. The system consisted of a desktop PC (xPC Target) with the installed measurement card NI-DAQ and master computer (xPC Host). The computers were interconnected through the Ethernet network. The supercapacitor was charged and discharged by (voltage-controlled) voltage source (Fig. 4b) of current efficiency up to ± 3 A. The measurement system was operated with the sampling frequency of 100 Hz, while all the measurements and analog control signals were processed with 16-bit resolution [25].
The main method for determining the dynamic properties of a system is based on the analysis of the step response [38]. In relation to the system model, this method allows estimation of its parameters. For this study, the step signal with various voltages (0.5/1.0/1.5/2.0/2.7 V) and constant duration (500 s) have been used (see Fig. 5 and Table 2). On the other hand, one of the typical applications of supercapacitors is the accumulation or delivery of energy into the power systems. In this case, the voltage change rate is rather small. To simulate it, the 400 mVpp and 0.03 rad/s signal with 2 V offset was used (Fig. 6). Additionally in order to examine the influence of the voltage and frequency changes on estimated parameters, various values of the latter were used (see Table 3).
There are several methods for estimation of model parameters. The main aim of the identification procedure in the time domain applied in this work was to estimate the vector of unknown parameters θ =[ α, C α , r S , r P ] of fractional model presented by (11). The least squares method was used to minimize the initial error. An optimization criterion involved minimization of the standard error (k) 2 2 , where where u C (k) is the output voltage measured from the tested system at moment k, whileû C (k) is the output voltage from the considered model for the input signal u(k). The identification problem is now reduced to finding a parameter vector θ ∈ ad that would minimize the square criterion J in such a manner that where ad denotes the set of admissible parameter values and N means the simulation time. There are many optimisation algorithms that can be used to solve the problem (13). The results presented in this paper were obtained by implementing the genetic algorithm in the Matlab environment.

Energy Calculation
A change in the energy stored in the supercapacitor depends on the power supplied to the capacitor per unit of time and can be described as follows By expressing the power supplied to the capacitor as a product of the current and voltage on capacitor terminals, the change in energy at given time t can be expressed as The total energy during the time interval [t 1 , t 2 ] can be obtained by integrating the energy changes over that time Accounting for Eq. (7), the total energy storage can be determined as Assuming t 1 = 0 and E t 1 = 0, the total energy stored in the supercapacitor during the time interval [ 0, t] is Note that for α = 1 Eq. (18) can be reduced to the classical one

Results and Discussion
Initially, the procedure for estimating the parameter vector of the supercapacitor model using the fractional calculus was performed. The estimation was performed based on the system presented in Fig. 2c, generating a voltage step or sinusoidal wave at its input. The model responses were calculated based on (11). The results obtained by the a b Fig. 5 Step responses for tested supercapacitor and its fractional model (a) and the model response error (b) two identification procedures are very similar, especially in the case of fractional capacitance C α and the fractional order α (see Table 1). Some differences in the estimates of series resistance r S may be a result of its dependence on frequency. The step signal consists of many high frequency harmonics while the sinusoidal wave only onethe 0.03 rad/s. The presented results were obtained for the commercial supercapacitor Samwha Green-Cap EDLC(DB), rated as 2.7 V with 100 F nominal capacitance and 8 m maximum equivalent series resistance (r S ) at 1 kHz. Figures 5a and 6a show the measured supercapacitor voltage and the calculated model responses, for step and sinusoidal signals, respectively, while Figs. 5b and 6b show the model response error.
All obtained results show high consistency between model responses and real measurements despite the fact that relatively simple models were proposed. Some discrepancies may result from the fact that model parameters should be estimated in the system of supercapacitor charged and discharged using the current source [25]. Also, very high estimates of r P may suggest that this resistance could be excluded from the supercapacitor model shown in Fig. 2c. Those very high estimates and their high discrepancies for different inputs indicate that the test signals used to estimate this parameter are not proper. The model (10) was used as the most general form. However, in order to accurately determine all its parameters, it was necessary to use other procedures and test signals. The value of r P characterizes the leakage current I L and should be determined using the constant voltage signal, but for a very long time -of order of several dozen hours.  Although the main goal of the study was to measure energy, various excitation conditions largely affected all the parameter estimates (see Table 2). For instance, the increase of the voltage step amplitude significantly changed the fractional integration order, as the result of increasing effect of the diffusion phenomena inside the supercapacitor. It can also be seen from Table 2 that the supercapacitor is quite nonlinear. As a result of the integration order changes, the variation of the fractional capacity is also observed. This also applies to sinusoidal excitation. The values of estimated parameters-especially α and C α -depend on the amplitude and frequency (see Table 3). For low frequencies, the amplitude value is important, while for higher frequencies the supercapacitor behaves like being excited with a constant voltage.

Energy Calculation
Figures 7a and 8a show measured values of the voltage and current of the supercapacitor for the configuration as presented in Fig. 4b. These values were used for calculation of the total energy stored in the capacitor (marked as E 1 in Figs. 7b and 8b) according to (16). Just as for parameter identification processes, the calculations were conducted both for voltage step and sinusoidal wave at the system input. The energy calculated in such a manner for each time t was compared with energy calculated based on the voltage and capacity in accordance with (19) (marked as E 3 in Figs. 7b and 8b) and energy calculated with fractional-order calculus (marked as E 2 in Figs. 7b and 8b) according to (18). For Eq. (19), a nominal value of supercapacitor was adopted (C n ), while in (18) the value obtained from the estimation process presented in Table 1 was used. Figure 7b shows results of measurements and energy calculations for voltage step, while Fig. 8b shows this same quantities for sinusoidal wave. Similar calculations were made for different voltage steps and sinusoidal excitations. Figure 9a and calculated energies for two voltage steps of 0.5 V and 2.7 V, respectively. Figure 10 shows the energy changes for a sinusoidal signal with the frequency of 0.03 rad/sec and different amplitudes of 0.1/0.25/0.5 and 0.7 V. It can be seen that the differences in the determined energy values correspond to differences in the estimated values of the fractional order α. The greater the difference from the value − 1, the greater is the difference in the calculated energies.

Discussion
The use of a porous material electrodes in supercapacitors in form of active carbon isolated by a very thin separator and use of charge accumulation mechanisms as so- order α of derivative/integral, one can precisely model phenomena and processes occurring inside the supercapacitor using simple mathematical models. Taking into account the real value of the accumulated energy determined by (16), the integer-order model with nominal parameters (19) under-estimates the amount of energy, while the fractional model (18) indicates almost the same value.
The performed tests and measurements were related to charging and discharging of the supercapacitor by a voltage source. Under industrial conditions, supercapacitors are usually charged and discharged by current sources. This can change the nature of the system because the capacitor is no longer an inertial system but becomes a typical integrating one. However, the measurements conducted by author also indicate the occurrence of diffusion processes in such cases. Anyway, usefulness of the Gründwald-Letnikov derivative/integral is confirmed here. Another issue is related to the implementation of the GL differo-integral operator as, e.g., the finite or truncated GL difference (6), which may be computationally burdensome. In future research, we will compare the Oustaloup [35] and FFLD [24,36,37] approximators to effectively solve the implementation issue.
The amount of energy storage in supercapacitor calculated only on the measured value of supercapacitor terminal voltage and using model (19) is not appropriate. The model (19) is only valid if the capacitor current is characterized by the integer order derivative of the capacitor voltage (i C (t) = du C (t)/dt). This is not true for supercapacitor as a consequence of its construction and used special materials. However, the same problem occurs with very large supercapacitors charged by current source. There are also quite new element as super-batteries. In all these applications, the current changes are not characterized by the integer-order derivative of the terminal voltage as a consequence of the specific properties of these elements.

Conclusions
In this paper, a new approach to estimation of an amount of energy accumulated in supercapacitors has been presented. The analysis has been conducted taking advantage of certain unique properties of fractional-order models. It has been shown that application of such sophisticated modeling leads to very accurate results, which can be obtained even though the models themselves are not of high complexity. This is due to a natural ability of noninteger order dynamics to model diffusion processes, just like charge redistribution in supercapacitors. The results of this paper have confirmed the fractional nature of supercapacitors.