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

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.


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.

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: 1 s ν and D ν (s) = s ν with 0 < ν < 1. (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: End for 3. Compute (5) 4. Define the fractional integrator (1) approximation in the frequency band [ω b , ω h ], by the transfer function As shown in [26], with this algorithm, the corner frequencies ω k and ω k (respectively the poles and zeros of the transfer function I ν N (s)) are geometrically distributed to obtain the required frequency behavior.
Figure 1 highlights this distribution and also compares the asymptotic Bode plots of I ν (s), I ν a (s) and I ν N (s).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: This leads to the following algorithm (Algorithm 2).

End for 3. Compute
For a fractional differentiation operator defined in the Laplace domain by the frequency truncation leads to For a fractional differentiation operator defined in the Laplace domain by the frequency truncation leads to and permits the following algorithm (Algorithm 3).
Algorithm 3. Fractional differentiator approximation-first method 1. Initialize End for 3. Compute C 0 with relation (5) 4. Define the fractional differentiator (12) approximation in the frequency band [ω b , ω h ], by the transfer function By analogy with Algorithm 2, the following algorithm (Algorithm 4) permits to get an integer differentiation like asymptotic behavior in low and high frequency.
End for 3. Compute C 0 with relation (5) 4. Define the fractional differentiator (12) approximation in the frequency band [ω b , ω h ], by the transfer function The set of algorithms previously established for a real fractional differentiation order can be extended to a complex fractional differentiation order ν = a + ib [21].

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: The set of algorithms previously established for a real fractional differentiation order can be extended to a complex fractional differentiation order ν = a + ib [21].

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: where 1  L denotes the inverse Laplace-Transform.The residue method for the computation of the inverse Laplace-Transform leads to [3]: The Laplace-Transform of ( 19) is thus given by: Using the change of variable Discretization of integral (22) leads to: in which 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: where previously established for a real fractional differentiation order can be ional differentiation order ν = a + ib [21].
onal Integrator by a Recursive Distribution of Poles and Zeros: An Analytical given in relation (1) using Algorithm 1, the analytical demonstration of ented in the previous section was done in [21] by considering the limit s 0 and ωh tends towards infinity.However, a simpler demonstration se of a fractional operator is now used.This impulse response is defined erse Laplace-Transform.The residue method for the computation of the leads to [3]: of ( 19) is thus given by: riable al (22) leads to: sampling interval.Relation ( 22) also highlights the link between a −1 denotes the inverse Laplace-Transform.The residue method for the computation of the inverse Laplace-Transform leads to [3]: The Laplace-Transform of ( 19) is thus given by: Using the change of variable x = e z and thus dx = e z dz, then: Discretization of integral ( 22) leads to: in which ∆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: roach For the approximation given in relation (1) using Algorithm 1, the analytical demonstration of graphical approach presented in the previous section was done in [21] by considering the limit where ωb tends towards 0 and ωh tends towards infinity.However, a simpler demonstration d on the impulse response of a fractional operator is now used.This impulse response is defined L denotes the inverse Laplace-Transform.The residue method for the computation of the rse Laplace-Transform leads to [3]: The Laplace-Transform of ( 19) is thus given by: Using the change of variable Discretization of integral ( 22) leads to: hich z  denotes the sampling interval.Relation ( 22) also highlights the link between a ional integrator impulse response and the Prony series as the inverse Laplace transform of (23) s to: . Relation (24) also highlights that the poles αk (always greater than 0) are linked by the ratio αk =e −Δz .Regarding the zeros, to the best of the author's knowledge there is no demonstration of a with the same ratio as that imposed by Algorithm 1.However, by applying a partial fraction mposition to relation (3) it can be written as: with α k = e k∆z and A k = e (1−ν)k∆z ∆z.Relation (24) also highlights that the poles α k (always greater than 0) are linked by the ratio 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: 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.
Algorithms 2018, 11, x FOR PEER REVIEW 7 of 17   with ' and ' ' 1 1 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 Ak 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.

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

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 6 ], with ν = 0.3.Using Algorithm 1 with N = 10 poles and zeros, the following factors are obtained: r = αη = 3.9811 η = 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  are now respectively given by: -3 -3 max ( ) 2.8689 10 max ( ) 2.5622 10  are now respectively given by: -3 -3 max ( ) 2.8689 10 max ( ) 2.5622 10  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 ] ,..., [ , to define these generalizations, let us define the ratio r such that: For a geometric distribution, the following relation holds: It can be generalized, among many others, by the following distributions: In Figure 6, these two distributions are compared with the geometric one.It shows that distribution D1 makes it possible to bring the first poles closer together while distribution D2 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: 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 [ω 0 , . . ., ω N−1 ], to define these generalizations, let us define the ratio r such that: For a geometric distribution, the following relation holds: It can be generalized, among many others, by the following distributions: In Figure 6, these two distributions are compared with the geometric one.It shows that distribution D 1 makes it possible to bring the first poles closer together while distribution D 2 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: With M = 2, to ensure that corner frequencies really belong in [ω b , . . ., ω h ], the following constraints must be imposed: P(0) = 0 and P(N) = N.This leads to imposing a 0 = 0 and a 1 N + a 2 N 2 = N, thus a 2 = (1 − a 1 )/N.By varying parameter a1, Figure 7 shows the large number of distributions that can be obtained in both cases.With M = 3, to ensure that corner frequencies really belong in [ω b , . . ., ω h ], the same constraints must be imposed leading to a 0 = 0 and a 1 + a 2 N + a 3 N 2 = 1.Moreover, if it is imposed that distribution D 3 fits the geometric distribution in the middle of the interval [ω b , . . ., ω h ], the following constraint must be added: P(N/2) = 1/2 and thus a 1 2 + a 2 N 4 + a 3 N 2 8 = 1/2.Finally, parameters a 2 and a 3 , meet the following system of equations: By varying parameter a 1 , Figure 7 shows the large number of distributions that can be obtained in both cases.By varying parameter a1, 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

Fractional Model Approximation and Stability Issues
We now consider a fractional order model described by the pseudo state space description, In comparison with Algorithm 1, these values are greatly reduced with Algorithm 5.

Fractional Model Approximation and Stability Issues
We now consider a fractional order model described by the pseudo state space description, where x(t) ∈

Fractional Model Approximation and Stability Issues
We now consider a fractional order model described by the pseudo state space description, 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: where matrix Ap is diagonal and can be written . Using the Laplace transform, the characteristic equations of model (45) are thus 0, [1.. ] with approximation ( 15) or ( 18), these characteristic equations become: n is the pseudo state vector, ν is the fractional order of the system and A ∈

Fractional Model Approximation
We now consider a fractional o where  

Fractional Model Approximation and Stability Issues
We now consider a fractional order model described by the pseudo state space description, 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: where matrix Ap is diagonal and can be written . Using the Laplace transform, the characteristic equations of model (45) are thus 0, [1.. ] with approximation (15) or ( 18), these characteristic equations become: nxm , C ∈

Fractional Model Approximation and Stability Issues
We now consider a fractional order model described by the pseudo state space description, 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: where matrix Ap is diagonal and can be written . Using the Laplace transform, the characteristic equations of model (45) are thus 0, [1.. ] with approximation (15) or ( 18), these characteristic equations become: pxn and D ∈

Fractional Model Approximation and Stability Issues
We now consider a fractional order model described by the pseudo state space description, 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: where matrix Ap is diagonal and can be written . Using the Laplace transform, the characteristic equations of model (45) are thus 0, [1.. ] with approximation (15) or (18), these characteristic equations become: 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:

Fractional Model Approximation and Stability Issues
We now consider a fractional order model described by the pseudo state space description, 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: where matrix Ap is diagonal and can be written . Using the Laplace transform, the characteristic equations of model ( 45) are thus 0, [1.. ] with approximation (15) or (18), these characteristic equations become: To explain how the stability domain is transformed by approximation (15) in Algorithm 3 or (18) in Algorithm 4, it 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: where matrix A p is diagonal and can be written Using the Laplace transform, the characteristic equations of model (45) are thus with approximation (15) or (18), these characteristic equations become: 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.
Algorithms 2018, 11, x FOR PEER REVIEW 13 of 17 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 ν .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 i V   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, 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 → V λ 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.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.
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.

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: To avoid such a situation, a careful analysis of the fractional model eigenvalues is required before approximation.

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.

4 .
Define the fractional integrator (1) approximation in the frequency band [ωb, ωh], by the transfer function

Figure 2 .
Figure 2. Comparison of coefficients A′k in relation (25) which tends towards coefficients Ak in (24) and zoom inside the figure.

Figure 2 .
Figure 2. Comparison of coefficients A k in relation (25) which tends towards coefficients A k in (24) and zoom inside the figure.

Figure 4 .
Figure 4. Gain plot of relative error and absolute error.Then, an optimization program aiming at reducing the relative error ) (s E rel by looking for the best location of the N poles and zeros in the transfer function

Figure 4 .
Figure 4. Gain plot of relative error and absolute error.Then, an optimization program aiming at reducing the relative error ) (s E rel by looking for the best location of the N poles and zeros in the transfer function

Figure 4 .
Figure 4. Gain plot of relative error and absolute error.

Figure 5 .
Figure 5. Gain plots of relative error and absolute error after optimization.

Figure 5 .
Figure 5. Gain plots of relative error and absolute error after optimization.

Figure 6 .
Figure 6.Two examples of pole distributions compared to the geometric one.

Figure 8 .
Figure 8.Comparison of relative (left) and absolute (right) errors between Algorithms 1 and 5 with distribution D3.

Figure 8 .
Figure 8.Comparison of relative (left) and absolute (right) errors between Algorithms 1 and 5 with distribution D 3 .

Figure 8 .
Figure 8.Comparison of relative (left) and absolute (right) errors between Algorithms 1 and 5 with distribution D3.

Figure 8 .
Figure 8.Comparison of relative (left) and absolute (right) errors between Algorithms 1 and 5 with distribution D3.

Figure 8 .
Figure 8.Comparison of relative (left) and absolute (right) errors between Algorithms 1 and 5 with distribution D3.

Figure 8 .
Figure 8.Comparison of relative (left) and absolute (right) errors between Algorithms 1 and 5 with distribution D3.

Figure 8 .
Figure 8.Comparison of relative (left) and absolute (right) errors between Algorithms 1 and 5 with distribution D3.

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

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 jω, ω ∈ [−∞, ∞] for approximation(15)   outside the domain on the right of the curve C 0 sN jω, ω ∈ [−∞, ∞],for approximation(18).
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 Ds defined by: 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 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