Supercapacitor Modeling: A System Identification Approach

Recently a great deal of attention has been given to supercapacitors (SC) due to their outstanding power densities and long cycling life. Their behavior has been extensively analyzed and tested through several modeling approaches. One common technique for modeling the dynamic operation of SCs is through an electrical circuit model (ECM). This article presents a new approach to identifying ECM parameters by applying subspace system identification (SSID) algorithms and incorporating coulombic efficiency. This novel application of SSID improves model accuracy by almost 50% in some cases compared to the literature procedures. This was done without manual tuning of the parameters, risk of non-convergence or any prior knowledge. The approach was validated at three different temperatures and with experimental data from an electric motorcycle. The resulting models are ready to be used as building blocks in a wide range of applications such as energy management systems, renewable power generation, electric vehicles, and microgrids.

I N OUR current energy landscape, energy storage systems (ESS) play a critical role not only in electric mobility and intermittent generation but also in power services, including frequency regulation and smoothing [1], [2], [3], [4]. These demanding applications require a variety of storage technologies, from pumped hydro storage (PHS) to Li-ion batteries (LIB), which differ in energy (kWh) and power (kW) capabilities. To maximize performance, lifetime, and cost and minimize their environmental impact, ESS systems may be equipped with many different storage technologies. These mixed ESS are called hybrid energy storage system (HESS) and they can be found in a wide range of applications [5]. One common feature of HESS is its high-power capabilities, which are usually achieved by using a SC. These components combine a superb cycling characteristic (approx. one half million cycles [6]) and outstanding power densities (500-5000 W/kg [7]). The applications of SCs include the regenerative braking of railways, electric buses, forklifts, pitch control systems in wind turbines and many more [8], [9].
There are 3 main types of SC: electrical double-layer capacitor (EDLC), pseudocapacitors and hybrid [10]. The EDLC is composed of two electrodes that store electrostatically adsorpted ions supplied by an electrolyte and isolated with a separator (non-Faradaic process). This adsorption is reversible and the voltagecharge characteristic is completely linear. Both electrodes are from the same material, usually activated carbon. These types of SCs are widely available on the market and a common example is Maxwell Tech. In the pseudocapacitor, the storage mechanism changes to pseudocapacitance, in which the charge is not stored due to the redistribution of ions but from an electron-transfer mechanism that occurs at the surface of the electrodes (Faradaic process). The electrodes are now redox-active molecules or metal oxides like RuO 2 that supplies electrons in the anode and receives protons in the cathode. They have greater charge 0885-8969 © 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
densities than EDLCs and their voltage-charge relationship is mostly linear inside the voltage operating window [11]. Finally, hybrid SCs combine both types by using an EDLC electrode and a pseudocapacitive/battery-like electrode. They can be symmetrical or asymmetrical depending on the combination of materials used for the electrodes. Their motivation is to enhance the energy density of the SC system, which they have successfully done increasing it up to 20-30 Wh/kg [11]. A common asymmetric hybrid is a combination of a non-faradaic activated carbon positive electrode with a faradaic Lithium-titanate (LTO) negative electrode. The hybridization of the SC device results in a non-linear relationship between open-circuit voltage (OCV) and stored charge. Their modeling and characterization has been extensively studied either via electrochemical impedance spectroscopy (EIS) techniques [12], [13], [14], [15], [16], [17] or time-domain response experimentation and simulation [18], [19], [20], [21]. Both approaches aim to build different equivalent circuit models which are crucial to describe the state of the cell during operation. Additionally, other data-driven methods and models based on artificial neural networks (ANN), fuzzy logic, or other so-called "intelligent methods" are found in the literature not only for SCs but for other storage devices like LIBs [17], [22], [23].
Circuits based on EIS [12], [13], [14] are parametrized through complex non-linear least squares (CNLS) optimization algorithms, such as Levenberg-Marquardt, Powell's or NP Simp. The procedure for fitting EIS spectra requires intensive interaction between humans and software, due to unfeasible local minimums, non-convergence, and other similar problems. Alternatively, step-response tests are also used to obtain ECMs [18], [19], [21]. These methods rely on manually tuning at least one time constant, which also restricts the full automatization of the fitting process. Thus, there is a need for a fully automatic parameter identification procedure with minimum human intervention.
Additionally, none of the above-cited references take into account the device's coulombic efficiency η t [24]. The references represent losses in the form of resistances, which account for the power losses due to the Joule effect. However, coulombic efficiency represents the ratio of charge output to charge input and is an important indicator of unwanted electrochemical processes. The discrepancy between the output and input charge comes from the presence of Faradaic processes in the device that consume charge during operation. Hence, η t represents how reversible the storage mechanism is, whereas the power losses are represented in the thermal efficiency rather than the coulombic efficiency. This link between electrochemical and electrical models is not present in the models previously mentioned and is specially important in aged SC.
These drawbacks complicate the implementation of a common framework for testing, modeling, simulating, and implementing ESS comprehensively and integrally. Such a framework should include device management systems (DMS) and energy management systems (EMS) that supervise the performance of the storage device-an SC in our case-as well as the whole system, including converter, grid interface, and other storage devices.
In this paper, we present a novel systematic approach to identify equivalent circuit models (ECMs) that describe the dynamic performance of SCs, following the system identification process by Ljung [26] and the algorithms in [27]. The main contribution of this paper is the application of SSID algorithms that naturally solve the shortcomings associated with the initialization of the algorithm, non-convergence, and manual tuning with improved accuracy in experimental validation. Since the identification is solved in a least-squares fashion (non-iterative) the resulting parameters are the optimal solution with respect to the gathered data. The algorithm also overcomes identifiability issues, because the mapping produced between data and model parameters is injective [26], [28]. The second contribution is the incorporation of the coulombic efficiency η t which yields a great improvement in model accuracy, especially for aged SCs. The third contribution is that the produced general workflow, an its resulting state-space model, can be easily integrated and implemented in high-level simulations and embedded systems. Including easy integration with functions such as State of Charge (SOC), State of Health (SOH), thermal management, and power balance.
This paper is organized as follows: Section II reviews previous equivalent circuit models presented in the literature and the tests and fitting procedures used, including their advantages and shortcomings; Section III presents our approach and algorithm; Section IV describes our tests, experimental setup, model fit, and validation; and finally Section V presents the conclusions and shortcomings of this approach.

