Solutions to the Sub-Optimality and Stability Issues of Recursive Pole and Zero Distribution Algorithms for the Approximation of Fractional Order Models

Abstract

This paper analyses algorithms currently found in the literature for the approximation of fractional order models and based on recursive pole and zero distributions. The analysis focuses on the sub-optimality of the approximations obtained and stability issues that may appear after approximation depending on the pole location of the initial fractional order model. Solutions are proposed to reduce this sub-optimality and to avoid stability issues.

1. Introduction

Fractional order models are generalizations of models described by differential equations or state space descriptions in which classical derivatives are replaced by fractional derivatives [1,2,3]. These models are now widely used to characterize the behavior of many systems in diverse areas. For example:

-

electrochemistry in which charge diffusion in batteries can be described by Randles’ models [4,5], or other kinds of fractional models [6];

-

thermal conduction where the exact solution of the heat equation in a semi-infinite medium links the heat rate to the surface temperature by a fractional differentiation of order 0.5 and applied to the thermal modeling of buildings [7,8];

-

biology for modelling complex dynamics in biological tissues [9];

-

mechanics with the dynamical property of viscoelastic materials and for wave propagation problems in these materials [10];

-

acoustics where fractional differentiation is used to model visco-thermal losses in wind instruments [11];

-

robotics through environmental modeling [12];

-

electrical distribution networks [13];

-

modelling of explosive materials [14].

The main reason for the use of fractional order models is their ability to generate long memory behaviors (time power law relaxations) in the same way as the systems previously mentioned. However, this interesting property comes at the price of several defects that have important consequences in the field of dynamical system analysis:

-

the physical consistency of the state concept is questionable, making it necessary to introduce the notion of a pseudo-state [15,16];

-

they are defined on an infinite time interval or on an infinite space domain (depending on the representation used) [15,16];

-

as a result, they are adapted to studying the input-output behavior of a system, but not its internal properties [17] (initialization [18], internal stability, exact controllability [19], exact observability, exact flatness);

-

the implementation of these models requires an approximation step [8,20].

This work looks at this approximation step. There are several methods in the literature that allow such an approximation. Among all these methods, a very well-known one because of its efficiency and the simple algorithm that it implements is based on the approximation of the fractional order differentiation or integration operator using a recursive (geometric) pole and zero distribution [21,22]. In this paper, this method is analyzed in terms of sub-optimality and a simple solution is proposed to improve the approximation accuracy. Some stability issues resulting from the approximation of a fractional model are also discussed. The results presented in this paper will contribute to obtaining more accurate and stable approximations of a fractional model, and above all will help the reader to understand that the geometric distribution of poles and zeros (also called “recursivity” in the literature) for the approximation of a fractional order integrator and differentiator is one among an infinity of other permitted distributions and cannot be interpreted as the physical reason for the observed long memory behaviors.

2. The Existing Algorithms Based on Pole and Zero Recursive (Geometric) Distribution

This section reviews the algorithms based on the recursive (geometric) pole and zero distribution found in the literature for the synthesis of a fractional order integrator or differentiator respectively described by the transfer functions:

$$I^{\nu }\left(s\right)=\frac{1}{s^{\nu }}\mathrm{and}{D}^{\nu }\left(s\right)=s^{\nu }\mathrm{with}0<\nu <1.$$

(1)

2.1. Approximation of a Fractional Integrator and Differentiator by a Recursive (Geometric) Distribution of Pole and Zeros: A Graphical Approach

The demonstration of the approximation of a fractional order integrator operator frequency response is here done graphically. This method appeared for the very first time in the literature in the 1960s in two studies, apparently carried out independently by the respective authors [23,24]. Some years later, this demonstration was used for the analog implementation of fractional integrators [25]. Subsequent work by other authors also contributed to this synthesis method [26,27,28].

For real orders, to obtain an approximation of a fractional integration operator, a solution consists in limiting the frequency band on which the fractional behavior is required. This first leads to the approximation:

$$I^{\nu }\left(s\right)\approx I_a^{\nu }\left(s\right)=C_0{\left(\frac{1+\frac{s}{{\omega }_h}}{1+\frac{s}{{\omega }_b}}\right)}^{\nu }\mathrm{with}{C}_0={\left(\frac{{\omega }_h}{{\omega }_b}\right)}^{\nu }{\left(\frac{{\omega }_b{}^2+1}{{\omega }_h{}^2+1}\right)}^{\frac{\nu }{2}}.$$

