A General Parameter Identiﬁcation Procedure Used for the Comparative Study of Supercapacitors Models

: Supercapacitors with characteristics such as high power density, long cycling life, fast charge, and discharge response are used in di ﬀ erent applications like hybrid and electric vehicles, grid integration of renewable energies, or medical equipment. The parametric identiﬁcation and the supercapacitor model selection are two complex processes, which have a critical impact on the system design process. This paper shows a comparison of the six commonly used supercapacitor models, as well as a general and straightforward identiﬁcation parameter procedure based on Simulink or Simscape and the Optimization Toolbox of Matlab ® . The proposed procedure allows for estimating the di ﬀ erent parameters of every model using a di ﬀ erent identiﬁcation current proﬁle. Once the parameters have been obtained, the performance of each supercapacitor model is evaluated through two current proﬁles applied to hybrid electric vehicles, the urban driving cycle (ECE-15 or UDC) and the hybrid pulse power characterization (HPPC). The experimental results show that the model accuracy depends on the identiﬁcation proﬁle, as well as the robustness of each supercapacitor model. Finally, some model and identiﬁcation current proﬁle recommendations are detailed.


Introduction
Energy storage systems are essential in the industrial, medical, renewable or transportation sectors, as well as other sectors. Some characteristics like high power density, reliability and safety are critical in those sectors, this is why the electrochemical double layer capacitor or the supercapacitor play an important role [1].
Many application areas in which supercapacitors are used can be mentioned like magnetic resonance imaging (MRI) that needs very short pulses with high current [2] or fuel cell supercapacitor hybrid bus, where the supercapacitor satisfy the dynamic power demand [3]. In addition, the supercapacitor can be used for the integration of a photovoltaic power plant [4], grid integration of renewable energies [5] and the improvement of energy utilization for mine hoist applications [6]. However, many applications are limited by the self-discharge behavior in wireless sensor network applications [7], where the new techniques of chemical modification to suppress this phenomenon are shown in reference [8] and reference [9].
In general, the supercapacitors models classify into three categories: electrochemical, mathematical, and electrical. Electrochemical models consist of a set of partial differential-algebraic equations with many parameters. The estimation of the electrochemical model is very accurate [10]. However, the simulation of these models consumes many resources. Mathematical models are an alternative based on three dimensional ordered structures [11]. It can get a good fitting with

Supercapacitors Models
In this section, six representative supercapacitor models are selected from the literature, which cover most of the typical applications. All of them are nonlinear models since this kind of models obtains better accuracy. The selected models are the Stern-Tafel Model [18], Zubieta Model [19], Series Model [20], Parallel Model [21], Transmission Line Model [22] and Thevenin Model [23]. In this section, the electrical equivalent circuit and the parameters of each model are reviewed.