II. REVIEW OF EQUIVALENT CIRCUIT MODELS
Because supercapacitors are still novel non-consolidated devices [29], the literature contains disparities about how to model them. Most of the models in the literature have different goals and come from different system identification processes. The two most common tests performed on the SC are EIS and time-domain current profiles (i.e. charge-discharge tests and step response tests). The first method is a frequency-domain analysis of the device's impedance, whereas the second characterizes the time-response of the system. Both methods have led to two different equivalent circuit models. The first approach is used in [12], [13], [14], [30], [31] to derive circuit models which contain non-linear components such as Warburg impedances, constant phase elements (CPE), Gerischer elements and so on. From the distributed circuital approximation of the pore (see Fig. 1) the pore impedance Z p is described by: where R el is the electrolyte's resistance and C dl is the capacitance of the double-layer, as per the distributed equivalent circuit shown in Fig. 1. The full ECM also includes a series resistance R 0 , and a high frequency inductor L hf that models the cables and connection to the instrument (see Fig. 2(a)). This model was extended in [25] and [30] to account for the ageing of the SC. This was done by replacing the distributed double-layer capacitances in the finite-length diffusion reflective boundary model with distributed CPEs. Thus, the pore impedance changes to: where γ accounts for the imperfection of the pore and is in the order of 10 −3 [25], and A dl substitutes the C dl in terms of the magnitude of the diffusion process, with units of [F.s −γ ]. The γ changes the 45 • angle in the Nyquist diagram. The second approach uses the time-domain response of the SC as the primary source of analysis [18], [19]. These authors do not appeal to electrochemical arguments to derive their model but instead use common electrotechnical methods to derive a non-linear 3-branch parallel circuit that describes the 2 timeconstants present in the system. Fig. 2(c) presents a simplified model that accounts for self-discharge with R leak , a fast dynamic mode represented with an intermediate branch with R 2 and C 2 and a non-linear capacitance C(v C ) that depends on its own voltage drop. This last C is the main capacitance of the supercapacitor.
Both models were synthesized in a unique dynamic model by Tironi [21], see Fig. 2(d), one especially aimed at power electronics applications for a frequency bandwidth up to a few kHz.
Finally, a "Basic model" comprising the main capacitor, a series resistance, and a leakage resistance is presented as a reference (see Fig. 2(b)). The model can be easily fit using the information from the manufacturer's datasheet and simple time-domain tests. In the following sections the model presented by Buller, Fig. 2(a), will be called "Series model," whereas the one presented by Faranda, Fig. 2(c), will be called "Parallel model". The model presented by Tironi, Fig. 2(d), will not be discussed any further in this work, because its complexity does not pay off in terms of performance and insights. At last but not least, it is worth noting that the presented Z p makes the series model non-linear in its dynamics, whereas the voltage dependent capacitor of the parallel model makes it non-linear in steady state. The first type of non-linearity is usually overcome by discretizing the number of poles [12] and the second one becomes more important in pseudocaps and hybrid SCs as we mentioned in Section I.
None of the three ECMs described thus far follow a systematic approach to identifying their descriptive parameters nor include the coulombic efficiency η t . Buller [12] presents a fit using EIS data, for which the most commonly used approach is solving CNLS problems using Levenberg-Marquardt or Gauss-Newton algorithms [32], which have the drawback of heavily depending on the initial conditions to converge on a suitable solution. Methods like genetic algorithms [33] can overcome this at the expense of a heavy computational burden, which complicates its implementation. Other works [18], [21] use a simpler approach but hand-tune one of the time constants and use manufacturer's data, which might not always be available. At last but not least, η t is an important electrochemical parameter that can be used to diagnose the SOH of the device because a drop in η t can be an indicator of unwanted side-reactions or inventory loss in the electrodes.