(2)

As shown by Figure 1, the Bode plots of relation (2) exhibit:

-

a gain whose slope is a fractional multiple of −20 dB/decade,

-

a constant phase whose value is a fractional multiple of −90°.

These Bode plots can be approximated using a distribution of poles and zeros. This leads to the following algorithm (Algorithm 1).

Algorithm 1. Fractional integrator approximation—first method

1. Initialize $$r=\alpha \eta ={\left(\frac{{\omega }_h}{{\omega }_b}\right)}^{\frac{1}{N}}\alpha =r^{\nu }\eta =r^{1-\nu }{\omega }_1={\eta }^{1/2}{\omega }_b{{\omega }^{\prime }}_1=\alpha {\omega }_1$$ (3) 2. For $i\in \left[1\dots N\right]$ do $${{\omega }^{\prime }}_{i+1}=r\omega {{}^{\prime }}_i{\omega }_{i+1}=r{\omega }_i$$ (4)    End for 3. Compute $$C{{}^{\prime }}_0=\frac{{\displaystyle \prod _{k=1}^N{\left(1+{\left(\frac{1}{{\omega }_k}\right)}^2\right)}^{\frac{1}{2}}}}{{\displaystyle \prod _{k=1}^N{\left(1+{\left(\frac{1}{\omega {{}^{\prime }}_k}\right)}^2\right)}^{\frac{1}{2}}}}$$ (5) 4. Define the fractional integrator (1) approximation in the frequency band [ωb, ωh], by the transfer function $$I_a^{\nu }\left(s\right)=C_0{\left(\frac{1+\frac{s}{{\omega }_h}}{1+\frac{s}{{\omega }_b}}\right)}^{\nu }\approx I_N^{\nu }\left(s\right)=C{{}^{\prime }}_0\frac{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{\omega {{}^{\prime }}_k}\right)}}{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{{\omega }_k}\right)}}.$$ (6)

As shown in [26], with this algorithm, the corner frequencies ${\omega }_k$ and $\omega {{}^{\prime }}_k$ (respectively the poles and zeros of the transfer function $I_N^{\nu }\left(s\right)$) are geometrically distributed to obtain the required frequency behavior.

Figure 1 highlights this distribution and also compares the asymptotic Bode plots of $I^{\nu }\left(s\right)$, $I_a^{\nu }\left(s\right)$ and $I_N^{\nu }\left(s\right)$. It also shows that the high and low asymptotic frequency behaviors are constant. In [22], to obtain an integer integral like asymptotic behavior at low and high frequency, relation (2) is replaced by:

$$I^{\nu }\left(s\right)\approx I_a^{\nu }\left(s\right)=C_0\frac{1}{s}{\left(\frac{1+\frac{s}{{\omega }_h}}{1+\frac{s}{{\omega }_b}}\right)}^{\nu -1}\mathrm{with}{C}_0={\left(\frac{{\omega }_h}{{\omega }_b}\right)}^{\nu -1}{\left(\frac{{\omega }_b{}^2+1}{{\omega }_h{}^2+1}\right)}^{\frac{\nu -1}{2}}.$$

(7)

This leads to the following algorithm (Algorithm 2).

Algorithm 2. Fractional integrator approximation—second method

1. Initialize $$r=\alpha \eta ={\left(\frac{{\omega }_h}{{\omega }_b}\right)}^{\frac{1}{N}}\alpha =r^{1-\nu }\eta =r^{\nu }\omega {{}^{\prime }}_1={\eta }^{1/2}{\omega }_b{\omega }_1=\alpha \omega {{}^{\prime }}_1$$ (8) 2. For $i\in \left[1\dots N\right]$ do $$\omega {{}^{\prime }}_{i+1}=r\omega {{}^{\prime }}_i{\omega }_{i+1}=r{\omega }_i$$ (9)    End for 3. Compute $C{{}^{\prime }}_0$ with relation (5) 4. Define the fractional integrator (1) approximation in the frequency band [ωb, ωh], by the transfer function $$I_a^{\nu }\left(s\right)=C_0\frac{1}{s}{\left(\frac{1+\frac{s}{{\omega }_h}}{1+\frac{s}{{\omega }_b}}\right)}^{\nu -1}\approx I_N^{\nu }\left(s\right)=C{{}^{\prime }}_0\frac{1}{s}\frac{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{\omega {{}^{\prime }}_k}\right)}}{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{{\omega }_k}\right)}}.$$ (10)