Stern-Tafel Model
The supercapacitor proposed in reference [24] and reference [25] uses the Stern-Tafel model to describe the nonlinear capacitance. This electrochemical model reproduces the double layer capacitance (C T ) related to the nonlinear diffusion dynamics. To do this, the supercapacitor model combines both the Helmholtz's capacitance (C H ) and Gouy-Chapman's capacitance (C GC ) [26], Being Energies 2019, 12,1776 3 of 20 where N p is the number of parallel supercapacitor cells, N s is the number of series of supercapacitor cells, N e is the number of layers of electrodes, d the molecular radius (m), c the molar concentration (mol.m −3 ), A i is the interfacial area between electrode and electrolyte (m 2 ), T is the operating temperature (K), F c is the Faraday constant (C/mol), R is the ideal gas constant (J/(K·mol)), ε is the relative permittivity of the electrolyte material (F/m), and ε 0 is the free space permittivity (F/m) [18]. The model equivalent circuit has a controlled voltage source and an internal resistance, as shown in Figure 1a. This model depends on several parameters where C n is the nominal capacitance (F), V max is the maximum supercapacitor voltage (V), R dc is the internal resistance (Ω), V T is the total voltage (V), and i sd is the self-discharge current (A) which is determined by the Tafel Equation (4) described in reference [27] as: where I f is the leakage current (A), V init is the initial voltage (V), α is the charge transfer coefficient and ∆V is the over-potential (V). The capacitance of the electrochemical model requires only a few data from manufacturer datasheet and universal constant as described in reference [28]. The Simulink implementation is shown in Figure 1b.
where Np is the number of parallel supercapacitor cells, Ns is the number of series of supercapacitor cells, Ne is the number of layers of electrodes, d the molecular radius (m), c the molar concentration (mol.m −3 ), Ai is the interfacial area between electrode and electrolyte (m 2 ), T is the operating temperature (K), Fc is the Faraday constant (C/mol), R is the ideal gas constant (J/(K·mol)), ε is the relative permittivity of the electrolyte material (F/m), and ε0 is the free space permittivity (F/m) [18].
The model equivalent circuit has a controlled voltage source and an internal resistance, as shown in Figure 1a. This model depends on several parameters where Cn is the nominal capacitance (F), Vmax is the maximum supercapacitor voltage (V), Rdc is the internal resistance (Ω), VT is the total voltage (V), and isd is the self-discharge current (A) which is determined by the Tafel Equation (4) described in reference [27] as: where If is the leakage current (A), Vinit is the initial voltage (V), α is the charge transfer coefficient and ∆V is the over-potential (V). The capacitance of the electrochemical model requires only a few data from manufacturer datasheet and universal constant as described in reference [28]. The Simulink implementation is shown in Figure 1b.

Zubieta Model
The proposed model in reference [19] includes a circuit with three parallel RC time constant, Figure 2a. The first branch, with the elements R 0 C 0 , and the voltage-dependent k c ·v c defines the response in seconds. The second branch R 1 C 1 provides the response in the range of minutes. The branch R 2 C 2 represents the response for a time longer than minutes. Finally, a resistor R lk reproduces the leakage resistance.
A simplified equivalent circuit with two branches is shown in reference [29], with a simplified parameter identification procedure through the differential equation of the circuit. Similar studies are proposed in reference [30] in which the model parameters are easily obtained when the supercapacitor is discharged with constant power. In addition, reference [31] proposes a multivariable minimization function to find the parameters, they are validated with a current profile of a hybrid electric vehicle.
The total capacitance and current of the voltage-controlled capacitance implemented in Simscape are shown in Figure 2b, which are defined by (5) and (6): where C 0 is the initial linear capacitance which represents the electrostatic capacitance of the capacitor, and k c a positive coefficient which represents the effects of the diffused layer of the supercapacitor.

Zubieta Model
The proposed model in reference [19] includes a circuit with three parallel RC time constant, Figure 2a. The first branch, with the elements R0C0, and the voltage-dependent kc·vc defines the response in seconds. The second branch R1C1 provides the response in the range of minutes. The branch R2C2 represents the response for a time longer than minutes. Finally, a resistor Rlk reproduces the leakage resistance.
A simplified equivalent circuit with two branches is shown in reference [29], with a simplified parameter identification procedure through the differential equation of the circuit. Similar studies are proposed in reference [30] in which the model parameters are easily obtained when the supercapacitor is discharged with constant power. In addition, reference [31] proposes a multivariable minimization function to find the parameters, they are validated with a current profile of a hybrid electric vehicle.
The total capacitance and current of the voltage-controlled capacitance implemented in Simscape are shown in Figure 2b, which are defined by (5) and (6): where C0 is the initial linear capacitance which represents the electrostatic capacitance of the capacitor, and kc a positive coefficient which represents the effects of the diffused layer of the supercapacitor.