III. IDENTIFICATION PROCEDURE
The drawbacks associated with CNLS problems and manual fitting can be overcome by applying subspace system identification [27]. In this approach, the main goal is to estimate the state-space model parameters of the system under study. This is done through several projections over measurement subspaces, created with the data collected during experimental tests. In our application, this is very useful given that electrical circuits have fairly straightforward representations in the state-space.
To carry out this approach several open-loop tests were performed (see Fig. 3(a)), in which the input u and output y signals were recorded. The sampled measurements were then processed through the SSID algorithm to produce a model. In order to correctly identify the model, the first step was to define a suitable excitation profile u (see Fig. 5(a)). This signal must excite all the relevant dynamic modes of the system [27], thus sweeping all the voltage operating window and containing a considerable amount of harmonics [21]. Secondly, the profile must follow particular constraints to correctly calculate the first set of model parameters, which is discussed in Section III-A.

A. Model Structure and First Output Estimate
The proposed algorithm is similar to the one presented by Plett for LIBs in [34], [35], [36], [37], since both devices can be modeled in the state-space. From a physical point of view, this is natural, given that the dominant transport mechanism in both storage devices is diffusion. Hence the appearance of Warburg impedances in both LIBs and SCs [12], [13]. Their main difference is that the SC models do not have open-circuit voltage sources or hysteresis states when compared to the so-called Enhanced Self-Correcting (ESC) model presented by Plett [37]. The general structure of the estimated state-space model has the following form in the continuous time domain: with x being the states vector and A, B, C and D the model matrices, w being the process noise, and v the sensor noise.   The system matrices change in the different ECMs. Table I presents the different structures used. The state-space model was discretized for implementation on Matlab.
Given that all three ECMs ("Series," "Parallel," and "Basic") have some basic observable behavior, they are grey-box models rather than black-box models. Thus, the first set of parameters, the main capacitance, and coulombic efficiency, is identified before applying SSID algorithms. These parameters vary depending on the model. Coulombic efficiency η t is identified in all three models, and the main capacitance is identified in all except the Basic model.
The identification algorithm is presented in Algorithm 1. The algorithm starts by taking the input data structure, the initial model structure, and the number of poles, np, to be identified. The first argument contains the time-series data of the tests, the second one is general information about the SC to be tested (rated voltage, for example) and the last argument varies depending on the type of model. Then, for all test temperatures T k , we proceed to identify a complete model.
The identification procedure starts by calculating both the coulombic efficiency, η t , and the main capacitance, C, having the following model for them: where Q dis,T k is the total charge discharged in all steps at test temperature T k , Q chg,T k is the total charge charged in all steps at test temperature T k , η t,T k is the coulombic efficiency at test temperature T k . This way of calculating η t takes into account the variability of our current since the charges Q chg/dis,T k come from integrating it along the time steps of the experiment. The only needed restrictions are that either an individual charge-discharge cycle has to be symmetrical in current [24], or by ensuring that the initial state of charge SOC(t = 0) is equal to the final SOC(t = t f ) at the end of the test. The latter is what we did in our experiments, since every time a train of pulses finished passing through our SC we then took the SC to the initial SOC(t = 0). Other methods for estimating η t have been presented for different electrochemical devices. For instance, for vanadium redox flow battery (VRFB) Xiong et al. [38] estimated it from OCV and leakage resistance whereas for LIBs Yang et al. [39] estimated η t by fitting life-cycle data in a least-squares fashion.
The main capacitance C is a linear function of the voltage v C and may be fitted with the methods presented in [18] or [15] which use simplifying assumptions in the time-domain and the EIS/frequency-domain respectively. In the procedure described here, the relationship was estimated by taking each charging half-cycle and computing a capacitance signalĈ as: where i chg is the charging current, which is fairly constant during each interval, v t is the terminal voltage of the SC, with the hat indicating estimation. This is because the voltage used is the terminal v t not the internal state v C which can not be measured directly. The estimated capacitance signalĈ is then mapped against the terminal voltage, as inĈ = f (v t | i=i chg ), and by using linear regression on the data we identify C 0 and k V . It is important to note that the voltage sweep must cover the whole voltage window of operation to correctly estimate both parameters.
Once η t and C are identified, the first output estimate y est1 is computed. In all models, except for the basic model, this is done by computing the voltage of the ideal capacitor C. This ideal voltage is the steady-state response of the SC, without accounting for the dynamics of the system. Next, the first output error, y err1 , is calculated by comparing y est1 against the output measurements y k .

