**1. Introduction**

Energy storage systems (ESSs) are among the most critical components concerning the full adoption of renewable energy sources and electric transportation [1]. In such applications, energy management systems (EMSs) are required for coordinating the operations of systems with multiple energy generation resources, as often is the case with microgrids and hybrid ESSs, which employ multiple complementary energy storage subsystems. For the energy management of ESSs, the most commonly desired objectives are: improved energy efficiency, extended lifetimes for storage elements and compliance with constraints of internal energy storage modules, e.g., no over-charging and no over-discharging [2]. In the case of rechargeable batteries, those goals typically require accurate monitoring; state of charge (SoC) and state of health (SoH) estimations; and charging control, which can be accomplished by using dedicated battery management systems (BMS). The methods employed for these operations range from experiment-based ones, such as incremental capacity [3], to data-driven ones, including machine learning [4]. State estimation and control algorithms in a BMS may require battery models capable of approximating the battery's response under given operating conditions [5].

Equivalent circuit models (ECMs) are preferred over electrochemical or empirical models, as they approximate the dynamic behaviour of a battery with relatively high accuracy [6], while offering simplified descriptions of the complex physical and chemical processes occurring within batteries, by representing them with a set of lumped elements, including resistors, capacitors and inductors. ECMs have been widely used as parts of battery state estimation and charging control strategies. The adoption of ECM-based optimal state estimation, by directly solving the model-constrained optimisation problem [7] or using techniques such as the Kalman filter and all its variants [8], has enabled the implementation of online SoC and SoH estimations in EMS frameworks. The implementation of

**Citation:** Ospina Agudelo, B.; Zamboni, W.; Monmasson, E. A Comparison of Time-Domain Implementation Methods for Fractional-Order Battery Impedance Models. *Energies* **2021**, *14*, 4415. https://doi.org/10.3390/en14154415

Academic Editor: Hugo Morais

Received: 1 May 2021 Accepted: 18 July 2021 Published: 22 July 2021


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

online optimal charging strategies, known for their high computational burden, has also been enabled by algorithms, including ECMs, by profiting from the accuracy–complexity trade off offered by such models [9].

The parameters of an ECM can be fitted from voltage and current data obtained during specific operating conditions [6]. However, these circuit components are often insufficient for modelling the dynamics of electrochemical processes such as charge and mass transfers and double layer capacitance in a battery, due to the spatial distribution of those processes [10]. This lack of physical significance may compromise the identification of relationships between health estimation and ECM parameters. Said drawback may be overcome by substituting capacitors in the ECM with constant phase elements (CPEs), defined in the frequency domain and analytically derived from the electrochemical principles of the diffusion processes. Such elements are often used to fit electrochemical impedance spectroscopy (EIS) data [11]. The reduced number of parameters in ECMs using CPEs instead of capacitors is mainly attractive for SoH estimation approaches based on the analysis of variations in the parameters associated with specific electrochemical processes [12].

The direct time-domain implementation of such models is particularly challenging, because they represent dynamic systems with non-integer-order derivative operations, hence the name fractional order models (FOM) [10]. The time-domain implementations of FOMs' electrical responses have been employed in applications such as online SoC estimation using fractional-order (FO) Kalman filters [13] and a time-domain characterisation of battery diffusion dynamics [14]. Three main implementation approaches have been identified: the multiple RC (mRC) circuit [15], the high-order integer transfer function [16] and the Grünwald–Letnikov (GL) fractional derivative [17] approximations. However, the literature is rather obscure about the motivation for the selection of one approach over the others. This work fills this gap by presenting a comparative study of the three approaches. The comparisons are made in terms of:


We used a simulation environment with a set of FOMs representing the middlefrequency range impedance of Li-ion batteries, normally represented by the so-called ZARC elements. Our study identifies which of the analysed approaches offers the best compromise between accuracy and computational burden for applications such as real-time simulations. The second goal was to understand for which of the three cases the timedomain identification of the FOM leads to a correct frequency-domain response, keeping the impedance model meaningful. This analysis may serve as a guide for the selection of implementation approaches for FOMs in EMS applications.

The paper is organised as follows. Section 2 introduces the battery impedance models based on constant phase elements, Section 3 presents, in detail, the three time-domain implementation approaches studied in this work. Then, the results of the comparative analysis in terms of accuracy and complexity are presented in Sections 4 and 5, respectively. Section 6 presents a time-domain identification study performed for the three implementation approaches. Finally, the conclusions are drawn in Section 7.

### **2. Battery Impedance Models**

EIS analyses the impedance of battery cells in a specific range of frequencies, drawing conclusions about internal electrochemical processes with different time constants. A typical EIS test is performed by injecting a sinusoidal voltage or current signal; acquiring the cell response; and computing the impedance through the amplitude ratio and the phase difference between voltage and current signals. That experiment is repeated for several

frequencies in a selected frequency range. Then, by calculating the impedance, a Nyquist plot, often called impedance plot, can be generated. If this is performed for multiple battery operating conditions (temperature, SoC and dc operating point), it is possible to visualise and quantify their influences over the frequency response. For a Li-ion battery cell, a typical impedance plot is shown in Figure 1 [18].

**Figure 1.** Qualitative impedance spectrum of a Li-ion cell and an ECM approximating it.

In the qualitative impedance plot shown in Figure 1, it is possible to identify four sections that can be associated with particular electrochemical processes [18]. In the first part, an inductive behaviour can be seen at high frequencies, related to the inductive reactances of metallic elements in the cell and wires. The presence of an ohmic resistance is revealed by the intersection with the real axis at a nonzero value. This corresponds to the sum of the current collectors, active material, electrolyte and separator resistances. The semi-circle-like section represents the double layer capacitance and charge transfer processes at the electrodes. Finally, at low frequencies, the main effect corresponds to the diffusion processes of the active material of the electrodes, which manifests as a section with a constant slope in the impedance plot. It is worth noting that measured spectra often show variations with respect to the qualitative curve presented. For example, the number of semi-circles can be increased to two or more, or the inductive part can exhibit a slope with an increasing real part. Moreover, the semi-circle section may tend to present a depression at its mid-point (non-constant radius) [19].

For a relatively accurate and meaningful reproduction of the impedance spectrum of a battery cell, the ECMs presented in Figure 1 can be considered [20]. The behaviour of the cell at high frequencies is represented by an ideal inductor *L* and the resistor *Rs* represents the ohmic resistance of the cell elements. The semi-circle section can be represented using a ZARC element, which corresponds to a parallel connection between a resistor and a CPE. The impedance of a CPE, *ZCPE*, is presented in:

$$Z\_{\rm CPE}(\omega) = \frac{1}{Q(j\omega)^{\Phi'}} \tag{1}$$

where *Q* corresponds to a generalised capacitance, *φ* to the depression factor, *ω* to the angular frequency and *j* to the imaginary unity [21].

The low-frequency response is represented by a Warburg impedance, which models a semi-infinite linear diffusion process. The Warburg impedance, namely, *ZW*, is characterised by a 45◦ phase, and therefore can be also represented by a CPE:

$$Z\_W(\omega) = \frac{A\_w}{\sqrt{\omega}}(1 - j) = \frac{\sqrt{2}A\_w}{\sqrt{j\omega}} = \frac{1}{Q\_w(j\omega)^{\Phi\_w}}\tag{2}$$

where *Aw* is the diffusion coefficient [21], and the equivalent generalised capacitance and depression factor for the Warburg element are given by *Qw* = 1/( <sup>√</sup>2*Aw*) and *<sup>φ</sup><sup>w</sup>* <sup>=</sup> 0.5, respectively.

CPEs are used instead of capacitors because when they are connected in parallel with a resistor forming a ZARC element, *φ* represents the depression often observed in the semi-circles in the impedance plot of a Li-ion cell. The factor *φ* can take values between zero and one; when *φ* = 0, the ZARC element represents only an ohmic resistor, and if *φ* = 1, the response of an RC element is obtained. The ZARC impedance is represented by:

$$Z\_{ZARC}(\omega) = \frac{R\_{p1}}{1 + R\_{p1}Q\_1(j\omega)^{\Phi\_1}},\tag{3}$$

where *Rp*1, *Q*<sup>1</sup> and *φ*<sup>1</sup> characterise the resistance and CPE parameters of the ZARC element. The overall impedance of the circuit presented in Figure 1 is then defined by:

$$Z\_{\rm bat}(\omega) = j\omega L + \mathbb{R}\_s + Z\_{\rm ZARC}(\omega) + Z\_{\mathcal{W}}(\omega). \tag{4}$$

For the sake of simplicity, the discussion that follows focuses on the computation of the voltage *v*(*t*) of the ZARC element for a given current signal *i*(*t*). The presented results can easily be extended to the Warburg element using similar principles, and from there to the whole battery's impedance.

### **3. Time-Domain Implementation of the ZARC Element Response**

In the framework of battery modelling, three main approaches for approximating the time-domain responses of elements with FO transfer functions have been proposed.


Flowcharts for the three implementation approaches are given in Figure 2. In the three approaches, the time response of the FOM is approximated by a set of discrete equations, corresponding to discrete state-space representations for the mRC and OU approaches and to a difference equation for the GL approach. These discrete implementations were evaluated in our simulation environment in terms of accuracy and complexity, as presented in the following sections. The main differences between the approaches lie in the parametrisation

processes. For the GL approach, the coefficients in the difference equation are computed directly from the impedance parameters using predefined mathematical relationships, which are introduced in Section 3.3. As illustrated in Figure 2, the other two approaches require an initial approximation of the FO transfer function with an integer-order system, for which a discrete state-space representation can be obtained. The main difference between the mRC and OU approaches is the structure of the integer-order transfer function and the method used for the approximation of the FO transfer function, which are introduced in Sections 3.1 and 3.2, respectively.

In this section we will focus on defining difference equations for the computation of the ZARC voltage *v*(*t*) for each one of the mentioned approaches. The resulting expressions have the form given by:

$$
\upsilon[k] = \mathbb{C}\mathfrak{x}[k] + \mathrm{Di}[k].\tag{5}
$$

Equation (5) allows one to compute the ZARC element voltage in the discrete instant *k* as a function of the input value *i* and the vector *x*. The vector *x* may be the states' vector or a set of previous values of *v*, depending on the implementation approach. The values of vector *C* and scalar *D* are functions of the parameters of each approximation, as is discussed in the following subsections.

**Figure 2.** A summary of the FOM implementation approaches considered.

### *3.1. Approximation 1: Multiple RC Approach*

For the approximation of the ZARC element impedance, series connections of multiple parallel RC circuits have been considered in the literature [15]. The values of the components for this kind of approximations are normally fitted from time measurements directly, by minimising the differences between experimental data and the model voltage output. However, for the case in which the initial point is an impedance model, the estimation process requires fitting the impedance spectra. The idea is to approximate the transfer function of the ZARC element with a set of RC parallel elements:

$$Z\_{ZARC}(\mathbf{s}) = \frac{R\_{p1}}{1 + R\_{p1}Q\_1\mathbf{s}^{\Phi\_1}} \approx \sum\_{h=1}^{n\_{R\mathcal{C}}} \frac{R\_h}{1 + R\_h\mathbf{C}\_h\mathbf{s}} = Z\_{R\mathcal{C}}(\mathbf{s}).\tag{6}$$

In (6), a transfer function representation is used for the impedance, using *s* as the Laplace complex variable. The approximation of the ZARC impedance *ZRC*(*s*) employs a set of *nRC* RC branches. The parameters *Rh* and *Ch* represent the resistance and capacitance of the *h*-th RC parallel branch in the approximation of the ZARC impedance presented in (6), accounting for 2*nRC* parameters to fit.

The values for the resistance and capacitance in (6) can be computed by minimising the difference between the ZARC impedance, given in (3), and the mRC approximation. One method for solving this minimisation problem is presented with detail in Appendix A. This method is employed from this point on during the parametrisation stage of the mRC-based implementations.

Once the values for *Rh* and *Ch* have been fitted, the continuous time response of the mRC circuit can be written in state-space representation:

$$
\dot{\mathbf{x}}\_{RC}(t) = A\_{RC}\mathbf{x}\_{RC}(t) + B\_{RC}i(t) \tag{7}
$$

$$
\boldsymbol{w}\_{\rm{RC}}(t) = \mathbb{C}\_{\rm{RC}} \mathbf{x}\_{\rm{RC}}(t) + D\_{\rm{RC}} i(t), \tag{8}
$$

where the states' vector *xRC* contains the RC elements voltages, *x*˙*RC* represents the states' derivatives and *vRC* is the approximation of the ZARC voltage using the mRC approach. Additionally, the matrices *ARC*, *BRC*, *CRC* and *DRC* are *nRC* × *nRC*, *nRC* × 1, 1 × *nRC* and 1 × 1 matrices, respectively, given by:

$$A\_{RC} = \text{diag}\left(-\frac{1}{R\_1 C\_1}, -\frac{1}{R\_2 C\_2}, \dots, -\frac{1}{R\_{nRC} C\_{nRC}}\right) \tag{9}$$

$$B\_{R\mathbb{C}} = \begin{bmatrix} \frac{1}{\mathbb{C}\_1} & \frac{1}{\mathbb{C}\_2} & \cdots & \frac{1}{\mathbb{C}\_{R\mathbb{C}}} \end{bmatrix}^\top \tag{10}$$

$$\mathcal{C}\_{\rm RC} = \begin{bmatrix} 1 & 1 & \cdots & 1 \end{bmatrix} \tag{11}$$

$$D\_{R\Gamma} = 0.\tag{12}$$

For the discretisation of this set of equations, the derivatives are approximated using the backward Euler approximation for stability reasons [28]. The obtained set of difference equations can be rewritten and used for estimating the value of the ZARC element voltage in a discrete instant *k*, *vRC*[*k*], as a function of the current *i*[*k*] and the RC elements voltages *xRC*[*k*]:

$$\mathbf{x}\_{\rm R\mathcal{C}}[k] = \left(I\_{\rm R\mathcal{C}} - T A\_{\rm R\mathcal{C}}\right)^{-1} \mathbf{x}\_{\rm R\mathcal{C}}[k-1] + T \left(I\_{\rm R\mathcal{C}} - T A\_{\rm R\mathcal{C}}\right)^{-1} B\_{\rm R\mathcal{C}} i[k] \tag{13}$$

$$
\omega\_{RC}[k] = \mathcal{C}\_{RC} \mathbf{x}\_{RC}[k] + D\_{RC} i[k]\_{\prime} \tag{14}
$$

in which *InRC* is the *nRC* × *nRC* identity matrix, and *T* is the sampling period. Equations (13) and (14) can be used for implementing the response of the RC approximation given a current signal and a set of initial values for the RC branch voltages.

### *3.2. Approximation 2: Oustaloup Approach*

This approach relies on approximating the transfer function of the CPE in the ZARC element using a transfer function with integer order *nOU* (the same number of zeros and poles) in a given frequency range. One of the most popular methods for the approximation of the CPE transfer function with a rational transfer function of odd order *nOU* was presented by Oustaloup et al. [16]. The approximation is valid in the frequency range *ω* ∈ [*ωl*, *ωh*], where *ω<sup>l</sup>* and *ω<sup>h</sup>* are the lower and higher frequency limits, respectively. Following this method, the transfer function of the CPE can be rewritten as:

$$Z\_{\rm CPE}(\mathbf{s}) = \frac{1}{Q\_1 s^{\theta\_1}} = K\_{\rm CPE} \left(\frac{\mathbf{s}}{\omega\_{\rm c}}\right)^{-\phi\_1} \approx K\_{\rm CPE} \left(\frac{\omega\_{\rm l}}{\omega\_{\rm c}}\right)^{-\phi\_1} \prod\_{h=-N}^{N} \frac{1 + \frac{s}{\omega\_{\rm c(CPE)h}}}{1 + \frac{s}{\omega\_{\rm p(CPE)h}}} = Z\_{\rm OII(\rm CPE)}(\mathbf{s}),$$

where *ZOU*(*CPE*)(*s*) is the integer-order approximation of the CPE impedance using the OU approach; *<sup>ω</sup><sup>c</sup>* <sup>=</sup> <sup>√</sup>*ωlω<sup>h</sup>* is the central frequency between the bounds of the range of interest; and *KCPE* and *N* are given by *KCPE* = <sup>1</sup> *<sup>Q</sup>*1*ωφ*<sup>1</sup> *c* and *N* = *nOU*−<sup>1</sup> <sup>2</sup> . In order to estimate the *nOU* zeros *ωz*(*CPE*)*<sup>h</sup>* and *nOU* poles *ωp*(*CPE*)*<sup>h</sup>* in (15), Oustaloup et al. proposed [16]:

$$
\omega\_{z(CPE)h} = \omega\_l \left(\frac{\omega\_h}{\omega\_l}\right)^{\frac{h + \frac{n\_{\text{OII}} + \phi\_1}{2}}} \qquad \qquad \omega\_{p(CPE)h} = \omega\_l \left(\frac{\omega\_h}{\omega\_l}\right)^{\frac{h + \frac{n\_{\text{OII}} - \phi\_1}{2}}} \cdot \tag{15}
$$

The transfer function can be used to approximate the frequency response of the CPE. Then, the approximation *ZOU*(*s*) of the whole ZARC impedance using the OU approach is:

$$Z\_{ZARC}(\mathbf{s}) = \frac{R\_{p1}}{1 + R\_{p1}Q\_1\mathbf{s}^{\Phi\_1}} \approx \frac{R\_{p1}}{1 + \frac{R\_{p1}}{Z\_{OL(CPE)}\{\mathbf{s}\}}} = Z\_{OL}(\mathbf{s}).\tag{16}$$

In order to obtain a time representation of the response of this approximation, the transfer function *ZOU*(*s*) needs to be rewritten in zero-pole-gain form, as expressed in:

$$Z\_{\rm OII}(\mathbf{s}) = K\_{ZAR\mathbb{C}} \prod\_{h=1}^{n\_{\rm OI}} \frac{\left(\mathbf{s} - \omega\_{zh}\right)}{\left(\mathbf{s} - \omega\_{ph}\right)} \tag{17}$$

where *KZARC*, *ωzh* and *ωph* represent the gain, zeros and poles of the *ZOU*(*s*) transfer function, respectively. A time-domain state-space representation for the OU approximation with the structure presented in (17) was introduced in [28]:

$$
\dot{\mathbf{x}}\_{\rm OUT}(t) = A\_{\rm OUT} \mathbf{x}\_{\rm OUT}(t) + B\_{\rm OUT} i(t) \tag{18}
$$

$$w\_{OL}(t) = \mathcal{C}\_{OL}\mathbf{x}\_{OL}(t) + D\_{OL}i(t),\tag{19}$$

where *xOU* and *x*˙*OU* represent the system states and their derivatives, and *vOU* is the approximation of the ZARC element voltage. The matrices *AOU nOU* × *nOU*, *BOU nOU* × 1, *COU* 1 × *nOU* and *DOU* 1 × 1, are given by:

$$\begin{aligned} A\_{\text{OU}} &= \begin{bmatrix} \omega\_{p1} & 0 & 0 & \cdots & 0 \\ (\omega\_{p2} - \omega\_{z2}) & \omega\_{p2} & 0 & \cdots & 0 \\ (\omega\_{p3} - \omega\_{z3}) & (\omega\_{p3} - \omega\_{z3}) & \omega\_{p3} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ (\omega\_{p\text{n}\_{\text{OU}}} - \omega\_{z\text{n}\_{\text{OU}}}) & (\omega\_{p\text{n}\_{\text{OU}}} - \omega\_{z\text{n}\_{\text{OU}}}) & \cdots & (\omega\_{p\text{n}\_{\text{OU}}} - \omega\_{z\text{n}\_{\text{OU}}}) & \omega\_{p\text{n}\_{\text{OU}}} \end{bmatrix} \\ B\_{\text{OU}} &= K\_{ZARC} \begin{bmatrix} (\omega\_{p1} - \omega\_{z1}) & (\omega\_{p2} - \omega\_{z2}) & \cdots & (\omega\_{p\text{n}\_{\text{OU}}} - \omega\_{z\text{n}\_{\text{OU}}}) \end{bmatrix}^{\top} \\ C\_{\text{OU}} &= \begin{bmatrix} 1 & 1 & \cdots & 1 \end{bmatrix} \\ D\_{\text{OU}} &= K\_{ZARC} \end{aligned} \tag{21}$$

Again, the backward Euler approximation was used to discretise the obtained statespace system. This discrete representation is obtained by replacing *ARC*, *BRC*, *CRC* and *DRC* with *AOU*, *BOU*, *COU* and *DOU* in (13) and (14).

### *3.3. Approximation 3: Grünwald–Letnikov Approach*

The ZARC element response can be approximated by adopting a FO derivative definition in the time domain. The FO differential equations representing the ZARC voltage can be obtained from transfer function (3) by replacing *ZZARC*(*s*) = *V*(*s*)/*I*(*s*) and rewriting the equation as:

$$s^{\phi\_1}V(s) = \frac{I(s)}{Q\_1} - \frac{V(s)}{R\_{p1}Q\_1} \, \tag{24}$$

where *V*(*s*) and *I*(*s*) represent the Laplace transform of the ZARC voltage and current, respectively. Then, by applying inverse Laplace transform, the FO differential equation is rewritten as:

$$
\partial^{\phi\_1} v(t) = \frac{i(t)}{Q\_1} - \frac{v(t)}{R\_{p1}Q\_1}.\tag{25}
$$

D*φ*<sup>1</sup> represents the derivative of FO *φ*1. Among the multiple definitions of the FO derivative, the GL one is of particular interest, as it allows one to directly obtain difference equations for the approximations of the time response of a FO system [17]. The definition of the considered equation is:

$$\partial^{a}f(t) = \lim\_{T \to \infty} \frac{1}{T^{a}} \sum\_{h=0}^{\lfloor \frac{t}{T} \rfloor} (-1)^{h} \binom{a}{h} f(t - hT), \tag{26}$$

where the derivative of FO *α* of the causal function *f*(*t*) is computed between 0 and *t*. In (26), *T* is the sampling period, *<sup>t</sup> T* represents the integer part of *t*/*T* and ( *α <sup>h</sup>*) represents the Newtonian binomial coefficients generalised to real numbers, computed as:

$$
\binom{a}{h} = \frac{a(a-1)(a-2)\cdots(a-h+1)}{h!} = \frac{\Gamma(a+1)}{\Gamma(h+1)\Gamma(a-h+1)},\tag{27}
$$

where Γ(·) stands for the gamma function, which works as a generalisation of the factorial operator for real numbers.

It is worth noting that, according to the GL definition in (26), the derivative of order *α* of the function at time *t* depends on all the values of that function in [0, *t*], which is due to the non-local property of fractional derivatives [29].

By fixing the value of *T* to an appropriately low value for the application, and adopting the discrete variable *k* instead of the continuous time *t*, it is possible to obtain the first order discrete approximation of the FO derivative:

$$\partial^{\mathfrak{A}}f[k] = \frac{1}{T^{\mathfrak{a}}} \sum\_{h=0}^{k} (-1)^{h} \binom{\mathfrak{a}}{h} f[k-h]. \tag{28}$$

This approximation may be used to obtain a difference equation for the numerical evaluation of function *f* [*k*]. The discrete version also requires all the data points of *f* [*k*] since *k* = 0 for the computation of the derivative approximation, which may imply large memory requirements for simulations using this approach. This drawback may be addressed by applying the short-memory principle reported in reference [17], taking into account the behaviour of the signal in only the recent past, in the interval [*k* − *L*, *k*], where *L* is the memory length. Applying this principle, we rewrite:

$$\partial^{\mathfrak{a}}f[k] = \frac{1}{T^{\mathfrak{a}}} \sum\_{h=0}^{L} (-1)^{h} \binom{\mathfrak{a}}{h} f[k-h]. \tag{29}$$

This short-memory approximation allows one to implement numerical difference equations in cases in which the required memory is a critical constraint. Obviously, this introduces some inaccuracy, mostly manifested in the form of static error [30].

By replacing (29) in (25), we obtain the difference equation:

$$\upsilon\_{GL}[k] = T^{\Phi\_1} \frac{R\_{p1}}{R\_{p1}Q\_1 + T^{\Phi\_1}} \mathbf{i}[k] - \frac{R\_{p1}Q\_1}{R\_{p1}Q\_1 + T^{\Phi\_1}} \sum\_{h=1}^{L} (-1)^h \binom{\Phi\_1}{h} \upsilon\_{GL}[k-h],\tag{30}$$

which allows one to compute the approximation of the ZARC element voltage using the GL approach *vGL*. For the sake of comparison with the other implementation approaches, (30) can be written in matrix form:

$$
\sigma\_{GL}[k] = \mathbb{C}\_{GL} \mathbf{x}\_{GL}[k] + D\_{GL} i[k], \tag{31}
$$

where *xGL*[*k*] is a *L* × 1 vector with the *L* previous values of the ZARC voltage:

$$\mathbf{x}\_{GL} = \begin{bmatrix} \upsilon\_{GL}[k-1] & \upsilon\_{GL}[k-2] & \cdots & \upsilon\_{GL}[k-L] \end{bmatrix}^\top. \tag{32}$$

The matrices *CGL* and *DGL*, having size *L* × 1 and 1 × 1, respectively, are given by:

$$\mathcal{L}\_{GL} = -\frac{R\_{p1}Q\_1}{R\_{p1}Q\_1 + T\Phi\_1} \begin{bmatrix} -\binom{\phi\_1}{1} & \binom{\phi\_1}{2} & \cdots & (-1)^h \binom{\phi\_1}{h} & \cdots & (-1)^L \binom{\phi\_1}{L} \end{bmatrix} \tag{33}$$

$$D\_{\rm GL} = T^{\Phi\_1} \frac{R\_{p1}}{R\_{p1}Q\_1 + T^{\Phi\_1}}.\tag{34}$$

It is worth noting that for the matrix *CGL*, the coefficients (−1)*h*( *φ*1 *<sup>h</sup>* ) for the previous samples can be precomputed for the implementation of the ZARC voltage.

### **4. Accuracy Comparison**

Ideally, the data required for the validation of approximations of the responses of FO battery impedance models must consider EIS and pulsed current tests both performed under similar SoC, SoH and temperature conditions. Due to the lack of availability of such data in the known datasets, and in order to perform accuracy comparisons between the analysed approximations, a reference model based on the analytical solution of FO differential equations was utilised for the generation of the reference data.

### *4.1. Reference Data Generation from the Analytical Solution of Fractional Differential Equations*

Assuming that the current signal can be written as a set of steps:

$$\dot{a}(t) = \sum\_{h=1}^{N\_u} \mathcal{U}\_h \mu(t - t\_{uh}),\tag{35}$$

where *u*(*t*) corresponds to the unit step function. Each one of the *Nu* current steps in *i*(*t*) is characterised by an amplitude *Uh* and application time *tuh*.

The ZARC element voltage can be expressed as [17]:

$$w(t) = \sum\_{h=1}^{N\_u} \frac{\mathcal{U}\_h}{Q\_1} (t - t\_{\text{uh}})^{\phi\_1} E\_{\phi\_1, \phi\_1 + 1} \left( -\frac{1}{R\_{p1} Q\_1} (t - t\_{\text{uh}})^{\phi\_1} \right) u(t - t\_{\text{uh}}), \tag{36}$$

where the function *Eα*,*<sup>β</sup>* is the two-parameter Mittag–Leffler function, defined by a series expansion presented in Appendix B, which also includes a derivation of (36).

Equation (36) was implemented in Matlab® for the generation of the reference data. For the Mittag–Leffler function, the implementation introduced in [31] was employed. It is worth mentioning that such an analytical voltage representation is not suitable for online implementation, due to the limitations imposed by the assumed input signal and to the iterative nature of the Mittag–Leffler function computation, which makes the evaluation of a single data point highly demanding from a computational point of view.

### *4.2. Analysis of the Voltage Approximation Signals*

The evaluation of the accuracy of the considered approximations requires one to perform an analysis of the differences between the responses of the reference model and the approximation of interest, for given ZARC element parameters and while using the same input current signal, as illustrated in Figure 3.

**Figure 3.** A schematic diagram representing the accuracy tests performed.

The test proposed in Figure 3 requires one to define an input current signal, a set of ZARC elements parameters and the order or memory length of the approximations. Tests such as this one were performed for six sets of ZARC element parameters, for currents generated with different sampling and dynamic characteristics and with variations of *nRC*, *nOU* and *L*. In order to generate the input currents and select the ZARC elements to be employed, first, we focused on the typical middle-frequency range of the dynamic response of Li-ion batteries, namely, the range between 0.01 and 200 Hz. This dynamic range is normally associated with the response of the double layer capacitance and the charge transfer resistance [18].

The considered input signals contained two stages: one aimed at evaluating the transient responses of the different approximations, and the second stage was for testing the steady state error. An input signal sample is presented in Figure 4. In this signal, the first stage has a total duration of 200 s, for which the amplitude and duration of each current pulse were selected randomly in the ranges [−1, 1] A and [0.5, 10] s, respectively. The second stage contains one single 0.5 A step with a duration of 500 s, with fixed 150 s rests before and after the step.

**Figure 4.** An example of the current profiles used during the accuracy tests.

The six sets of ZARC parameters we used are presented in Table 1. The ZARC parameters were selected for obtaining characteristic frequencies *ω*<sup>0</sup> = (1/(*Rp*1*Q*1))1/*φ*<sup>1</sup> covering the frequency range of interest, with *ω*<sup>0</sup> = 2*π f*0, while keeping the parameters' values inside typical ranges—namely, *Rp*<sup>1</sup> <sup>∈</sup> [0.1, 100] <sup>m</sup>Ω, *<sup>Q</sup>*<sup>1</sup> <sup>∈</sup> [1, 1000] F s*φ*1−<sup>1</sup> and *φ*<sup>1</sup> ∈ [0.5, 0.9] [32,33].


**Table 1.** Values of the parameters for the ZARC elements considered during the accuracy tests.

The evaluation of accuracy was performed for:


Figures 5 and 6 present some examples of the voltage computed by mRC, OU and GL approximations for ZARC 4 employing *T* = 0.01 s. The plots in Figure 5a present the reference voltage and the voltages obtained for three mRC-based approximations during the first 20 s of the random stage of the accuracy test. In general, all the implementations were able to approximate the dynamics of the reference signal, with only appreciable differences for the approximation being of the lowest order. The error signals, presented in Figure 7a for this set of implementations, show spikes always under a few hundred μV during all the current steps, which reduce in magnitude as the order increases. Similar results can be observed for the OU-based implementations, as shown in Figure 5b and the first 200 s of the error signal in Figure 7b. The low-order OU approximations resulted in higher error magnitudes than the mRC ones. In the case of the GL approximations, as seen in Figure 5c, the lower *L* values caused higher offset errors. This is further illustrated by the error signal in Figure 7c during the first 200 s on which the spikes, at least for the lower memory lengths, seem to be wider than those in the mRC and OU approximations. An increase in *L* led to decreases in the magnitudes of the error signal, showing that adding terms to the sum in (30) leads to a better approximation of the analytical response, as expected.

**Figure 5.** Examples of ZARC 4 voltages during the dynamic stage of the accuracy test. (**a**) mRC; (**b**) OU; (**c**) GL.

**Figure 6.** Examples of ZARC 4 voltages during the static stage of the accuracy test. (**a**) mRC; (**b**) OU; (**c**) GL.

**Figure 7.** Examples of ZARC 4 voltage errors during the accuracy test. (**a**) mRC; (**b**) OU; (**c**) GL.

As is shown in the static state results presented in Figure 6a for three mRC-based implementations, increasing the order of the approximation reduced the errors during this stage. A higher value of *nRC* allowed a better approximation of the distribution of time constants represented by the FO element, leading to an extension of the validity of the approximation over a wider frequency range. Similar considerations apply for the OU approximations, presented in Figure 6b, except for the worse performance at low orders, under nine, with respect to the mRC case. The offset error for the GL approximations is more evident during the static stage of the test, as shown in Figures 6c and 7c. The static error obtained for GL implementations with *L* under a few thousand for ZARC 4 is considerably higher than the errors obtained for the other approximations. This highlights the main drawback of the GL approximations using the short memory principle: by reducing the number of previous samples that are used for the computation, some level of inaccuracy appears, particularly in static state. For the sake of completeness, it is worth mentioning that Podlubny [17] proposed a relationship for estimating a suitable memory length for the approximation of the FO derivative presented in (28), given an expected error level.

### *4.3. Effects of the Approximation Order and Memory Length on The Accuracy*

In order to analyse the effects of the parameters in each implementation approach, namely, *nRC*, *nOU* and *L*, on the approximation accuracy, the mean relative errors during the dynamic and the static stages were used as indicators. For ZARCs 1, 3 and 6, the mean relative errors during the dynamic and the static stages are presented in Figures 8a and 9a, respectively. The very similar results for the other ZARC elements are not reported for the sake of brevity.

**Figure 8.** Average relative errors during the dynamic stage of the test. (**a**) mRC; (**b**) OU; (**c**) GL.

**Figure 9.** Average relative errors during the static stage of the test. (**a**) mRC; (**b**) OU; (**c**) GL.

The error during the static stage of the tests for the RC approximations remained under 1%, even for the approximations with fewer RC branches. This result was expected, as in the fitting procedure adopted in this work (Appendix A), the sum of the resistances in the mRC approximation was set to match the value of *Rp*1, leading to similar voltage drops in the response after the initial transitory. On the other hand, even if the values of the relative error during the dynamic stage tended to decrease with the number of employed RC elements, considerable improvements were only obtained up to *nRC* = 9. After that point, further improvements could be achieved by decreasing the value of *T*. The observed behaviour at the highest orders may be associated with the variability introduced by the fitting procedure required for the computation of the parameters of the mRC case from the ZARC parameters. Nevertheless, for cases over *nRC* = 5, the mean relative errors were always under 5%.

In the case of the OU-based approximations, the mean relative errors during the dynamic and the static stages are presented in Figures 8b and 9b, respectively. In general, for the OU approximations during the dynamic stage, orders *nOU* higher than nine are required for reaching average relative errors under 5%. It is worth noting that, compared with the mRC approach, similar average relative errors were achievable in general, but with approximations of higher order. Again, the average relative errors for the static tests were almost always below 1%. For both stages of the accuracy tests, a monotonic reduction in the errors could be observed with increases in the approximation order, highlighting the advantage of computing the integer order approximation of the FO transfer function with a set of predefined equations instead of performing a fitting. This behaviour can be useful when trying to select the approximation order by analysing the accuracy–complexity trade off.

Then, for a set of GL-based approximations, with *L* values between 5 and 10,000 samples (which correspond to time windows between 0.05 s and 100 s), a similar accuracy analysis was performed. Figures 8c and 9c show the results obtained for the dynamic and static average relative errors obtained for this set of approximations. Regarding the results in the dynamic stage, only the higher memory lengths, over 500 samples, allowed us to reach values in the same order of magnitude as the ones obtained for the other approaches. In the case of the static stage, the average relative error is always about

one order for magnitude higher, and it is evident that higher memory lengths or sampling times need to be considered for reducing the static error to a similar range. The errors were higher for the ZARC elements with slower dynamics, showing that slower systems require longer memory lengths to reach an acceptable accuracy level.

### **5. Computational Burden Comparison**

A battery ECM can be used in an EMS as a part of the BMS state estimation structure or for battery simulation purposes during validation of energy management algorithms, particularly in real-time simulation scenarios. In both cases, considering that normally middle to low-end processing devices are often favoured due to budgetary restrictions, care needs to be taken in regard to the computational requirements of the battery model implementation. Here, we analyse those requirements in a general sense, by addressing the sizes of the matrices and the number of operations for each FOM implementation approach as indicators of required memory and computational complexity, respectively, in an eventual deployment.

For the three approaches, the time implementation relies on a set of matrix additions and multiplications. For mRC and OU-based approaches, the implementations consist of sets of equations in the form of (37) and (5) for the state and output equations, respectively. The discrete state equation established for the mRC approach can be generalised as:

$$\mathbf{x}[k] = A\_d \mathbf{x}[k-1] + B\_d i[k],\tag{37}$$

where *Ad* and *Bd*, namely, the discrete state and input matrices, are given by:

$$A\_d = \left(I\_{n\_{RC}} - TA\_{RC}\right)^{-1} \tag{38}$$

$$B\_d = T\left(I\_{\text{fl}\_{\text{RC}}} - TA\_{\text{RC}}\right)^{-1}B\_{\text{RC}}.\tag{39}$$

Conversely, the GL-based implementation relies only on a difference equation with the structure of (5), but in which *x*[*k*] does not represent the system states vector but a vector with *L* previous values of the ZARC voltage.

Table 2 presents the sizes of the matrices and vectors used in the three approaches. Even if the mRC and OU approaches seem to be equivalents in terms of memory requirements, the simpler structure of the mRC can be exploited for the reduction of its memory requirements. Additionally, even if the number of arrays required for a GL implementation is lower and its size dependency given *L* seems simpler than those for the matrices in the other approaches, it is worth keeping in mind that in general for a given accuracy level *L* will take values in the range of hundreds or thousands, while *nRC* or *nOU* will be under 20.


**Table 2.** Sizes of the matrices in the state and output equations in the analysed implementations.

Table 3 summarises the numbers of additions and multiplications required by each implementation approach. The specific structures of the matrices can be also exploited in the mRC and OU-based implementations, to refine the results presented in Table 3 by skipping the multiplications by zero and expressing the multiplication of column vectors by row vectors full of ones as an addition. The new operations count with this considerations is presented in Table 4.


**Table 3.** Numbers of additions and multiplications required for the analysed implementations.

**Table 4.** Numbers of additions and multiplications with simplifications.


Similarly to what was concluded for the array dimensions discussion, the expressions in Table 4 show that for the same order, a mRC implementation will require fewer operations than an OU one. It is worth mentioning that for the GL approach, even if there is not dependence on the square of *L* in the expressions for the number of operations required, the value of *L* needs to be considerably higher than the order for the other approaches to reach a given accuracy level.

The number of multiplications required for the evaluation of each type of implementation was used for assessing the computational burden in each case. Figure 10 plots the accuracy against the computational burden in terms of the number of multiplications for ZARC 4. The errors in static and dynamic stages are shown in Figure 10a,b, respectively.

**Figure 10.** Average relative errors vs. numbers of multiplications for mRC, OU and GL. (**a**) Dynamic stage; (**b**) static stage.

The curves for *T* = 0.01 s in Figure 10a show clearly how for a fixed mean relative error level in the dynamic stage of the accuracy tests, the number of required multiplications is always lower for the mRC approach, followed by the OU one. It is worth mentioning that the three approaches converge to values in the same order of magnitude for the mean relative error when increasing the complexity of the implementation; this may be an effect of the local truncation error due to the discretisation process. The asymptotic values of the analysed errors are comparable for a fixed sampling time. This can be ascribed to the fact that both the backward Euler and the GL derivative approximations are firstorder approximations, leading to an *O*(*T*) local truncation error, using big O notation [17]. To show the dependency on the sampling time of the identified error asymptotic values, the accuracy test was repeated using the current signal in Figure 4, but downsampled using *T* = 0.1 s. Those results are also summarised in Figure 10a, showing how the limit in the mean relative error is reached at a higher value, confirming the relationship between the sampling time and the maximum achievable accuracy.

The results presented in Figure 10b show that an increase in complexity has a more pronounced effect on the accuracy under static conditions, which tracks back to the requirement of higher orders or memory lengths for covering a wider time-constant range in the response approximation. In this case the previously observed oscillatory behaviour for higher orders in the mRC implementations is more evident for the two considered values of *T*. The effect of the local truncation error due to the discretisation process is also observed in this case, even if for lower values of the mean relative error.

From the results of the accuracy against complexity analysis, it can be clearly observed that the best accuracy–complexity compromise is offered by the mRC approach, allowing one to reach low error levels with a small computational burden. It is worth mentioning that this is the case provided good fitting of the ZARC element in the frequency domain can be performed. This is true for applications in which the model is used for battery simulation, but it is not the case when identification of impedance models from time measurements is required. In such instances, a good relationship between the parameters fitted from time-domain experiments with the frequency response is required. Thus, it is worth checking the suitability of the implementations for the time-domain identification of impedance parameters.

### **6. Analysis of the Suitability of Time-Domain Identification of the Battery Impedance**

One of the main reasons for adopting battery FOMs is the capability of accurately approximating the voltage response while requiring a low number of parameters, which is of interest for tasks such as state estimation and battery characterisation. For this reason, in order to illustrate the applicability of the considered FOM implementation approaches in the framework of battery model identification using time-domain measurements, a set of time-domain FOM fitting tests were performed.

For the fitting tests, the reference voltage data *v*[*k*] correspond to the random stage of the accuracy tests as the reference voltage signals presented in Figure 5. For all the ZARC elements and the implementations considered in the accuracy analysis, the associated parameters were fitted by minimising the mean square error between the reference voltage *v*[*k*] and *vx*[*k*], which corresponded to *vRC*[*k*], *vOU*[*k*] or *vGL*[*k*] depending on the evaluated approximation. For all cases, the minimisation problem was solved in Matlab® using the default particle swarm optimisation (PSO) algorithm, implemented by the Matlab® function "particleswarm". The default PSO algorithm employs a number of particles automatically selected as the minimum between 10 times the number of parameters to be fitted and 100 particles; a function tolerance of 10<sup>−</sup>6; and a maximum iteration number of 200 times the number of parameters to find [34].

The set of identified parameters changes depending on the approach evaluated. In the case of the GL-based implementations, the three parameters of the ZARC element, namely, *Rp*1, *Q*<sup>1</sup> and *φ*1, can directly be identified due to the nature of this implementation, where the time response of the FO element is directly approximated, as presented in Section 3.3. This can be observed in the schematic of the identification procedure presented in Figure 11c, where the parametrisation process takes the identified values of the ZARC parameters as inputs for generating the vectors required for the time-domain implementation. Similarly, for the OU approach, the direct identification of the ZARC element parameters from time measurements is possible due to the direct relationship between *Rp*1, *Q*1, *φ*<sup>1</sup> and the poles and zeros of the implemented integer-order transfer function, as introduced in Section 3.2. Figure 11b shows the implementation based on the OU approach, which was not modified for the fitting tests; only the source of the ZARC element parameters' changes, now being generated by the minimisation algorithm.

For the mRC approach, an initial frequency-domain fitting of the FO element impedance is required. For identification using time-domain measurements, this step cannot be performed. The resistance and capacitance values needed to be fit directly, as represented in the flowchart in Figure 11a.

On the one hand, for the OU and the GL approaches, the number of parameters to be identified was always 3, independently of the order *nOU* or memory length *L*. On the other hand, the number of parameters to identify with the mRC approach increased with the order, being equal to 2*nRC*. This highlights the main drawback of the mRC approach: when fitting the time response of the FO element, overfitting issues may arise due to the high number of parameters.

For all the fitting tests, the mean of the relative error between the reference voltage and the response of the fitted approximation was computed as an indicator of the goodness of the time-domain fit. Additionally, the indicator of how close the obtained impedance is to the expected one in the range from 0.01 to 20 Hz is the following:

$$\delta\_{\rm Z}(\omega) = \frac{\left| Z\_{\rm ZARC}(\omega) - Z\_{app}(\omega) \right|}{\left| Z\_{\rm ZARC}(\omega) \right|}. \tag{40}$$

Here, *δ<sup>Z</sup>* corresponds to the relative distance between the impedance of the fitted approximation, called *Zapp*, and the one of the original ZARC element, *ZZARC*.

The results obtained for ZARC 3 are summarised in Figure 12 as a plot of the relative error in time against the one in frequency. Even if relative errors in time under 5% were obtained for multiple instances of each approach, similar results in frequency were only reached for OU and GL approximations. In the case of mRC-based approximations, a good fitting in time is not necessarily translated to a low distance between the approximation impedance and the original one. It is worth mentioning that the minimisation algorithm was not optimised, but despite this, good accordance between the original and approximated impedance was reached for the OU and the GL approximations and not for the mRC case. Similar results were reached for all the other ZARC elements, illustrating how for the mRC approach, the lack of a direct relationship between the approximation parameters—namely, the resistance and capacitance values—and the ZARC impedance parameters reduces the possibility of identifying an approximation in time domain that also has a frequency response close to the real one.

An on-board-oriented implementation of these identification methods should address additional issues. One of them is the noise affecting the current and voltage measurements. Such noise contributes to causing a bias in the identified parameters, which should be compensated for with advanced fitting algorithms, as shown, for instance, in reference [35]. The noise compensation is worthy to be a matter of future study.

**Figure 12.** Fitting ZARC 3: average relative voltage error vs. average relative impedance distance.

### **7. Conclusions**

In this work, the three main approaches adopted in the literature for the implementation of the time-domain response of battery FOMs were introduced and compared in terms of accuracy, computational requirements and suitability for the time-domain identification of battery impedance. The study was performed in a simulation framework with six different ZARC elements, which are normally used for the approximation of battery impedance in the middle-frequency range. The comparison was performed in a simulation environment where the reference solution was an analytical expression for the response of a ZARC element under a multiple-step current. The proposed expression, obtained using FO calculus theory, was used for generating the reference data required for the accuracy analysis of the considered implementations. Even if the discussion focused on ZARC elements, the results can be extended to the Warburg element and to the total battery impedance.

The primary results of the study can be summarised as follows.


The selection of the FOM implementation method depends on the application requirements. On the one hand, if the interest is only in the battery response simulation, the mRC approach offers the best accuracy–complexity compromise, which is desirable for real-time simulations oriented towards the validation of energy management algorithms. On the other hand, if the application requires accurate identification of the impedance parameters from time-domain measurements, the OU approach offers the best compromise among identifying the impedance model parameters, the complexity and the accuracy requirements.

**Author Contributions:** Conceptualization, B.O.A. and W.Z.; methodology, B.O.A. and W.Z.; software, B.O.A.; validation, B.O.A.; formal analysis, B.O.A.; investigation, B.O.A. and W.Z.; resources, W.Z.; data curation, B.O.A.; writing—original draft preparation, B.O.A. and W.Z.; writing—review and editing, B.O.A., W.Z. and E.M.; visualization, B.O.A.; supervision, B.O.A. and W.Z.; project administration, W.Z. and E.M.; funding acquisition, W.Z. and E.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by funds from the the project "Holistic approach to EneRgy-efficient smart nanOGRIDS—HEROGRIDS" (PRIN 2017 2017WA5ZT3) within the Italian MUR 2017 PRIN programme; by the "Ex-WISCH e D-CODE" and FARB funds of the University of Salerno; and by funds of a public grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (reference: ANR-16-IDEX-0008).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

### **Abbreviations**

The following abbreviations are used in this manuscript:


### **Appendix A. Fitting of the Multiple RC Approaches' Parameters**

The values of *Rh* and *Ch*, from the mRC approximation of the ZARC impedance introduced in Equation (6), can be computed by solving the minimisation problem:

$$\min\_{\mathcal{Q}\_{\mathcal{KC}}} \left\{ \left| Z\_{ZAR\mathcal{C}}(j\omega) - Z\_{\mathcal{RC}}(j\omega) \right|^2 \right\},\tag{A1}$$

where the array of parameters is given by Θ*RC* = [*R*1, ··· , *RM*+1, *C*1, ··· , *CM*+1]. This optimisation problem has been previously solved in the literature. The most widely accepted solution is the one proposed by [15], which presents the approximations of ZARC elements for given values of *φ*<sup>1</sup> using 1, 3 and 5 RC branches. That approximation relies on a number of optimisation parameters precomputed in [15] only for a few discrete values of *φ*1. Here, we present a generalisation of this method for an odd number of RC branches *nRC* = *M* + 1 and any positive value of *φ*<sup>1</sup> < 1.

For the procedure, we take as starting point the parameters of the ZARC element, namely, *Rp*1, *Q*<sup>1</sup> and *φ*1. We want to approximate the ZARC frequency response with *M* + 1 RC branches, with the values of the resistances and capacitors sorted as follows: *RRC* = [*R*1, *R*2, ··· , *RM*/2+1, ··· , *RM*, *RM*+1] and *CRC* = [*C*1, *C*2, ··· , *CM*/2+1, ··· , *CM*, *CM*+1]. The generalisation can be obtained by defining the values of the resistance "in the middle", defined by the index *M*/2 + 1, according to:

$$R\_{M/2+1} = \frac{K\_1 R\_{p1} \sin\frac{\phi\_1 \pi}{2}}{1 + \cos\frac{\phi\_1 \pi}{2}}.\tag{A2}$$

For all the remaining resistances, except those with indexes 1 and *M* + 1, and by taking advantage of the the symmetry of the ZARC element impedance spectra, the values are computed as fractions of the difference between *Rp*<sup>1</sup> and *RM*/2<sup>+</sup>1, as expressed in:

$$R\_h = R\_{M+2-h} = \frac{K\_h \left(R\_{p1} - R\_{M/2+1}\right)}{2}, \quad 2 \le h \le \frac{M}{2}.\tag{A3}$$

Finally the most "external" resistances are computed as:

$$R\_1 = R\_{M+1} = \frac{R\_{p1} - \sum\_{h=2}^{M} R\_h}{2}.\tag{A4}$$

Then, the values of the capacitors can be estimated using the resistance and the characteristic frequency *ω*<sup>0</sup> = (1/(*Rp*1*Q*1))1/*φ*<sup>1</sup> values as follows:

$$\mathbb{C}\_{M/2+1} = \frac{1}{\omega \eta R\_{M/2+1}}\tag{A5}$$

$$\mathcal{C}\_{h} = \frac{1}{K\_{M/2+h}\omega\_{0}\mathcal{R}\_{h}}, \quad 1 \le h \le \frac{M}{2} \tag{A6}$$

$$\mathcal{C}\_{h} = \frac{1}{\frac{\omega\_{0}}{\overline{\mathcal{K}\_{3M/2+2-h}}} R\_{h}^{\prime}} \quad \frac{M}{2} + 1 \le h \le M + 1. \tag{A7}$$

While using this method for rewriting the optimisation problem presented in (A1), we can replace the parameters vector for Θ*RC* = [*K*1,..., *KM*], reducing the number of unknown parameters from 2*M* + 2 (the number of resistors plus the number of capacitors) to *M*. The solution of the optimisation problem presented in (A1), and the set of Equations (A2)–(A7), were implemented in Matlab®, allowing us to obtain the parameters of the RC circuit that approximate the frequency response of an ideal ZARC element.

### **Appendix B. ZARC Element—Analytical Voltage Expression**

For a ZARC element, the voltage *<sup>v</sup>*(*t*) = <sup>L</sup> <sup>−</sup>1{*V*(*s*)} can be computed as [17]:

$$v(t) = \mathcal{L}^{\ell-1} \left\{ \frac{\frac{1}{Q\_1}}{s^{\Phi\_1} + \frac{1}{R\_{p1}Q\_1}} I(s) \right\} = \frac{1}{Q\_1} \left[ t^{\Phi\_1 - 1} E\_{\Phi\_1 \Phi\_1} \left( -\frac{1}{R\_{p1}Q\_1} t^{\Phi\_1} \right) \right] \* i(t)$$

$$= \frac{1}{Q\_1} \int\_0^t (t - \tau)^{\Phi\_1 - 1} E\_{\Phi\_1 \Phi\_1} \left( -\frac{1}{R\_{p1}Q\_1} (t - \tau)^{\Phi\_1} \right) i(\tau) d\tau,\tag{A8}$$

where <sup>L</sup> <sup>−</sup>1{*X*(*s*)} is the inverse Laplace transform of *<sup>X</sup>*(*s*), *<sup>E</sup>α*,*<sup>β</sup>* is the two-parameter Mittag–Leffler function, the symbol ∗ stands for convolution and *τ* is the integration variable. The two-parameter Mittag–Leffler function is defined by the following series expansion:

$$E\_{a,\beta}(t) = \sum\_{h=0}^{\infty} \frac{t^h}{\Gamma(ah+\beta)},\tag{A9}$$

and it can be seen as a generalisation of the exponential function *e<sup>t</sup>* , which can be considered as a particular case of the *Eα*,*β*(*t*) function with *α* = *β* = 1 [17].

Equation (A8) can be used for computing the voltage response of the ZARC element from the current signal. The main limitation of this solution is that it requires an expression for the current in order to be evaluated. Additionally, this analytical expression is very demanding computationally, as the computation for a given time value requires all the previous values. Simplifications of this analytical expression can be achieved by considering a step current signal *i*(*t*) = *Uu*(*t* − *tu*), where *U* and *tu* are the step amplitude and time, and *u*(*t*) represents the unit step function. By taking (A8), and using *I*(*s*) for the Laplace transform of the current step, it is possible to obtain an expression that does not require the evaluation of an integral for the computation of the ZARC voltage:

$$v(t) = \mathcal{L}^{-1}\{V(s)\} = \mathcal{L}^{-1}\left\{\frac{\frac{1}{Q\_1}}{s^{\Phi\_1} + \frac{1}{R\_{p1}Q\_1}}\frac{\mathcal{U}}{s}s^{-t\_u s}\right\}
$$

$$= \frac{\mathcal{U}}{Q\_1}(t - t\_u)^{\Phi\_1} E\_{\Phi\_1\Phi\_1 + 1}\left(-\frac{1}{R\_{p1}Q\_1}(t - t\_u)^{\Phi\_1}\right)u(t - t\_u). \tag{A10}$$

Using the superposition principle, it is possible to extend these results for current signals composed by combinations of sets of *Nu* step signals, meaning that they can be written as:

$$\dot{u}(t) = \sum\_{h=1}^{N\_w} \mathcal{U}\_h u(t - t\_{uh})\_\prime \tag{A11}$$

which is true for inputs such as square or pseudo-random binary sequence signals. If such an input current signal is considered, the resulting voltage signal for the ZARC element is given by:

$$v(t) = \sum\_{h=1}^{N\_{\rm u}} \frac{lI\_h}{Q\_1} (t - t\_{\rm uh})^{\Phi\_1} E\_{\Phi\_1, \Phi\_1 + 1} \left( -\frac{1}{R\_{p1} Q\_1} (t - t\_{\rm uh})^{\Phi\_1} \right) u(t - t\_{\rm uh}). \tag{A12}$$

### **References**