For a fractional differentiation operator defined in the Laplace domain by

$$D^{\nu }\left(s\right)=s^{\nu }\mathrm{with}0<\nu <1$$

(11)

the frequency truncation leads to

$$D^{\nu }\left(s\right)\approx D_a^{\nu }\left(s\right)=C_0\left(\frac{1+\frac{s}{{\omega }_b}}{1+\frac{s}{{\omega }_h}}\right)\mathrm{with}{C}_0={\left(\frac{{\omega }_b}{{\omega }_h}\right)}^{\nu }{\left(\frac{{\omega }_h{}^2+1}{{\omega }_b{}^2+1}\right)}^{\frac{\nu }{2}}$$

(12)

and permits the following algorithm (Algorithm 3).

Algorithm 3. Fractional differentiator approximation—first method

1. Initialize $$r=\alpha \eta ={\left(\frac{{\omega }_h}{{\omega }_b}\right)}^{\frac{1}{N}}\alpha =r^{\nu }\eta =r^{1-\nu }\omega {{}^{\prime }}_1={\eta }^{1/2}{\omega }_b{\omega }_1=\alpha \omega {{}^{\prime }}_1$$ (13) 2. For $i\in \left[1\dots N\right]$ do $$\omega {{}^{\prime }}_{i+1}=r\omega {{}^{\prime }}_i{\omega }_{i+1}=r{\omega }_i$$ (14)    End for 3. Compute $C{{}^{\prime }}_0$ with relation (5) 4. Define the fractional differentiator (12) approximation in the frequency band [ωb, ωh], by the transfer function $$D_a^{\nu }\left(s\right)=C_0{\left(\frac{1+\frac{s}{{\omega }_b}}{1+\frac{s}{{\omega }_h}}\right)}^{\nu }\approx D_N^{\nu }\left(s\right)=C{{}^{\prime }}_0\frac{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{\omega {{}^{\prime }}_k}\right)}}{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{{\omega }_k}\right)}}.$$ (15)

By analogy with Algorithm 2, the following algorithm (Algorithm 4) permits to get an integer differentiation like asymptotic behavior in low and high frequency.

Algorithm 4. Fractional differentiator approximation—second method

1. Initialize $$r=\alpha \eta ={\left(\frac{{\omega }_h}{{\omega }_b}\right)}^{\frac{1}{N}}\alpha =r^{1-\nu }\eta =r^{\nu }{\omega }_1={\eta }^{1/2}{\omega }_b\omega {{}^{\prime }}_1=\alpha {\omega }_1$$ (16) 2. For $i\in \left[1\dots N\right]$ do $$\omega {{}^{\prime }}_{i+1}=r\omega {{}^{\prime }}_i{\omega }_{i+1}=r{\omega }_i$$ (17)    End for 3. Compute $C{{}^{\prime }}_0$ with relation (5) 4. Define the fractional differentiator (12) approximation in the frequency band [ωb, ωh], by the transfer function $$D_a^{\nu }\left(s\right)=C_{0}s{\left(\frac{1+\frac{s}{{\omega }_b}}{1+\frac{s}{{\omega }_h}}\right)}^{\nu -1}\approx D_N^{\nu }\left(s\right)=C{{}^{\prime }}_{0}s\frac{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{\omega {{}^{\prime }}_k}\right)}}{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{{\omega }_k}\right)}}.$$ (18)

The set of algorithms previously established for a real fractional differentiation order can be extended to a complex fractional differentiation order ν = a + ib [21].

2.2. Approximation of a Fractional Integrator by a Recursive Distribution of Poles and Zeros: An Analytical Approach

For the approximation given in relation (1) using Algorithm 1, the analytical demonstration of the graphical approach presented in the previous section was done in [21] by considering the limit case where ω_b tends towards 0 and ω_h tends towards infinity. However, a simpler demonstration based on the impulse response of a fractional operator is now used. This impulse response is defined by:

$$I^{\nu }\left(t\right)={\mathcal{L}}^{-1}\left\{\frac{1}{s^{\nu }}\right\}$$

(19)

where ${\mathcal{L}}^{-1}$ denotes the inverse Laplace–Transform. The residue method for the computation of the inverse Laplace–Transform leads to [3]:

$$I^{\nu }\left(t\right)=\frac{\sin \left(\nu \pi \right)}{\pi }{\displaystyle \underset{0}{\overset{\infty }{\int }}\frac{e^{-tx}}{x^{\nu }}}dx0<\nu <1.$$

(20)

The Laplace–Transform of (19) is thus given by:

$$I^{\nu }\left(s\right)=\frac{\sin \left(\nu \pi \right)}{\pi }{\displaystyle \underset{0}{\overset{\infty }{\int }}\frac{1}{x^{\nu }\left(s+x\right)}}dx.$$

(21)

Using the change of variable $x=e^z$ and thus $dx=e^{z}dz$, then:

$$I^{\nu }\left(s\right)=\frac{\sin \left(\nu \pi \right)}{\pi }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\frac{e^z}{e^{\nu z}\left(s+e^z\right)}}dz=\frac{\sin \left(\nu \pi \right)}{\pi }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\frac{e^{-\nu z}}{\left(\frac{s}{e^z}+1\right)}}dz.$$

(22)

Discretization of integral (22) leads to:

$$I^{\nu }\left(s\right)=\frac{\sin \left(\nu \pi \right)}{\pi }{\displaystyle \sum _{k=-\infty }^{\infty }\frac{e^{-\nu k\Delta z}}{\frac{s}{e^{k\Delta z}}+1}\Delta z}=\frac{\sin \left(\nu \pi \right)}{\pi }{\displaystyle \sum _{k=-\infty }^{\infty }\frac{e^{\left(1-\nu \right)k\Delta z}}{s+e^{k\Delta z}}}\Delta z,$$

(23)

in which $\Delta z$ denotes the sampling interval. Relation (22) also highlights the link between a fractional integrator impulse response and the Prony series as the inverse Laplace transform of (23) leads to:

$${\mathcal{L}}^{-1}\left\{I^{\nu }\left(s\right)\right\}=I^{\nu }\left(t\right)=\frac{\sin \left(\nu \pi \right)}{\pi }{\displaystyle \sum _{k=-\infty }^{\infty }{A}_k{e}^{-{\alpha }_{k}t}}$$

(24)

with ${\alpha }_k=e^{k\Delta z}$ and $A_k=e^{\left(1-\nu \right)k\Delta z}\Delta z$.

Relation (24) also highlights that the poles α_k (always greater than 0) are linked by the ratio α_k+1/α_k = e^−Δz. Regarding the zeros, to the best of the author’s knowledge there is no demonstration of a link with the same ratio as that imposed by Algorithm 1. However, by applying a partial fraction decomposition to relation (3) it can be written as:

$$I_N^{\nu }\left(s\right)=C_0\frac{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{\omega {{}^{\prime }}_k}\right)}}{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{{\omega }_k}\right)}}={\displaystyle \sum _{k=1}^N\frac{A{{}^{\prime }}_k}{s+\alpha {{}^{\prime }}_k}}\mathrm{with}\alpha {{}^{\prime }}_k={\omega }_k\mathrm{and}A{{}^{\prime }}_k=C_0\frac{{\displaystyle \prod _{i=1}^N{\omega }_k\left(1+\frac{{\omega }_k}{\omega {{}^{\prime }}_i}\right)}}{{\displaystyle \prod _{\begin{array}{l}i=1\\ i\ne k\end{array}}^N\left(1+\frac{{\omega }_k}{{\omega }_i}\right)}}.$$

(25)

By choosing α′_k = α_k, it can be shown using a graphical representation (see Figure 2) that coefficient A′_k in relation (25) tends towards coefficient A_k in (24), under the same hypothesis as the one used in [21]:

-

a number N of poles that tends towards infinity,

-

a ratio r or Δz that tends towards 1.

3. Sub-Optimality of Algorithms 1–4 and Beyond Geometric Distribution

The algorithms 1 to 4 described in Section 2 have been used in many applications, due to the simplicity of their implementation. However, it can be demonstrated by an example that Algorithm 1, for instance, is sub-optimal (a similar demonstration can be made for Algorithms 2 to 4). Let us assume that the aim is to compute the approximation of transfer function (1) on the frequency band [ω_b, ω_h] = [1, 10⁶], with ν = 0.3. Using Algorithm 1 with N = 10 poles and zeros, the following factors are obtained:

$$r=\alpha \eta =3.9811\eta =2.6303$$

(26)

which enable all the corner frequencies of relation (3) to be computed using relations (4) and (5). The Bode plots of relations (1) to (3) are represented by Figure 3. To evaluate the efficiency of Algorithm 1, the relative error and absolute error respectively computed by the relations