B. Identifying A RC
Once this first identification is made, the order of the model can be reduced. The reduced model describes the different relaxation times that account for the dynamics of the SC. The information is contained in the A RC matrix, a reduction of the original A matrix in the state-space model; the same reduction applies to the other model matrices. By using y err1 and u we can identify A RC and its companion matrices B RC , C RC and D RC . In this particular case, we only identify A RC because we have proposed different ECMs that relate the different matrices between them. Thus, all the parameters to be identified are contained inside the state transition matrix A RC .
To obtain this matrix, we follow the combined deterministicstochastic problem presented in [27], which starts with: The output equation (9) shows how to generate the future outputs sequence Y f through the future states sequence X f = X d f + X s f , which is composed of a deterministic and stochastic component, the sequence of future inputs U f , process noise N w f and sensor noise N v f . Subscript p refers to past, f to future and superscripts d and s to deterministic and stochastic respectively. To process the states, the extended observability matrix O i is required, with i being the number of block rows. Finally the extended Toeplitz matrices built with Markov's parameters Ψ i and Ψ w i process both the input and the process noise. These last 3 matrices (O i , Ψ i and Ψ w i ) are built with the system matrices A RC , B RC , C RC and D RC (see Appendix). The implemented algorithm aims to estimate a state sequence and approximate the system matrices with it. The algorithm goes as follows: 1) Calculate the following projections: where W p is the combined past input-output Hankel matrix (see Appendix) and the / symbol denotes the projection operation. Z i is the optimal prediction of the output data Y f taking into account the future inputs U f and the past W p . The noise terms in (9) are eliminated through orthogonal projections Z i and Z i+1 on to the combined input-output subspace formed by W p . Finally the oblique projection ξ i is used to clear out the future inputs term in (9). 2) The oblique projection ξ i holds the following property: where O i is the extended observability matrix andX i is the Kalman filter state sequence. Incorporating weighting matrices W 1 and W 2 , its singular value decomposition is: Different methods exist to estimate O i depending on the weighting matrices to be used. Our proposed algorithm follows the robust implementation in [27] which uses only one weighting matrix: where Π U ⊥ f is the orthogonal projection onto the future input subspace U f . The SVD of the weighted oblique projection is computed as: 3) Use the model order np to appropriately reduce U and Σ into U 1 and Σ 1 , which have np rank. 4) Compute the observability matrices as:

5) Estimate A RC withÂ RC from:
and recompute the observability matrices O i and O i−1 .
To ensure the feasibility of the model the eigenvalues of A RC have to be real and less than 1. 6) Solve forB andD from the minimization problem: 7) Finally we compute the residuals ρ w , ρ v and covariances of the output estimates with (17) and (20): For our case study we only solve this algorithm until step 6 because matrices B and D will be identified in the next subsection with a more simple circuit analysis approach.

C. Remaining Parameters and Final Output Estimate
After the model dynamics are identified we still need to decompose the information of the system's poles into the different resistors and capacitors of the ECM. To do that, the x states of the models are simulated using the reduced order model. Then the output equation, (5), has to be reorganized taking into account the missing information of the model. For example, in the series model, Fig. 2(b), (5) takes the form of: since the first output estimate is y est1 = v C in (21) turns into: where i and i R j are known and the resistances are unknown.
Since the problem is handled in the discrete time-domain (22) can be solved as a non-negative least squares (NNLS) problem [40]. The NNLS problem is: whereŴ is the parameter estimator and H † is the pseudoinverse matrix containing the simulated states x and the input measurements u. In the case of the series model they are: For the rest of the models the output equation, (5), was restructured and solved in the same fashion.
Finally the last estimate and error are: (25) y err2 = y k − y est2 (26) Summing up, to identify the whole model, the first estimation y est,1 is done with a classical model with minimum assumptions. Then the remaining parameters are fitted by solving (9) with 3 projections that clear out the noise and then the perturbation. Once the observability is estimated, identifying A RC comes naturally. The x states have to be simulated to restructure (5) and implement an NNLS problem that estimates the remaining parameters. With them, the final estimate y est,2 comes straightforward.
Notice that the method does not have any of the disadvantages seen in the reference models. A minimum number of assumptions is used to estimate the parameters. The proposed approach does not depend on hand tuning any parameter, it is not susceptible to non-convergence nor any of the shortcomings mentioned in Section II. In particular, for our problem the solutions of the least squares problems (A RC estimation and resistance estimation) have to be feasible. This means that the identified poles need to be real and within the unit circle and the resistances must be real and positive. When these constraints are active the global optimality is compromised. For the unconstrained case the achieved solution is always a global optimum with respect to the measured data since the least-square problems presented are linear. For infinitely long datasets this is demonstrated by Van Overschee [27], Ch. 4 and Appendix A. The impact of finite datasets on subspace identification methods is discussed by Markovsky [41]. In the authors' experience the performance of the used algorithm is sufficient and modifications to address the sample size are not necessary.