Series Model
The series model is an equivalent circuit obtained through the AC impedance approach, which consists of two parallel RC circuit compound by R 1 (vsc), C 1 (vsc), R 2 (vsc), C 2 (vsc), connected in series with another RC circuit compound by R s and C s (vs), as described in references [20,32,33]. This equivalent circuit shows in the first branch of Figure 3a. In reference [34] a modified version of this circuit was presented, which includes the model proposed by Buller and Zubieta, in order to represent a complete model for a full frequency range. This complete model includes three branches in a parallel compound by R 3 and C 3 , R 4 and C 4 , and the leakage resistance R lk , as shown in Figure 3a. Figure 3b shows the Simscape implementation of the modified series model.

Series Model
The series model is an equivalent circuit obtained through the AC impedance approach, which consists of two parallel RC circuit compound by R1(vsc), C1(vsc), R2(vsc), C2(vsc), connected in series with another RC circuit compound by Rs and Cs(vs), as described in references [20], [32] and [33]. This equivalent circuit shows in the first branch of Figure 3a. In reference [34] a modified version of this circuit was presented, which includes the model proposed by Buller and Zubieta, in order to represent a complete model for a full frequency range. This complete model includes three branches in a parallel compound by R3 and C3, R4 and C4, and the leakage resistance Rlk, as shown in Figure 3a. Figure 3b shows the Simscape implementation of the modified series model.

Parallel Model
The basic parallel model with constant values is described in reference [35] and reference [36]. Reference [37] describes an approximation to calculate the parameters without data acquisition, only using the information provided by a supercapacitor datasheet, as well as the main basic equations to obtain the constant parameters using this information. A modified four parallel RC networks with voltage-dependent parameters are presented in reference [21], and it is shown in Figure 4a. This model is more complex, but it achieves better accuracy. Figure 4b shows the implementations of the modified parallel model in Simscape.

Parallel Model
The basic parallel model with constant values is described in reference [35] and reference [36]. Reference [37] describes an approximation to calculate the parameters without data acquisition, only using the information provided by a supercapacitor datasheet, as well as the main basic equations to obtain the constant parameters using this information. A modified four parallel RC networks with voltage-dependent parameters are presented in reference [21], and it is shown in Figure 4a. This model is more complex, but it achieves better accuracy. Figure 4b shows the implementations of the modified parallel model in Simscape.

Transmission Line Model
Transmission line model is composed of nRC branches in order to reproduce the supercapacitor frequency response from 10 mHz to 1 kHz. This model was proposed for hybrid and electric vehicles, and it was described in reference [38] and reference [39]. This model consists of four parallel networks based on R1, C1(v1), R2, C2(v2), R3, C3(v3) and R4, C4(v4), and a parallel leakage resistance Rlk, as shown in Figure 5a. Reference [22] describes a procedure to estimate the parameters through time response and the equations of the circuit. Also, this model is used to evaluate the supercapacitor physical aging process [40], by estimating the uncertainties of the parameters. Reference [41] uses a different number of networks according to the simulation time step. Figure 5b shows the model implemented in Simscape with the described Equations (5) and (6). (a)

Transmission Line Model
Transmission line model is composed of nRC branches in order to reproduce the supercapacitor frequency response from 10 mHz to 1 kHz. This model was proposed for hybrid and electric vehicles, and it was described in reference [38] and reference [39]. This model consists of four parallel networks based on R 1 , C 1 (v 1 ), R 2 , C 2 (v 2 ), R 3 , C 3 (v 3 ) and R 4 , C 4 (v 4 ), and a parallel leakage resistance R lk , as shown in Figure 5a. Reference [22] describes a procedure to estimate the parameters through time response and the equations of the circuit. Also, this model is used to evaluate the supercapacitor physical aging process [40], by estimating the uncertainties of the parameters. Reference [41] uses a different number of networks according to the simulation time step. Figure 5b shows the model implemented in Simscape with the described Equations (5) and (6).

