Impulse response of commensurate fractional-order systems: multiple complex poles

The impulse response of a fractional-order system with the transfer function sδ/[(sα-a)2+b2]n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s^{\delta }/{[(s^{\alpha }-a)^2+b^2]^n}$$\end{document}, where n∈N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n \in \mathbb {N}$$\end{document}, a∈R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a \in {\mathbb {R}}$$\end{document}, b∈R+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b \in {\mathbb {R}}^+$$\end{document}, α∈R+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \in {\mathbb {R}}^+$$\end{document}, δ∈R\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta \in {\mathbb {R}}$$\end{document}, is derived via real and imaginary parts of two-parameter Mittag-Leffler functions and their derivatives. With the aid of a robust algorithm for evaluating these derivatives, the analytic formulas can be used for an effective transient analysis of fractional-order systems with multiple complex poles. By some numerical experiments it is shown that this approach works well also when the popular SPICE-family simulating programs fail to converge to a correct solution.


Introduction
Evaluation of impulse or transient response of a linear fractional-order system from its transfer function is often part of designing and verifying the operation of analog circuits, especially frequency filters [1], oscillators [2], and generalized PID controllers (i.e., controllers including non-integer-order elements) [3].
In principle, this calculation can be performed in three ways: 1) in the environment of a standard analog circuit simulation program, as for instance SPICE; 2) using a universal program for scientific and technical calculations, e.g. MATLAB; 3) using customized software, e.g. SNAP, for the analysis of fractional-order circuits [4].
SPICE offers the possibility of a direct time response evaluation from a transfer function specified via controlled Laplace sources. However, serious numerical errors may arise because SPICE is not primarily designed for analyzing systems with dynamics of non-integer order. Another solution may be to use known efficient algorithms (or develop new ones) for transient analysis of fractional-order systems and implement them, for example, in MATLAB. An analog designer may prefer the analysis in a different software environment that would allow working in a similar style to SPICE, but with added capabilities such as the so-called (semi)symbolic analysis of the transfer function from the circuit diagram [4] or software-generated analytical formula of the impulse response using Mittag-Leffler (ML) functions [3].
The computation of the impulse response of a linear commensurate fractional-order system whose transfer function K (s) contains the s operator with non-integer powers of the form α k = p k /q k , p k ∈ N 0 , q k ∈ N, can then proceed in the following steps: -converting the transfer function to the domain of the operator w = s α with α = 1/q and q the least common multiple of the q k coefficients [5]; -decomposing the transfer function K (w) into partial fractions; -converting each partial fraction from Laplace to time domain by means of ML functions and their derivatives [3].
For the case of real simple and multiple w-domain poles, commonly used correspondences between the Laplace and time domains are available, which can be summarized by the correspondence [6] α,β (z) denotes the m-th order derivative of the ML function (see Section 2).
The case of complex poles corresponding to the denominators of w-domain transfer functions with second-order polynomials has received marginal attention in the literature so far. One possible reason is that the corresponding impulse responses can be obtained indirectly also using (1.1) if decomposition into partial fractions is performed for the corresponding conjugate pairs of complex poles, i.e. with the denominators in the form of binomials with complex coefficients instead of trinomials with real coefficients. One of the few attempts at a direct solution is given in [7] for simple poles, where a correspondence is found to exist between the Laplace transform of the form b/(s 2α + as α + b) and the corresponding impulse response formula expressed via the ML function. In [8], impulse responses for the transfer function 1/(s 2α +as α +b) γ are analyzed, but only for 0 < γ < 2. Moreover, the impulse response formula is found only for a special, not a general, constellation of parameters. It is also known that the impulse response corresponding to the "three-term transfer function" or "three-term differential equation" can be expressed using an infinite series of ML functions [6]; however, using series of ML functions is not very efficient for numerical computations.
There are two main contributions in this work: -finding the missing Laplace-time correspondence for multiple complex poles; -suggesting a reliable algorithm for evaluating two-parameter ML functions and their derivatives.
This builds a basis for an algorithmic generation of the analytic formulas of the impulse or step responses, involving the basal functions (ML functions and their derivatives as functions of time) associated with arbitrary transfer functions K (w). By evaluating these formulas, numerical results are then obtained.
Since the relation (1.1) is used to derive the above missing correspondence, its validity is first analyzed in Section 3.1 with respect to the range of values of the relevant parameters a, α and β. Then, in Sections 3.2 and 3.3, the relationships between Laplace transforms and time responses are derived for a conjugate pair of simple and also multiple complex w-domain poles. Section 4 shows examples of the numerical analysis of time responses, using the derived relations and algorithms for evaluating two-parameter ML functions and their derivatives. A comparison with the outputs of three types of SPICE simulators illustrates the robustness of the proposed procedure even in cases of a fatal failure of the SPICE algorithms.

The Mittag-Leffler function and its derivatives
For α ∈ C with Re(α) > 0, β ∈ C, and any complex argument z, the ML function is defined by means of its series expansion and it is an entire function of order 1/ Re(α) and type 1, [9]. Its derivatives can be easily determined after derivation of each term in the series as where m ∈ N 0 . For the purpose of analyzing fractional-order systems, one can restrict the definition domain of parameters α and β to the set of real numbers. Then (2.1) and (2.2) hold for α ∈ R + , β ∈ R and m ∈ N 0 .
Derivatives of the two-parameter ML function (2.2) can be more conveniently represented in terms of the three-parameter ML function (also known as the Prabhakar function [10,11]) and since for m ∈ N 0 it is Γ (m + 1) = m!. Specific numerical methods [12][13][14][15][16][17][18][19] have been devised to evaluate ML functions and their derivatives. The availability of corresponding MATLAB implementations [20][21][22][23] allows to apply the outlined procedure to the transfer function and determine the impulse and step response.

Impulse responses in terms of ML functions
To provide a proper representation of the impulse response of a linear commensurate fractional-order system we examine different cases according to the nature of the poles of the transfer function.

Multiple real poles
After substituting (2.2) into the right-hand side of (1.1) and applying the Laplace transform, we get The above infinite series converges absolutely, and for sufficiently large |s|, namely when Hence, after rearranging some terms, one easily obtains which is the relation (1.1) holding for arguments a ∈ C as well. It follows from the above derivations that, when condition (3.1) is satisfied, and if we restrict (just for simplicity) to real parameters α ∈ R + and δ ∈ R, the relations (1.1) and (3.2) hold in general for any k ∈ N and δ < kα. Thus, (1.1) can be also used for calculations of time responses for complex w-domain poles. Parameters α ∈ R + and δ ∈ R can be however chosen independently, without the additional coupling condition δ < kα. This means, for instance, that (3.2) also holds for δ ≥ kα when the order of the numerator of the transfer function is greater than or equal to the order of the denominator, which is an indication that the corresponding time response will contain derivatives of the Dirac impulse of non-integer order.
The q-th order derivative of the Dirac impulse can be indeed written in the form [24,25] and these derivatives can be directly included in the time function on the right-hand side of (1.1) when δ > kα.
This procedure can be illustrated by a simple example for k = 1 and δ ≥ α via the following arrangement: The first term in the time response represents the derivative of the Dirac impulse of order δ − α. If the inequality δ − α ≥ α (or, equivalently, δ ≥ 2α) holds, then this would mean that the second term could be used to split off another derivative of the Dirac impulse of order δ − 2α.
The same procedure can be repeated until the time response is decomposed into a sequence of derivatives of Dirac impulses and a residual function whose corresponding Laplace transform has a pseudo-polynomial in the denominator of higher order than in the numerator. Then the ML function occurring in the residual function has the second parameter greater than 0.

Simple complex poles
To analyze transfer functions with simple complex poles, let us consider a conjugate pair of complex w-domain poles with multiplicity 1, i.e.
and the corresponding transfer function with its decomposition into partial fractions is After substituting w = s α , and imposing the condition After multiplying the Laplace transform (3.4) by the term s δ and applying the relation (1.1), we get For any z ∈ C, by writing z = ρ cos θ + iρ sin θ , with ρ = |z| and θ = arg(z), it is elementary to observe that and, hence, a minor rearrangement yields the correspondence which is the basic formula for the correspondence between the Laplace and time domains for the case of simple complex poles. Since numerical algorithms for the computation of ML functions [21,23] usually make use of complex arithmetic, the relation (3.6) can be directly used for efficient time response computations.

Multiple complex poles
For the analysis of transfer functions with multiple complex poles of multiplicity n ∈ N, we consider again the conjugate pair of complex poles (3.3) and the transfer function where (3.8) Each of the terms in the right-hand side of (3.7) can be converted to the time domain thanks to the relation (1.1), after replacing the real pole with a complex one and considering the convergence condition (3.5), i.e.
By using (2.2), and following the same reasoning in the previous Section, it is possible to see that where, again, we exploited the representation z = ρ cos θ + iρ sin θ for z ∈ C. Therefore, we are able to evaluate α,α−δ (ct α ) k odd thanks to which formula (3.9) can be rearranged into the series whose last n-th term contains the n-th derivative of the two-parameter ML function of the real or the imaginary part of the argument, depending on whether n is even or odd. Coefficients A k are derived from the P k coefficients (3.8) as follows For convenience, the above series can be written in the more compact form with [x] denoting the integer part of x.

Numerical simulations
The relationships (3.6) and (3.11) between the Laplace and the time domains for complex poles will be used in this section to test the accuracy and robustness of the numerical transient analysis of fractional-order systems with different types of transfer functions, with both simple and multiple w-domain poles.
To exploit the relations (3.6) and (3.11) numerical algorithms are needed to evaluate two-parameter ML functions, and their derivatives, with complex arguments and for different sets of real parameters α and β.
Despite the special role played by ML functions in analysis and computation of fractional-order systems (essentially, the same role of the exponential function for integer-order systems), programming languages and computational packages usually do not provide built-in functions for the computation of ML functions and their derivatives; therefore it is necessary to rely on routines devised in the context of other research works.
One of the first available codes is the Matlab mlf.m function by I. Podlubny and M. Kacenak [23]; this is a code based on a mix of asymptotic expansions and numerical integration of an integral representation of the ML function, in accordance with the analysis previously provided in [14]. It is of historical importance since it has been the first available Matlab tool for computing ML functions and it made possible to begin investigations based on this function. However the code [23] is not designed to evaluate derivatives of the ML function and therefore it is not suitable to handle systems with multiple poles.
In [3] it was devised an algorithm to evaluate derivatives of the ML function by truncating the series representation (2.2); this method, however, can suffer from unexpected numerical cancellation due to the possible presence of series terms with large absolute value and alternating sign.
The problem of specifically computing derivatives of the ML functions was faced in [13] as well; essentially this work concerns with the computation of ML functions with matrix arguments, by exploiting the Schur-Parlett algorithm [26], which, in turns, demands for the computation of (possibly high-order) derivatives of the ML function in arguments related to the eigenvalues of the matrix argument, and, hence, in the whole complex plane. The algorithm proposed in [13] is based on the numerical inversion of the Laplace transform of the derivative of the ML function, thus generalizing previous algorithms devised in [12] and [19]. To improve stability of the algorithms, when high-order derivatives are requested, the algorithm in [13] exploits a formula for the representation of higher-order derivatives in terms of lower-order derivatives. The Matlab function ml_deriv(z, α, β, k), extracted from the code devised in [13], and evaluating the k-th order derivative of the ML function with two parameters α and β at each entry of the vector z, is therefore used in the numerical experiments and the corresponding results are simply denoted as ML.
The results are hence compared with those obtained from the transient analysis in popular SPICE-family simulating programs, namely PSPICE, LTSpice, and Micro-Cap (and labelled, respectively, as PS, LT and MC), which enable a direct specification of the transfer function of a fractional-order system using the controlled Laplace source. SPICE evaluates the respective time waveforms via the convolution of the driving signal and the impulse response, the latter being computed from samples of the corresponding frequency response. The transient analysis of classical integer-order systems then leads to known errors associated with specific types of transfer functions, leading to a "no causal" impulse response.
As a benchmark for the transient analysis, we consider the transfer function with the parameters a = −628.319, b = 6252 chosen so that the transfer function (4.1) corresponds for n = 1 and α = 1 to a 2-nd order system with a characteristic frequency of 1 kHz and a quality factor of 5. The time waveforms corresponding to (4.1) for this constellation of parameters are well known and the results of the transient analysis are therefore easy to verify.
Errors in the transient analysis of more complex transfer functions can also be detected by testing the limiting values of the response for t → 0 + and t → ∞. The denominator in (4.1) has all s-roots in the closed left-hand complex half-plane if Then, according to the final value theorem, the time waveform f (t), corresponding to the Laplace transform (4.1), for t → ∞ converges to 0 when δ > −1, and to 1 when δ = −1. For δ < −1 this limit is improper since the term s δ then represents a multiple pole at zero.
In the transient analysis, the fractional-order system modeled by the transfer function (4.1) was driven by a unit step. Then, to obtain the impulse response, it is necessary to multiply the transfer function by the s operator. Since different modifications of the corresponding algorithms are implemented in the above SPICE-family programs, it was interesting to compare the results obtained.
The comparisons are shown in Figures 1 and 2 in a uniform coordinate system, namely the response f (t) vs. time t.
The results (a)-(f) in Figure 1 for n = 1 are obtained for the fixed parameters α, a, b while varying the parameter δ. They confirm the stable behavior of procedures based on the computation of the ML function. The case (a) represents a double pole at zero, which leads to unbounded growth of the response over time. From among the SPICE programs, only Micro-Cap correctly evaluates this unstable behavior. PSpice generates a fundamentally wrong response and the simulation in LTSpice does not run at all (the Laplacian is singular at DC).
The simulations (b) and (c) satisfy the test for the limiting values f (0 + ) = 0 and f (∞) = 1. For the case (d), i.e. δ = 0.8, the limiting value changes at the origin to f (0 + ) = a 2 + b 2 = 3.948 × 10 7 , which is confirmed for all the computational methods. However, the results show that all SPICE simulators introduce an error in the steady state calculation.
The simulation for the parameters in the plot (e) for δ = 1 is handled by SPICE programs only with large errors, while the ML algorithm based on the Mittag-Leffler functions gives accurate results. When the parameter δ is increased further (see the plot (f)), the SPICE programs start to fail fatally. Only the ML method exhibits problemfree behavior.
Plot (g) shows the simulation results for a = 990, which is just below the 990.22 threshold between the stable and unstable response as given in (4.2). The corresponding oscillating response is confirmed by the ML method as well as by the PS and MC simulations. The PS simulation fails in this case. Finally plot (h) shows a more complex case of unstable behavior, simultaneously induced by a conjugate pair of complex poles in the instability region for a = 1500 and a double pole at zero. The response is correctly calculated only by the ML and MC methods. Figure 2 shows examples of testing the methods of transient analysis for the more complex case of multiple poles, in this case n = 3.
Plots (a)-(g) are again ordered for increasing values of δ, which implies a successive differentiation of the responses with respect to time. All tested methods show satisfactory results for δ = −1 and δ = 0. Outside this interval errors can be successively observed for all SPICE programs without exception. When δ = 4.4 is reached (see plot (e)), the limiting value of f (0 + ) goes from 0 to (a 2 + b 2 ) 3 = 5.972 × 10 22 , which is confirmed by the ML, LT and MC methods. The PSPICE does not converge for these parameters and even for higher values of δ. Plots (f) and (g) exhibit a high sensitivity of the results from LT and MC to δ increasing above the critical value of 4.4. Then only the ML method generates usable outputs.
Figure (h) finally shows the simulation of system (4.1) for n = 3 when it is set to be close to the instability limit by the parameter a = 990. Only the ML and MC methods passed this challenging test successfully.
The outputs of these simulations for gradually increasing values of δ also illustrate the robustness of the balancing derivatives algorithm implemented in [13] and used by the ML method; this method indeed generates correct results even for high values of δ, when instead SPICE programs fail.

Concluding remarks
The dictionary of current Laplace-time correspondences for fractional-order systems is completed by adding the missing correspondences (3.6) and (3.11) for single and multiple w-domain poles. This addition now allows the automated generation of analytical formulas of impulse responses using two-parameter ML functions and their derivatives for any constellation of poles, i.e., real and complex, simple and with arbitrary multiplicity. It can also be potentially used for developing useful software tools for the symbolic analysis of analog networks containing fractional-order elements.
The MATLAB code for evaluating the derivatives of two-parameter ML functions, based on the algorithm described in [13], is used and its robustness is demonstrated on benchmark simulations.

Declarations
Competing interests The authors have no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.