$$E_{rel}(\omega )={\left|\frac{I_a^{\nu }\left(s\right)}{I_N^{\nu }\left(s\right)}-1\right|}_{s=j\omega }{E}_{abs}(\omega )={\left|I_a^{\nu }\left(s\right)-I_N^{\nu }\left(s\right)\right|}_{s=j\omega }$$

(27)

are represented by Figure 4. The maximum values of $E_{rel}(\omega )$ and $E_{abs}(\omega )$ are respectively given by:

$$\underset{\omega }{\max }\left|E_{rel}(\omega )\right|=1.2762\times {10}^{-2}\underset{\omega }{\max }\left|E_{abs}(\omega )\right|=1.1382\times {10}^{-2}.$$

(28)

Then, an optimization program aiming at reducing the relative error $E_{rel}(s)$ by looking for the best location of the N poles and zeros in the transfer function $I_N^{\nu }\left(s\right)$ is implemented under a Matlab environment. The optimization program is initialized with the pole and zero values obtained with Algorithm 1. The resulting relative error and absolute error are shown in Figure 5. The maximum values of $E_{rel}(\omega )$ and $E_{abs}(\omega )$ are now respectively given by:

$$\underset{\omega }{\max }\left|E_{rel}(\omega )\right|=2.8689\times {10}^{-3}\underset{\omega }{\max }\left|E_{abs}(\omega )\right|=2.5622\times {10}^{-3}$$

(29)

and are almost five times smaller than the ones given by Algorithm 1.

To improve the accuracy of Algorithms 1–4, it is necessary to go beyond the geometric distribution of poles and zeros, but without introducing great complexity. With N poles $[{\omega }_0,\dots ,{\omega }_{N-1}]$, to define these generalizations, let us define the ratio r such that:

$$\frac{{\omega }_{N-1}}{{\omega }_0}=r^N.$$

(30)

For a geometric distribution, the following relation holds:

$$\log \left({\omega }_k\right)=\log \left({\omega }_b\right)+\phi \log \left(r\right)\phi \in \left[0\dots N\right].$$

(31)

It can be generalized, among many others, by the following distributions:

$$D_1:\log \left({\omega }_k\right)=\log \left({\omega }_b\right)+\frac{{\phi }^2}{N}\log \left(r\right)\phi \in \left[0\dots N\right],$$

(32)

$$D_2:\log \left({\omega }_k\right)=\log \left({\omega }_b\right)+\sqrt{N}\sqrt{\phi }\log \left(r\right)\phi \in \left[0\dots N\right].$$

(33)

In Figure 6, these two distributions are compared with the geometric one. It shows that distribution D₁ makes it possible to bring the first poles closer together while distribution D₂ makes it possible to increase the distance between them. However, no additional tuning parameters are associated to these distributions. It is thus proposed to replace the nonlinear function of N in relations (32) and (33) by a polynomial. The following distribution is thus obtained:

$$D_3:\log \left({\omega }_k\right)=\log \left({\omega }_b\right)+P\left(\phi \right)\log \left(r\right)\mathrm{with}P\left(\phi \right)={\displaystyle \sum _{m=0}^M{a}_m{\phi }^m}.$$

(34)

With M = 2, to ensure that corner frequencies really belong in $[{\omega }_b,\dots ,{\omega }_h]$, the following constraints must be imposed: P(0) = 0 and P(N) = N. This leads to imposing a₀ = 0 and a₁N + a₂N² = N, thus a₂ = (1 − a₁)/N.

With M = 3, to ensure that corner frequencies really belong in $[{\omega }_b,\dots ,{\omega }_h]$, the same constraints must be imposed leading to a₀ = 0 and $a_1+a_{2}N+a_3{N}^2=1$. Moreover, if it is imposed that distribution D₃ fits the geometric distribution in the middle of the interval $[{\omega }_b,\dots ,{\omega }_h]$, the following constraint must be added: P(N/2) = 1/2 and thus $\frac{a_1}{2}+\frac{a_{2}N}{4}+\frac{a_3{N}^2}{8}=1/2$. Finally, parameters a₂ and a₃, meet the following system of equations:

$$\left[\begin{array}{cc}N& N^2\\ \frac{N}{4}& \frac{N^2}{8}\end{array}\right]\left[\begin{array}{c}{a}_2\\ a_3\end{array}\right]=\left[\begin{array}{c}1-a_1\\ \frac{1-a_1}{2}\end{array}\right].$$