Transmission Line Model
Transmission line model is composed of nRC branches in order to reproduce the supercapacitor frequency response from 10 mHz to 1 kHz. This model was proposed for hybrid and electric vehicles, and it was described in reference [38] and reference [39]. This model consists of four parallel networks based on R1, C1(v1), R2, C2(v2), R3, C3(v3) and R4, C4(v4), and a parallel leakage resistance Rlk, as shown in Figure 5a. Reference [22] describes a procedure to estimate the parameters through time response and the equations of the circuit. Also, this model is used to evaluate the supercapacitor physical aging process [40], by estimating the uncertainties of the parameters. Reference [41] uses a different number of networks according to the simulation time step. Figure 5b shows the model implemented in Simscape with the described Equations (5) and (6). (a)

Thevenin Model
The equivalent electric circuit of the Thevenin model, which includes several parallel RC and a nonlinear state-of-charge (SOC) voltage-dependent source is described in reference [42]. The SOC is calculated by coulomb counting using (7): with Qinit being the initial supercapacitor charge, QT being the total supercapacitor charge and i(τ) as the supercapacitor current.
In this paper, three RC branches are used to get a better accuracy, where OCV represents the open circuit voltage, R0 represents the internal resistance, and three parallel networks based on R1, C1, R2, C2, R3, and C3 reproduce the supercapacitor dynamic, as shown in Figure 6a. All parameters are state-of-charge dependent. The proposed model applied to a hybrid storage system for an electric vehicle gives a better agreement for a simulated vs. experimental response when 3-branches are used in the model [23]. Figure 6b shows the Simscape implementation.

Thevenin Model
The equivalent electric circuit of the Thevenin model, which includes several parallel RC and a nonlinear state-of-charge (SOC) voltage-dependent source is described in reference [42]. The SOC is calculated by coulomb counting using (7): with Q init being the initial supercapacitor charge, Q T being the total supercapacitor charge and i(τ) as the supercapacitor current.
In this paper, three RC branches are used to get a better accuracy, where OCV represents the open circuit voltage, R 0 represents the internal resistance, and three parallel networks based on R 1 , C 1 , R 2 , C 2 , R 3 , and C 3 reproduce the supercapacitor dynamic, as shown in Figure 6a. All parameters are state-of-charge dependent. The proposed model applied to a hybrid storage system for an electric vehicle gives a better agreement for a simulated vs. experimental response when 3-branches are used in the model [23]. Figure 6b shows the Simscape implementation.

Thevenin Model
The equivalent electric circuit of the Thevenin model, which includes several parallel RC and a nonlinear state-of-charge (SOC) voltage-dependent source is described in reference [42]. The SOC is calculated by coulomb counting using (7): with Qinit being the initial supercapacitor charge, QT being the total supercapacitor charge and i(τ) as the supercapacitor current.
In this paper, three RC branches are used to get a better accuracy, where OCV represents the open circuit voltage, R0 represents the internal resistance, and three parallel networks based on R1, C1, R2, C2, R3, and C3 reproduce the supercapacitor dynamic, as shown in Figure 6a. All parameters are state-of-charge dependent. The proposed model applied to a hybrid storage system for an electric vehicle gives a better agreement for a simulated vs. experimental response when 3-branches are used in the model [23]. Figure 6b shows the Simscape implementation.

Parameters Estimation Procedure
Parametric models explicitly contain differential equations, transfer function or block diagrams. The parameters update could be offline or online. For obtaining the parameters, in the offline mode, the data are stored to later process, on the other hand, in the online mode, the procedure is executed in parallel to the experiment [43]. In the literature, there are many proposed procedures to obtain the model parameters such as e.g., the unscented Kalman filter [44] or the Luenberger-style technique [17].
Taking into account the literature, this paper focuses on the proposal of a practical, interactive, simple and enough general offline procedures for estimating the model parameters. Figure 7a shows the proposed identification procedure block diagram. This procedure can be divided into several steps, shown and described in Table 1.