IV. EXPERIMENTAL SETUP & RESULTS
Experimental validation was carried out in two parts. First, the SC models were identified with the reference procedures [12], [18] and SSID. All the models were identified using the profile of Fig. 5(a) except for the reference Series model which was fitted following the procedure presented in [12]. Thus, additional EIS experiments using standard electrochemical equipment (Biologic VMP3) were performed. Measurement was done with a bidirectional power supply (PSB 9000 3 U, from Elektro-Automatik), a datalogger (Yokogawa DL850EV) and a computer were used for processing the data. Communication between different devices was done using a USB memory. More sophisticated communication solutions, such as Ethernet, are up to the user (factory, laboratory, etc.). The term "automatic process" here refers to the fact that the SSID algorithm does not need human intervention to achieve a suitable solution. This a crucial difference with CNLS [12] optimization and the hand-tuning needed in [18].
All six identified models were then simulated against the identification profile, see Fig. 5(b). The identification was done for several test temperatures, but since the resulting identification  plot of the models is mostly the same for every temperature only the result at 25 • C is shown. Once this was done the identified models were validated against data gathered from a real electric motorcycle [42] (see Fig. 6).
The SCs under study were BCAP3000 cylindrical SCs with activated-carbon electrodes and organic electrolyte from Maxwell Technologies. These cells can be regarded as EDLCs.

A. Model Identification
The experimental setup, scheme and equipment for model identification are shown in Figs. 3(a) and 4(a). A controlled current source was used to inject the current profile shown in Fig. 5(a), and to produce the voltage drift we injected a train of 10 A pulses into a node with an SC in parallel with a shunt resistive load (see Fig. 4(a)). Every time the current source went to zero, the circuit breaker opened and the current flow was let through the loop.
For the models presented in Table I the resulting model identification is presented in Table II. For the series model the number of RC parallels was limited to np = 2 in the reference fit and to np = 1 in the SSID algorithm. This is because the contribution of the other poles is indistinguishable from the noise, thus the SSID algorithm could not identify more significant poles. Given that the pore impedance has been discretized the parameters  Table II are numbered following the τ j poles, with j = 1, 2. . .np.
One important comment is that for 0 • C and 35 • C the reference methods resulted in unfeasible parameter sets for the parallel model. Thus, the model is not included during validation. The SSID algorithm had no problem identifying feasible solutions for those temperatures.
At T = 25 • C, the SSID approach had difficulties finding feasible solutions for the last resistance R 2 of the parallel model. This last element represents the leakage resistance of the cell. This is because the time constant of the self-discharge is so big, that to identify it correctly we would need to let the SC rest for at least 72hs according to the manufacturer. Thus, the final output estimate has a reasonable root mean squared error (RMSE) but an unsuitable solution.
The resulting fits are presented in Fig. 5(b) for the T = 25 • C case. The reference fits and the SSID parallel fit drift away from the measurements as the test progresses, while the series and basic models with SSID fits represent the real system much better, as we can see in the inset graph of Fig. 5(b). This is reflected in the total RMSE (Table III), in which our proposed SSID approach beats its competitor by almost one order of magnitude. This has serious implications during implementation, as the computational burden on the DMS is lessened due to the less frequent update of the parameters. The rest of the identification temperatures show the same trend, so there is no use on showing them.