(35)

By varying parameter a₁, Figure 7 shows the large number of distributions that can be obtained in both cases.

In the two cases, it must be pointed out that only one parameter needs to be tuned to obtain a wide variety of distributions.

For the geometric distribution, according to Algorithm 1, the poles, and zeros in relation (3) are defined by relation (34) with

$$\phi =k-1+\frac{1-\nu }{2},\mathrm{for}\mathrm{pole}{\omega }_k$$

(36)

and

$$\phi =k-1+\frac{1+\nu }{2}\mathrm{for}\mathrm{zero}\omega {{}^{\prime }}_k.$$

(37)

To define the poles and the zeros using the distribution introduced, the same values of φ are used, but the values of ${\omega }_k$ and $\omega {{}^{\prime }}_k$ are given by relation (34) leading to the following algorithm (Algorithm 5).

Algorithm 5. Fractional integrator approximation—improved method

1. Initialize $$P\left(\phi \right)={\displaystyle \sum _{m=0}^M{a}_m{\phi }^m}r={\left(\frac{{\omega }_h}{{\omega }_b}\right)}^{\frac{1}{N}}$$ (38) 2. For $i\in \left[1\dots N\right]$ do $$\log \left({\omega }_k\right)=\log \left({\omega }_b\right)+P\left(\phi \right)\log \left(r\right)\mathrm{with}\phi =k-1+\frac{1-\nu }{2}$$ (39) $$\log \left(\omega {{}^{\prime }}_k\right)=\log \left({\omega }_b\right)+P\left(\phi \right)\log \left(r\right)\mathrm{with}\phi =k-1+\frac{1+\nu }{2}$$ (40)    End for 3. Compute $C{{}^{\prime }}_0$ with relation (5) 4. Define the fractional integrator (1) approximation in the frequency band [ωb, ωh], by the transfer function $$I_a^{\nu }\left(s\right)=C_0{\left(\frac{1+\frac{s}{{\omega }_h}}{1+\frac{s}{{\omega }_b}}\right)}^{\nu }\approx I_N^{\nu }\left(s\right)=C{{}^{\prime }}_0\frac{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{\omega {{}^{\prime }}_k}\right)}}{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{{\omega }_k}\right)}}.$$ (41)

The efficiency of this algorithm is highlighted on the example used at the beginning of Section 3, but with a distribution D₃ and M = 3. The parameter a₁ is computed to reduce the relative error given in relation (27). Given that only one parameter has to be found, it is done by a gridding of parameter a₁ on the interval [−0.1, 3]. A comparison of relative and absolute errors obtained with Algorithms 1 and 5 is given in Figure 8. The maximum values of these absolute and relative errors are given by:

$$\underset{\omega }{\max }\left|E_{rel}(\omega )\right|=1.0133\times {10}^{-2}\underset{\omega }{\max }\left|E_{abs}(\omega )\right|=3.1903\times {10}^{-3}.$$

(42)

In comparison with Algorithm 1, these values are greatly reduced with Algorithm 5.

4. Fractional Model Approximation and Stability Issues

We now consider a fractional order model described by the pseudo state space description,

$$\{\begin{array}{c}\frac{d^{\nu }}{d{t}^{\nu }}x\left(t\right)=Ax\left(t\right)+Bu\left(t\right)\\ y\left(t\right)=Cx\left(t\right)+Du\left(t\right)\end{array},$$

(43)

where $x\left(t\right)\in {\mathsf{R}}^n$ is the pseudo state vector, $\nu$ is the fractional order of the system and $A\in {\mathsf{R}}^{nxn}$, $B\in {\mathsf{R}}^{nxm}$, $C\in {\mathsf{R}}^{pxn}$ and $D\in {\mathsf{R}}^{pxm}$ are constant matrices. A solution for the implementation of such a model consists in the approximation of the fractional order derivative by the approximation given by Algorithm 3 or Algorithm 4. However, such an approximation has an impact on model stability. It is well known that model (43) is stable if and only if all the eigenvalues of matrix A belong to the stability domain D_s defined by:

$${\mathcal{D}}_s=\left\{s=\rho e^{i\theta }\in \mathsf{C}/\left|\theta \right|>\frac{\pi }{2}\right\}.$$

(44)

To explain how the stability domain is transformed by approximation (15) in Algorithm 3 or (18) in Algorithm 4, it is assumed that all the eigenvalues of model (43) are distinct. Thus, there exists a change of variable z(t) = P x(t), such that model (43) can be written as:

$$\{\begin{array}{c}\frac{d^{\nu }}{d{t}^{\nu }}z\left(t\right)=A_{P}z\left(t\right)+B_{P}u\left(t\right)\\ y\left(t\right)=C_{P}z\left(t\right)+D_{P}u\left(t\right)\end{array},$$

(45)

where matrix A_p is diagonal and can be written $A_P=diag({\lambda }_i)$, $i\in \left[1\dots n\right]$. Using the Laplace transform, the characteristic equations of model (45) are thus

$$s^n-{\lambda }_i=0,i\in [1\dots n].$$

(46)

with approximation (15) or (18), these characteristic equations become:

$$C_0\frac{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{\omega {{}^{\prime }}_k}\right)}}{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{{\omega }_k}\right)}}-{\lambda }_i=0\mathrm{and}{C}_{0}s\frac{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{\omega {{}^{\prime }}_k}\right)}}{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{{\omega }_k}\right)}}-{\lambda }_i=0.$$

(47)

After approximation, model (43) is stable if characteristic Equation (47) has no root in the right half part of the complex plane. To evaluate the location of these roots, the Cauchy argument principle can be used. For that, the path Г in Figure 9 can be used.

Images of path Г by approximations (15) and (18) are represented by Figure 10. They are compared with the image of this path by s^ν.

From these curves, the stability domain of model (43) can be deduced. The image of path Г by relation (47) is the image of path Г by relations (15) and (18) translated by the vector ${\overrightarrow{V}}_{{\lambda }_i}$ that appears in Figure 10 and that links the point of the coordinate (Re(λ_i), Im(λ_i)) to the origin of the complex plane. Thus, after this translation the origin of the complex plane is located at the point (Re(λ_i), Im(λ_i)). As the denominator of relations (47) has no root inside the path Г, if the origin (after translation) of the complex plane is inside one of the images, then according to argument principle, the corresponding characteristic equation has one root inside the path Г and the system (43) is unstable. This permits the following theorem.

Theorem 1.

The stability domain in the complex plane, for the eigenvalues of matrix A_P (or for the poles) of model (45) or equivalently for model (43) is

-

outside the domain defined by,$C_0\frac{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{\omega {{}^{\prime }}_k}\right)}}{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{{\omega }_k}\right)}}$,$s=j\omega$,$\omega \in \left[-\infty ,\infty \right]$for approximation (15)

-

outside the domain on the right of the curve$C_{0}s\frac{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{\omega {{}^{\prime }}_k}\right)}}{{\displaystyle \prod _{k=1}^N\left(1+\frac{s}{{\omega }_k}\right)}}$,$s=j\omega$,$\omega \in \left[-\infty ,\infty \right]$, for approximation (18).

These stability domains are illustrated by Figure 11.

This analysis also highlights stability issues that can occur after approximation. Considering Figure 12, which is an enlargement of Figure 11, it can be seen that some instability domains for model (43) after approximation are inside the stability domains for model (43) (before approximation). Some of these areas are marked with crosses on this figure. As a conclusion, approximation of a stable model using approximations (15) or (18) can produce an unstable model if some eigenvalues are located outside the intersections of the various stability domains.

To avoid such a situation, a careful analysis of the fractional model eigenvalues is required before approximation.

5. Conclusions

Geometric (or recursive) pole and zero distribution is a solution often found in the literature for the approximation of fractional order models. Some writers even go so far as to say that “recursivity (in the sense of geometric distribution) is fundamental in the non-integer differentiation synthesis” ([29], p. 165). Such an assertion is questionable. It is indeed shown in this paper that this distribution and the associated algorithms for its computation:

-

result in the discretization of the impulse response of a fractional model,

-

are sub-optimal,

-

are one among an infinity of other permitted distributions.

The paper thus proposes several other distributions for the approximation of fractional operators and fractional models. It also shows that the geometric distribution found in the literature leads to stability issues after approximation that can be avoided by an analysis of the pole location (in the complex plane) of the fractional model before approximation.

For all these reasons, it cannot be said that this “recursivity” is the physical reason for the observed long memory behaviors as is sometime claimed for the modeling of capacitors [29] or other systems, and the author is currently exploring stochastic reasons.