Parameters Estimation Procedure
Parametric models explicitly contain differential equations, transfer function or block diagrams. The parameters update could be offline or online. For obtaining the parameters, in the offline mode, the data are stored to later process, on the other hand, in the online mode, the procedure is executed in parallel to the experiment [43]. In the literature, there are many proposed procedures to obtain the model parameters such as e.g., the unscented Kalman filter [44] or the Luenberger-style technique [17].
Taking into account the literature, this paper focuses on the proposal of a practical, interactive, simple and enough general offline procedures for estimating the model parameters. Figure 7a shows the proposed identification procedure block diagram. This procedure can be divided into several steps, shown and described in Table 1.

Steps Description
1 Apply the identification current profile to obtain supercapacitor current and voltage waveforms (identification data) from the experimental test. E.g., as shown in Section 4.2: current profiles and supercapacitor voltage response (a), (b) and (c).
2 Select and build the equivalent circuit model in Simulink or Simscape through a block diagram or circuit. E.g., as shown in Section 3: Figure 7b.
3 Create a new experiment in Simulink and to import the identification data. Simulate the model with the initial parameters and the identification current profile to obtain the simulation data. E.g., as shown in Section 3: Figure 7c. 4 Choose the variables and their limits to estimate their value. E.g., as shown in Section 3: Figure  7c.

Steps Description
1 Apply the identification current profile to obtain supercapacitor current and voltage waveforms (identification data) from the experimental test. E.g., as shown in Section 4.2: current profiles and supercapacitor voltage response (a), (b) and (c).
2 Select and build the equivalent circuit model in Simulink or Simscape through a block diagram or circuit. E.g., as shown in Section 3: Figure 7b.
3 Create a new experiment in Simulink and to import the identification data. Simulate the model with the initial parameters and the identification current profile to obtain the simulation data. E.g., as shown in Section 3: Figure 7c. 4 Choose the variables and their limits to estimate their value. E.g., as shown in Section 3: Figure 7c.

5
Set up optimization options (optimization method, algorithm, and parameter and function tolerance). E.g., as shown in Section 3: Figure 8.

6
Run the parameter estimation process applying the selected optimization solver (E.g., sum-squared error) to match the identification data with the simulation data. E.g., as shown in Section 3: Figure 7c. If the error is not small enough, return to step 1 ( 1 ); or change the identification method and return to step 3 ( 2 ); or modify the current profile and return to step 2 ( 3 ), Figure 7a. 7 Once the model parameters have been obtained from the identification data, the next step is to verify the model response using the application current profile and the application data. For that, it is necessary to compare the application data with the new simulated data, using the obtained parameters in step 6, E.g., as shown in Section 4.2: Figure 9d,e. If the error is not small enough, return to step 1 ( 1 ); or change the identification method and return to step 3 ( 2 ); or modify the current profile and return to step 2 ( 3 ), Figure 7a.
In step 5, the optimization method has to be selected. This paper uses an offline parameters estimation based on the error minimization between the measured and simulated supercapacitor voltage. The iterative procedure tunes the supercapacitors model parameters (p) to get a simulated response (V s ) that tracks the measured response (V m ), with a finite number of samples (n). To do that, the solver minimizes the next cost function for each current profile: where p varies between zero and infinity (e.g., 0 to 10 10 ).
The minimization problem is carried out with Simulink ® Design Optimization™ of Matlab (Version R2018b, MathWorks, Natick, MA, USA). This toolbox provides an interactive interface that helps to minimize the square of the error between the measured and simulated supercapacitor voltage, using the nonlinear least squares method for parameters estimation. This method is selected in the user interface as shown in Figure 8. Set up optimization options (optimization method, algorithm, and parameter and function tolerance). E.g., as shown in Section 3: Figure 8.

6
Run the parameter estimation process applying the selected optimization solver (E.g., sumsquared error) to match the identification data with the simulation data. E.g., as shown in Section 3: Figure 7c. If the error is not small enough, return to step 1 (); or change the identification method and return to step 3 (); or modify the current profile and return to step 2 (), Figure 7a. 7 Once the model parameters have been obtained from the identification data, the next step is to verify the model response using the application current profile and the application data. For that, it is necessary to compare the application data with the new simulated data, using the obtained parameters in step 6, E.g., as shown in Section 4.2: Figure 9d and Figure 9e. If the error is not small enough, return to step 1 (); or change the identification method and return to step 3 (); or modify the current profile and return to step 2 (), Figure 7a.
In step 5, the optimization method has to be selected. This paper uses an offline parameters estimation based on the error minimization between the measured and simulated supercapacitor voltage. The iterative procedure tunes the supercapacitors model parameters (p) to get a simulated response (Vs) that tracks the measured response (Vm), with a finite number of samples (n). To do that, the solver minimizes the next cost function for each current profile: where p varies between zero and infinity (e.g., 0 to 10 10 ) . The minimization problem is carried out with Simulink® Design Optimization™ of Matlab (Version R2018b, MathWorks, Natick, MA, USA). This toolbox provides an interactive interface that helps to minimize the square of the error between the measured and simulated supercapacitor voltage, using the nonlinear least squares method for parameters estimation. This method is selected in the user interface as shown in Figure 8. This method uses the Simulink function named as lsqnonlin, that requires at least (2k + 1) simulations per iteration, where k is the number of parameters to be estimated [46]. The required CPU time and memory increase as a function of the numbers of parameters and their initial values. The offline runtime estimation is in the order of minutes.
If runtime estimation has to be reduced, other techniques based on the layered technique to break the global optimization into a smaller task [47], or based on differential mutation strategy [48], or based on genetic programming [49], among others, could be used, although the flexibility and simplicity provided by the Simulink user interface could be affected. This method uses the Simulink function named as lsqnonlin, that requires at least (2k + 1) simulations per iteration, where k is the number of parameters to be estimated [46]. The required CPU time and memory increase as a function of the numbers of parameters and their initial values. The offline runtime estimation is in the order of minutes.
If runtime estimation has to be reduced, other techniques based on the layered technique to break the global optimization into a smaller task [47], or based on differential mutation strategy [48], or based on genetic programming [49], among others, could be used, although the flexibility and simplicity provided by the Simulink user interface could be affected.
On the other hand, the algorithm selected is the Trust-Region-Reflective, which is based on a gradient process with a trial step by solving a trust region. Specific details of the algorithm can be found in reference [45]. Additional information is detailed in reference [50], in which the process of how to import, analyze, prepare and estimate model parameters in Simulink is described.
Using the proposed procedure, based on Simulink ® Design Optimization™ of Matlab, the most model can be built, from a practical point of view. Nevertheless, this procedure is limited by the optimization methods and algorithms included in Simulink.

Supercapacitor Testing System
The experimental setup includes a supercapacitor, a data acquisition system, a power source, and an electronic load, as shown in Figure 9. The supercapacitor used to develop the test has been the Maxwell BCAP3000. An equivalent bidirectional current source compound of the electronic load and the power source, connected in parallel, emulates the current profile. This equivalent current source includes the typical regenerative breaking present in automotive applications. The experimental current profile and the data acquisition system are conducted using the following set of equipment listed in Table 2: On the other hand, the algorithm selected is the Trust-Region-Reflective, which is based on a gradient process with a trial step by solving a trust region. Specific details of the algorithm can be found in reference [45]. Additional information is detailed in reference [50], in which the process of how to import, analyze, prepare and estimate model parameters in Simulink is described.
Using the proposed procedure, based on Simulink® Design Optimization™ of Matlab, the most model can be built, from a practical point of view. Nevertheless, this procedure is limited by the optimization methods and algorithms included in Simulink.

Supercapacitor Testing System
The experimental setup includes a supercapacitor, a data acquisition system, a power source, and an electronic load, as shown in Figure 9. The supercapacitor used to develop the test has been the Maxwell BCAP3000. An equivalent bidirectional current source compound of the electronic load and the power source, connected in parallel, emulates the current profile. This equivalent current source includes the typical regenerative breaking present in automotive applications. The experimental current profile and the data acquisition system are conducted using the following set of equipment listed in Table 2:  All these elements have been synchronized with a computer running to manage the data logging and supervisory control using LabVIEW® software.

Supercapacitor Test Schedule
The parameter identification procedure uses three different current profiles. The current profile i1 is a current step, Figure 10a; the current profile i2 are repetitive charging current steps applied until to reach the maximum supercapacitor voltage, Figure 10b; and the current profile i3 is a dynamic charge-discharge current step modulated in amplitude and time applied until the middle value of the supercapacitor voltage range, [51], Figure 10c. From the modeling perspective, the validation current profile must be more dynamic in amplitude and frequency than the identification current profile, as shown in Figure 10d and Figure 10e.  All these elements have been synchronized with a computer running to manage the data logging and supervisory control using LabVIEW ® software.

Supercapacitor Test Schedule
The parameter identification procedure uses three different current profiles. The current profile i 1 is a current step, Figure 10a; the current profile i 2 are repetitive charging current steps applied until to reach the maximum supercapacitor voltage, Figure 10b; and the current profile i 3 is a dynamic charge-discharge current step modulated in amplitude and time applied until the middle value of the supercapacitor voltage range, [51], Figure 10c. From the modeling perspective, the validation current profile must be more dynamic in amplitude and frequency than the identification current profile, as shown in Figure 10d,e. These identification current profiles apply to those models aforementioned in section 2 to obtain their parameters. The current profile applied in every model is shown in Table 3. Table 3. Identification Current Profiles Used to Supercapacitor Parameters Estimation.

i1 i2 i3
Stern-Tafel ✓ -- These identification current profiles apply to those models aforementioned in Section 2 to obtain their parameters. The current profile applied in every model is shown in Table 3.
The robustness and accuracy of the supercapacitor models are evaluated by means of different standardized test profiles, which include the Hybrid Pulse Power Characterization (HPPC) test and European Urban Driving Cycle (ECE15) for long-time responses. Figure 9d shows the HPPC test that is described in the Freedom Car Battery Manual [52]. The ECE15 test, described in reference [53], is a more dynamic current profile, as shown in Figure 10e.

Experimental Results, Comparison, and Discussion
After obtaining the parameters for each model, detailed in Appendix A in Tables A1-A9, using the procedure described in Section 3 and identification current profiles described in Section 4, the output voltage accuracy and robustness analysis for the six supercapacitor models described in Section 2 is performed based on statistical metrics, such as relative error and root-mean-square (RMS) error.
Comparative results with identification current profile i 1 are illustrated in Figure 11a-d for the HPPC test and Figure 11e-h for the ECE15 test. Figure 11a,e show the experimental supercapacitor voltage and the voltages provided by the Stern-Tafel and Zubieta models. Figure 11b,f show the relative error between these models and the experimental data. The robustness and accuracy of the supercapacitor models are evaluated by means of different standardized test profiles, which include the Hybrid Pulse Power Characterization (HPPC) test and European Urban Driving Cycle (ECE15) for long-time responses. Figure 9d shows the HPPC test that is described in the Freedom Car Battery Manual [52]. The ECE15 test, described in reference [53], is a more dynamic current profile, as shown in Figure 10e.

Experimental Results, Comparison, and Discussion
After obtaining the parameters for each model, detailed in Appendix A in Table A1-A9, using the procedure described in Section 3 and identification current profiles described in Section 4, the output voltage accuracy and robustness analysis for the six supercapacitor models described in Section 2 is performed based on statistical metrics, such as relative error and root-mean-square (RMS) error.
Comparative results with identification current profile i1 are illustrated in Figure 11a-d for the HPPC test and Figure 11e-h for the ECE15 test. Figure 11a and Figure 11e show the experimental supercapacitor voltage and the voltages provided by the Stern-Tafel and Zubieta models. Figure 11b and Figure 11f show the relative error between these models and the experimental data.  Figure 11c and Figure 11g represent the relative error in percentage. Figure 11d and Figure 11h show the RMS error in mV. It shows that the Stern-Tafel model has lower error values in comparison with the Zubieta model. In any case, the relative error tendency with the time increase in both models, therefore the accuracy of both models identified with the i1 current profile is not proper.
Similar information is shown when current profile i2 is used to obtain the model parameters. Figure 12a-d depicted the obtained result for the HPPC test and Figure 12e-h for the ECE15 test. This current profile is applied to five out of the six models, with the exception of the Stern-Tafel model. In this case, the Series model is the best one, since it presents a reduced relative error that maintained with the time.   Figure 11d,h show the RMS error in mV. It shows that the Stern-Tafel model has lower error values in comparison with the Zubieta model. In any case, the relative error tendency with the time increase in both models, therefore the accuracy of both models identified with the i 1 current profile is not proper.
Similar information is shown when current profile i 2 is used to obtain the model parameters. Figure 12a-d depicted the obtained result for the HPPC test and Figure 12e-h for the ECE15 test. This current profile is applied to five out of the six models, with the exception of the Stern-Tafel model. In this case, the Series model is the best one, since it presents a reduced relative error that maintained with the time. Finally, the result obtained with the current profile i3, which is the most dynamic current profile, is depicted in Figure 13a-d for the HPPC test and Figure 13e-h for the ECE15 test. This current profile has been applied to the same models as current profile i2. Again, the Serie Model has the best performance, and even the obtained relative error is lower than using the previous current profiles. Nevertheless, the Parallel model, Transmission Line model and Thevenin model get good behaviors. The main conclusions obtained from these results are the following: • The greater complex identification current profile i3 gets greater accuracy for every model in which it can be applied.

•
In most cases, the Series model provides the minimum relative error. Finally, the result obtained with the current profile i 3 , which is the most dynamic current profile, is depicted in Figure 13a-d for the HPPC test and Figure 13e-h for the ECE15 test. This current profile has been applied to the same models as current profile i 2 . Again, the Serie Model has the best performance, and even the obtained relative error is lower than using the previous current profiles. Nevertheless, the Parallel model, Transmission Line model and Thevenin model get good behaviors. Finally, the result obtained with the current profile i3, which is the most dynamic current profile, is depicted in Figure 13a-d for the HPPC test and Figure 13e-h for the ECE15 test. This current profile has been applied to the same models as current profile i2. Again, the Serie Model has the best performance, and even the obtained relative error is lower than using the previous current profiles. Nevertheless, the Parallel model, Transmission Line model and Thevenin model get good behaviors. The main conclusions obtained from these results are the following:

13
• The greater complex identification current profile i3 gets greater accuracy for every model in which it can be applied.

•
In most cases, the Series model provides the minimum relative error. The main conclusions obtained from these results are the following: • The greater complex identification current profile i 3 gets greater accuracy for every model in which it can be applied.

•
In most cases, the Series model provides the minimum relative error.

•
If a simple and basic supercapacitor model has to be used, the best option is to use Zubieta model identified with the current profile i 3 .

Conclusions
This paper describes a parameter identification general procedure with a flexible and interactive interface used to build supercapacitor models in Simulink or Simscape. This procedure enables estimating the different models parameters based on the use of the Optimization Toolbox of Matlab ® . Once, the procedure steps are explained, the procedure is used to develop a comparative study of six commonly used supercapacitor models. In addition, the procedure enables using different identification current profiles, providing the possibility of analyzing the influence of three different identification current profiles in the accuracy and robustness of every model.
The experimental results obtained from the six models and three different identification current profiles, used to develop the study, show that both the model and the identification current profile are critical to obtaining good accuracy and robustness, which must be maintained over time.
From the comparison between the experimental results and the simulation results obtained using the model, it can be concluded that the greater complexity of the current identification profile, the greater accuracy and robustness of the model. In this case, the most complex identification current profile i 3 gets the best accuracy for every model in which it can be applied.
In a short simulation period, most models provide enough accuracy results. However, in a long simulation period the differences among models as well as among the current identification profiles