B. Validation -Electric Motorcycle
To validate our models, they were tested with a real-life application. In [42] data was gathered from a full-electric  motorcycle equipped with a HESS composed of SCs and LIBs (see Fig. 4(b)). The HESS splits the current demanded by the electric motor into the battery current and the supercapacitors' current. The HESS aims to smooth the battery current so the AC variable current is provided by the SCs. The signals for validating the proposed technique are the current and voltage of the SCs. During this test, the current of the branch and the terminal voltage of the stack were measured with different sampling frequencies depending on the temperature (from 1 Hz to 5 kHz). At T = 25 • C a stack of SCs was used whereas for the rest of the validation temperatures an individual cell was used. When the stack was used the terminal voltage was assumed to be equally divided between all series cells. Fig. 6 shows the simulations and the experimental data. Following the results of Section IV-A the reference Parallel model was not included in for T = 0 • C and T = 35 • C. Moreover, the reference Series model was also discarded for T = 35 • C because its voltage estimation was out of range.
The robustness of our approach can also be seen in Fig. 7 and Table III. In short, all models decrease their metrics when comparing validation against identification, which is reasonable given that they were fitted with a different signal. For 25 • C, the SSID parallel model substantially improves its performance and the SSID series model keeps beating its reference counterpart, but the basic ECM with our fitting approach loses against the reference fit. Even so, the ref. series model is not consistent across the different temperatures. At 0 • C it is the best performer, but its voltage error increases tremendously with temperature. This adds to the case of the ref. parallel model, which is also inconsistent for the data gathered at different temperatures. On the contrary the SSID method maintains tight first quartiles below 50 mV for all models.
In the authors' experience, the time difference between identifying a single ECM and multiple ones with SSID is not significant once the software has been developed. Hence, to find the appropriate model, our suggestion is to identify and validate all of the models and pick the one with the lowest RMSE during the validation test. For the model to be accurate during operation, both the identification and validation profiles should correspond to the final application. This means that the pulses and electric motorbike test profiles, identification and validation respectively, will change depending on the desired final application. The final model type depends strongly on the application and the type of SC used (EDLC, pseudocap or hybrid SC).
Our proposed technique is more robust than the referenced traditional identification procedures for a broader range of operating conditions. Additionally, it does not require prior knowledge of the system since the number of poles np is the only parameter necessary to run the identification. Hence, in the case of missing data from the manufacturer, an accurate model of the storage device can still be obtained. The SSID models consistently preserve their accuracy, ensuring good confidence intervals without any prior knowledge.
Nevertheless, our strategy still has its drawbacks, such as the numerical conditioning of the NNLS problem for the parallel model. This can be overcome by applying different adaptive control techniques during operation, such as dual extended Kalman filtering [36] or by actually solving step 6 of Section III-B. Van Overschee [27] even proposes two more algorithms to estimate the system matrices. The major drawback is that a significant amount of engineering hours must be spent on software development/adaptation of current libraries.
We can thus state that our approach is solid enough and can accurately describe the behavior of a real-life application.

V. CONCLUSION
In this work, we synthesize different modeling strategies in one integral approach by applying subspace system identification. This new application of a theoretically proven and industry-tested method has shown an important improvement in accuracy and flexibility with respect to the methods present in the literature. We have successfully overcome the disadvantages described in Section II by implementing an automatic grey box model based on simple electrochemical considerations.
The resulting model can be easily integrated into device management system algorithms for several applications in an almost automatic fashion. This yields considerable advantages in developing the HESS and devices themselves by being a bridge between the electrochemical/physical models and the implementation of the final product's control systems.
The proposed approach remains to be implemented in hybrid SCs, which exhibit different non-linearities not present in EDLCs.

APPENDIX NOMENCLATURE
This appendix serves a reference on the nomenclature used in Section III.
The model states are a linear combination of deterministic, x d k , and stochastic, x s k ,components: and the output signal is: The state-space model in the discrete time has 2 components, the deterministic one: y d k = C.x d k + D.u k (30) and the stochastic one: The input u k block Hankel matrices U 0|2i−1 , U p , U f , U + p , U − f are divided into past and future blocks indicated by sub-indices p and f defined by the number of block rows i. This last index i is a user defined parameter. Superscripts + and − stand for "add/delete one block row". These matrices are computed as: The output y k block Hankel matrices are defined similarly. For the input noise w k and output noise v k the Hankel matrices N w are also defined in the same way. The input-output Hankel matrices are defined as: The state sequences are defined as: where i denotes the subscript of the first element of the state sequence. Hence, we note that the past and future state sequences are: The extended (i > n) observability matrix is: (39) and the controllability matrix is: The lower block triangular Toeplitz matrices are build with Markov's parameters. For the deterministic part: and for the stochastic part of the problem: (42) Finally the covariance and cross covariance matrices are defined as: