Optimal approximation of analog PID controllers of complex fractional-order

Complex fractional-order (CFO) transfer functions, being more generalized versions of their real-order counterparts, lend greater flexibility to system modeling. Due to the absence of commercial complex-order fractance elements, the implementation of CFO models is challenging. To alleviate this issue, a constrained optimization approach that meets the targeted frequency responses is proposed for the rational approximation of CFO systems. The technique generates stable, minimum-phase, and real-valued coefficients based approximants, which are not always feasible for the curve-fitting approach reported in the literature. Stability and performance studies of the CFO proportional-integral-derivative (CFOPID) controllers for the Podlubny’s, the internal model control, and the El-Khazali’s forms are considered to demonstrate the feasibility of the proposed technique. Simulation results highlight that, for a practically reasonable order, all the designs achieve good agreement with the theoretical characteristics. Performance comparisons with the CFOPID controller approximants determined by the Oustaloup’s CFO differentiator based substitution method justify the proposed approach.


Introduction
Applications of fractional calculus, the branch of mathematics that is regarded as a generalization of classical calculus, are ever-increasing across all domains [37,38]. Fractional-order (FO) differentiation/integration (differintegration) generalizes the traditional derivation from an integer to any real or complex number [58]. Various definitions of FO derivatives are available in the fractional calculus literature. Caputo's definition of fractional derivative allows a similar representation of initial conditions to that of the integer-order derivative because initial conditions are determined by integer derivatives. Hence, it provides an interpretation of initial conditions to be used in real-world problems more easily than the Riemann-Liouville definition [36]. However, many researchers dispute the ability of Caputo's definition to take properly into account initial conditions if used to define or simulate fractional-order systems. In synthesis, there is a discussion on the respective benefits and drawbacks of the two mentioned (or other) definitions. In detail, the Caputo and the Riemann-Liouville derivatives of a function f (t) are respectively defined as (1.2) where n is an integer, n − 1 ≤ α < n, and (·) denotes the gamma function [17]. From the input-output control point of view, initial conditions are usually set to zero such that transfer functions are considered, then the problem of initial conditions is not significant. The dynamical systems can be represented by FO differential equations; hence, the linear time-invariant models can be represented by using fractional-order transfer functions (FOTFs). Due to the presence of additional tuning knobs, the FO models can more effectively capture the dynamics of real-world systems. Applicability of FO systems is reported in control systems, signal processing, circuit theory, biomedical engineering, etc. [26,56].
The majority of the literature has dealt with the design and implementation of FOTFs comprising real-valued exponent of s. The simplest transfer function in this regard is H (s) = s α , which can be obtained by taking the Laplace transform of (1.1) and (1.2) with zero initial conditions. Assuming α ∈ + , we may write F(s) = s p s α− p = s p s γ , where p is an integer such that p = α . Therefore, γ = α − p implies γ ∈ (0, 1). Fitting, numerical, optimization, and continuedfraction expansion techniques help to approximate the frequency responses of s γ using a finite-order rational transfer function for a specified bandwidth [1,19,21,43,44,48,54]. From the perspective of signal processing, F R (s) = s γ represents a realorder FO differentiator (RFOD) for γ ∈ (0, 1), whereas γ ∈ (−1, 0) corresponds to its inverse function [33]. The transfer function of the classical proportional-integralderivative (PID) controller has been generalized as C R (s) = K p + K i s −λ + K d s δ , where 0 ≤ λ, δ ≤ 1; K p , K i , and K d are the proportional, integral, and derivative gains [55]. The traditional filters, including the special types, are also extended to the FO domain [6]. The FOTF can be practically implemented by replacing one or more traditional capacitor/inductor with the constant phase element (CPE), such as the FO capacitor/inductor [15]. Several works experimentally demonstrated the performances of FO circuits and systems (controllers, filters, etc.) employing the CPEs [16,28]. The CPEs are not yet available as commercial devices in the market; though, proposals to implement the FO systems for industrial applications are drawing attraction [50,66]. The impedance characteristics of the CPE can be emulated using passive or active emulator circuits [30,31]. However, the hardware complexity of the emulators is substantial compared to a single element fractor.
An alternative approach that avoids using CPEs in realizing the FO systems involves determining the integer-order approximation of the FOTF. Such models may be obtained through direct substitution of the rational approximant of s γ operator in the FOTF [67]. For the non-commensurate type or more complicated FOTFs, the substitution method incurs large overhead. Hence, numerical or optimal approaches are the preferred approximation tools [11,42].

Complex fractional-order controller
The generality introduced by fractional calculus allows complex FOTF-based system modeling, i.e., the exponent of s may be a complex number. For instance, the RFOD is a subset of the complex FO differentiator (CFOD) F C (s) = s α+ jβ , where α ∈ [0, 1] and β ∈ [68]. The complex FOTFs provide more degrees-of-freedom in modeling as compared to their real-order counterparts. The complex-order derivative of a sine function at steady-state is expressed as where j = √ −1 and ω is the angular frequency variable expressed in radians per second (rad/s) [40].
The third-generation CRONE controller is a seminal work that deals with the design and application of complex-order models [35]. Other application areas of complexorder derivatives can be found in the locomotion control of robots [51,61], chaotic oscillators [52,53], system identification [5,29], filter modeling [2], viscoelasticity [8], optimization algorithms [7,39] and neural networks [32]. The effectiveness of complex fractional-order PID (CFOPID) controller was demonstrated on time-delay process [3]. The improved time-domain responses yielded by the genetic algorithm based complex FO controllers as compared to the classical and real-order FO controllers were justified [64]. Tuning rules for the CFOPID and CFOPI controllers were proposed in [23,59]. Modeling of CFODs and CFOPIDs in the discrete-time domain was also carried out [9,40,57]. Stability analysis of complex-order fractional difference equation was reported in [10], while [47] investigated the geometrical and physical interpretation of the complex-order fractional integral.

Motivations for this work
The effectiveness of complex FOTF modeling has not been experimentally validated since implementing such transfer functions requires the complex-order fractance elements; realizing such elements has not yet received much attention in the literature [22,60]. Therefore, the rational approximation of the complex FOTFs remains a feasible choice. The motivation for carrying out this research is based on the limitations of the existing methods, as discussed below.

Limitations of fitting routines
Curve-fitting techniques were employed for the approximation of analog CFOD [14], CFOPID [13], and complex-order filter [12]. The fitting technique is based on the Sanathanan-Koerner algorithm with Levy's cost function [14]. The fifth-order approximation of the CFOPID controller reported in (28) of [14] is given by Vector Fitting (VF) [24,25] is another popular technique to fit the frequencydomain responses with rational functions. The VF method is applied for a bandwidth of [0.001, 1000] rad/s with 1000 sample points to generate the fourth and fifth-order approximations of the theoretical CFOPID controller given in (1.4). The corresponding approximants are given by   However, the presence of two right-half plane zeros (occurring in complex-conjugate form) for both the approximants highlight non-minimum phase behavior. Such systems may lead to stability issues in a feedback control mechanism. The magnitude and phase plots for the VF-based approximants are presented in Fig. 1. It can be seen that both the approximants fail to attain good conformity with the theoretical response.

Limitations of a built-in MATLAB function
Using the MATLAB function invfreqs with minimization of the norm of the gradient (with parameter iter set as 30 and considering 1000 sample points), the fourth and fifth-order approximations of the theoretical CFOPID controller function given by (1.4)   Although stable and minimum-phase approximations are obtained, this method gives an ill-conditioned function (coefficients with extremely large or small values), such that its implementation will lead to impractical values of passive components. The magnitude and phase responses of the invfreqs-based fourth and fifth-order approximants are presented in Fig. 2. These plots highlight a significant deviation from the theoretical responses for a wide range of the considered bandwidth.

Applicability of the Oustaloup's recursive algorithm
An alternative technique to approximate the theoretical CFOD characteristics was reported by Oustaloup et al. [48]. According to the Oustaloup's method, the rational approximant of a CFOD comprises of complex poles and zeros. The CFOD approximant for s v , where v = a + jb, based on this technique is given by where M is a positive integer, μ = ω h ω b , and the bandwidth of approximation is in the interval [ω b , ω h ] rad/s. The values of ω k and ω k are respectively given by . The Oustaloup's CFOD approximant may be re-written as a ratio of complex functions as given by where P (s), P(s), Q (s), and Q(s) are real coefficients based polynomials. Alternatively, (1.12) may be written as A French patent [49] employing operational amplifiers for the circuit implementation of (1.13) was reported three decades ago. However, it may be more straightforward to practically implement a real-valued function as compared to (1.13), which is a complex function. Moreover, as shown in Sect. 3, the Oustaloup's approximation even with M = 1 may lead to a CFOPID controller of high order, with coefficients having large values.
Therefore, it is worth investigating the following question: Within the bandwidth approximation limits where six (or less) decades may be sufficient (in general) for process control applications, is it possible to achieve reasonable accuracy using a realvalued rational transfer function with smaller orders (even and odd) to approximate the characteristics of the CFOPID controller?

Contributions of this work
In this paper, the approximation of the CFOPID controller is presented wherein the approximant is strictly a real function. The real-valued coefficients of the rational approximant are determined optimally, while enforcing the poles and zeros to lie on the left-half s-plane. Thus, the proposed work helps in alleviating the issues associated with the existing approaches. The primary contributions of this work are the following: 1. A constrained optimization approach based on a state-of-the-art evolutionary method is presented that satisfies the conditions of stability, minimum-phase, and achieves real-valued coefficients real rational approximants of both even and odd orders for the CFOPID controller. 2. The feasibility and effectiveness of the proposed approach is validated on three different variants of the CFOPID controller, namely the Podlubny's form [55], the internal model control (IMC) form [65], and the El-Khazali's form [20]. Effects of the approximation order on the approximation performance are evaluated using various frequency response error metrics. 3. The accuracy of the proposed technique is compared with the Oustaloup's method. Results demonstrate that the proposed approach provides a competitive performance when compared against the much higher order of the competing designs. While the Oustaloup's CFOD approximants are substituted in the theoretical CFOPID transfer function to obtain the corresponding approximation, the proposed method directly determines the CFOPID controller coefficients.
In the rest of the paper, Sect. 2 presents the proposed optimization problem formulation and the solution methodology. In Sect. 3, design performance studies are carried out, and the simulation results are discussed. The paper concludes in Sect. 4.

Problem formulation
The generalized transfer function of a real FOTF is defined as Replacing the real-valued exponents of s in (2.1) with complex numbers, such as λ i = r i + j y i (i = 0, 1, . . ., m) and μ k = w k + j z k (k = 0, 1, . . ., n), leads to the complex FOTF, as defined below, which provides further generalization and more degrees-of-freedom than (2.1): Any classical continuous-time system can be represented as a ratio of polynomials in s, as defined by where a k (k = 0, 1, . . ., N ) and b k (k = 0, 1, . . ., N − 1) are the real-valued coefficients of the numerator and denominator polynomials, respectively, of H N P (s); and N is the order of the system. H N P (s) is bounded-input bounded-output (BIBO) stable if all the poles reside in the left-half s-plane; if all the zeros lie on the left-half s-plane, then the system attains minimum-phase behavior. In the case of complex conjugate poles or zeros, their real part must be negative. A necessary condition for BIBO stability and minimum-phase response is that all the coefficients of Q N (s) and P N (s), respectively, must possess the same sign with no missing terms.
The frequency-domain transfer functions of (2.2) and (2.3) can be obtained by substituting s with jω, as respectively defined below (2.5) Table 1 shows the transfer function and frequency domain expressions for the different CFOPID controllers considered for approximation in this work. The design problem can be solved in the optimal sense, where the purpose of the search routine is to minimize the frequency response error between the theoretical (targeted) model F C (s) and the rational approximant H N P (s), subject to satisfying the stability and minimum-phase criteria. The objective function f of the proposed constrained optimization (minimization) problem can be formulated as where L is the total number of log-spaced data points in the bandwidth of approximation [ω min , ω max ] rad/s; X represents the decision variables vector, i.e., ; σ z and σ p denote the real part of the roots of P N (s) and Q N (s), respectively. The total number of decision variables for this optimization problem is 2N + 1, which represents the dimension (D).

Solution methodology
A metaheuristic algorithm facilitates the integration of various constraint handling strategies to generate a feasible solution [34]. Although (2.6) may be minimized using Table 1 CFOPID models considered for approximation [20] El-Khazali any constrained optimization tool, this work employs the constrained composite differential evolution (C 2 oDE) algorithm [69] due to its effectiveness in solving similar problems for real FO systems [41]. Moreover, C 2 oDE is a real-parameter optimizer, implying that the decision variables will always be real numbers. Hence, the drawback of [12] concerning the generation of complex coefficients for the rational approximants can be eliminated. A brief discussion of the metaheuristic search process of C 2 oDE to solve the devised problem is presented below. C 2 oDE is a constrained evolutionary optimization technique that strives to maintain the trade-off between the diversity of solution versus convergence and satisfy the constraints without compromising the quality [69]. Under the evolutionary framework of the composite differential evolution (CoDE) algorithm [70], two popular and efficient constraint handling schemes, along with a restart mechanism, are integrated in C 2 oDE. The major phases of the algorithm are outlined as follows.

Initialization Phase
Similar to all population-based metaheuristic algorithms, C 2 oDE comprises NP search agents (also known as target vectors), x i (i ∈ 1, . . . , NP). For the proposed optimization problem, the initial population of C 2 oDE can be formulated as (2.7) All the agents are initialized with randomly chosen real numbers from a uniform distribution within the search limits in the D-dimensional hyperspace. The lower bound of the design variables is selected as 10 −5 (i.e., a positive number) to satisfy the same-sign condition for the coefficients of P N (s) and Q N (s). A small upper bound value may not provide sufficient diversity in the population. In order to limit the large spreading of coefficients, initially, the upper bound for the design variables is selected as 1000. In case the obtained coefficient(s) value gets stuck at the upper bound, then the limit is increased by a factor of 10 and the search process is repeated.

Search Algorithm
Three mutation techniques are incorporated in the strategy pool of C 2 oDE to create a mutant vector (v t i ), where i = 1, 2, …, NP, and t denotes the generation number. The reason for incorporating several mutation schemes instead of a single one is to provide more diversity and also enhance the algorithm's convergence speed. These selected strategies exhibit complementary characteristics since DE/current-to-rand/1 promotes diversity, whereas modified DE/rand-to-best/1/bin and DE/current-to-best/1/bin improve the convergence characteristics. The expressions for the three mutation schemes are given by , and x t r 3 (r 1 = r 2 = r 3 = i) are three randomly chosen target vectors; F denotes the scale factor; and rand ∈ U (0, 1), with U uniform distribution; where the agent with the smallest degree of constraint violation is indicated by x t Gbest and x t r 4 (r 1 = r 2 = r 3 = r 4 = i) denotes a randomly chosen target vector; and where x t f best is the individual with the best (minimum) values of the objective function. The expression for the degree of constraint violation for x i is given by where g j (x) < 0 represents the inequality constraints to satisfy the stability criteria. After undergoing mutation, the algorithm executes the binomial crossover on x t i and v t i (obtained in (2.9) and (2.10)) to create the trial or offspring (u t i ) where i ∈ (1, 2, . . . , NP); j ∈ (1, 2, . . . , D); rand j ∈ U (0, 1); jrand ∈ randint(1, D); and CR is the probability of crossover. Note that DE/currentto-rand/1 operation incorporates arithmetic crossover within the formula itself. So, for this mutation technique, a binomial crossover is not used. C 2 oDE incorporates two pools, F pool = {0.6, 0.8, 1.0} and CR pool = {0.1, 0.2, 1.0}, where the values of F and CR are available. Every time the algorithm needs the F and CR values for performing the mutation and crossover operations; it picks up a value from the corresponding pools (i.e., F pool and CR pool ) randomly. This idea in C 2 oDE is different from that of the CoDE algorithm, where three specific pairs of {F, CR} combinations were employed. Thus, u t i1 , u t i2 , and u t i3 , are the three trial vectors generated subsequent to the three mutation rules and crossover.

Handling of Constraints
At first, Deb's constraint handling strategy [18] is employed to determine the best solution (u t i,best ) among the trial vectors {u t i1 , u t i2 , u t i3 }. Then, the ε-constrained technique [63] is invoked to select between x t i and u t i,best . Given any two agents x i and x j , the ε-constrained method performs a pair-wise comparison and favors selecting x i over x j according to the rules framed as Here, the initial threshold ε 0 = max G(x 1 1 ), G(x 1 2 ), . . . , G(x 1 NP ) ; t max denotes the maximum number of generations; p controls the degree of exploitation of the objective function information andλ helps avoiding an infinite value for cp if the condition ε 0 = 0 occurs (i.e., no constraint violation occurs in the initial population).

Restart Mechanism
For highly complicated constraints, a restart mechanism is also integrated in C 2 oDE. This activation is triggered if no feasible solution exists or the standard deviation values for the degree of constraint violation/objective function are smaller than a specific threshold (μ). The random initialization phase is repeated if any of the aforementioned conditions exist. The steps of the C 2 oDE algorithm are presented in Table 2.
The flowchart of the proposed design technique is shown in Fig. 3. The inputs to the optimization routine are the parameters of the complex FOTF (α, β, γ , λ, θ, μ, φ, K p , K i , K d , K C , T i ) as applicable, the desired approximation order (N), L, and bandwidth of approximation ([ω min , ω max ] rad/s). Due to the stochastic nature of the search process of C 2 oDE, multiple runs of the algorithm are required to determine the best solution quality. The proposed work performs maxrun number of independent trial runs of C 2 oDE for each of the considered design examples. Each run of the algorithm terminates after completing the FE max number of objective function evaluations. For each trial run of the algorithm, the best feasible solution vector in that run and its fitness value are stored. Finally, at the end of all the runs, the best among all the stored best solutions obtained in each run, i.e., the one achieving the least value of the objective function, is selected as the near-global optimal and feasible decision variables vector X * .

Simulation results and discussion
In this work, the internal control parameters of C 2 oDE are chosen as per the recommendations provided in [69], whereas the basic parameter NP is selected as 10D as specified by a rule of thumb [62]. For demonstration purposes, the values of L, ω min , and ω max are chosen as 1000, 10 −3 rad/s, and 10 3 rad/s, respectively. The optimization Table 2 Steps of C 2 oDE algorithm Step 1. Set NP = 10D,μ = 10 −8 ,λ = 6, p = 0.5, t = 1, FE max = 10000D.
Step 3. Randomly initialize the search agents in the chromosome population, . , x t NP }, within the lower and upper bounds.
Step 6. For the ith (i = 1, 2, . . ., NP) agent, Gbest , x t r 1 , x t r 2 , x t r 3 , x t r 4 , F ∈ F pool , and CR ∈ CR pool . Acquire v t i2 as per (2.9). Perform binomial crossover on v t i2 and x t i to generate u t i2 . Select x t f best , x t r 1 , x t r 2 , F ∈ F pool , and CR ∈ CR pool . Acquire v t i3 as per (2.10). Perform binomial crossover on v t i3 and , u t i3 } using Deb's constraint handling rules. Select between u t i,best and x t i using the ε-constrained technique. FE = FE + 3.
Step 7. Initiate the restart mechanism (if applicable).
Step 9. Stop if FE ≥FE max ; else, repeat from Step 5. routine is implemented in MATLAB (software version: 2014a). To measure the accuracy of the approximants, the absolute magnitude error (AME) and absolute phase error (APE) metrics, as defined below, are used:

Performance analysis of the proposed controllers
Three different CFOPID controllers, such as those of the Podlubny's form (cases I and II) [55], IMC-based (case III) [65], and El-Khazali's form (case IV) [20], are shown in Table 3. The optimal coefficients of the proposed controllers for two different values of N for each case, are also presented in   6 }, respectively, based on the trial-and-error method mentioned previously. Note that increasing N should improve the approximation by reducing the errors in the Bode gain and phase plots. A higher order N is needed to represent the long memory of the system very accurately, but the computational burden increases. On the contrary, a lower N produces more ripples in gain and phase plots, but is advantageous, because the changes in coefficients due to tolerances of components (in an analog implementation) or due to the limitations of microprocessor words and the quantization effects (for a digital implementation) can be contained [45]. Hence a low sensitivity to parameter variations is obtained. Moreover, an accurate implementation over a wide frequency range is difficult for a conventional hardware realization due to computational complexity [27] and hardware costs can be limited if the accuracy is required in a restricted frequency range. Finally, from the practical point of view, approximations are not required over a very large frequency range. Namely, in systems and control engineering applications, approximations should work within a bounded frequency range, where a few decades are usually satisfactory. The stability and minimum-phase of the designs can be confirmed from Table 4, which shows that all the poles and zeros reside on the left-half s-plane. This is an important property since the occurrence of non-minimum-phase zeros in a controller may lead to stability issues in a closed-loop control system, including approximants of such controller.
The frequency responses of the designed controllers are compared with the theoretical models for the four cases, as shown in Figs. 4(a)-(d). It may be observed that: (i) for the case I, the sixth-order design achieves lower magnitude error than the fifthorder model in the frequency range [0.001, 0.028] rad/s. The fifth-order approximant outperforms the sixth-order one in terms of phase response accuracy in the interval [0.001, 0.145] rad/s. Frequency responses of both the approximants are similar at high frequency regions. Increasing the order of the proposed approximant further slightly reduces the error in phase approximation but does not significantly affect the error in magnitude approximation; (ii) for case II with both the approximation orders, a similar accuracy in magnitude is obtained in the frequency range [0.181, 6.884] rad/s, whereas the same accuracy for phase characteristic is obtained between 0.774 and 3.4 rad/s. Overall, an improved accuracy in magnitude and phase is respectively attained using Table 4 Locations of zeros and poles of the proposed optimal approximants of the CFOPID controllers  Table 5 Comparison of error metrics for the proposed optimal approximants of the CFOPID controllers (best performance is highlighted in boldface) the designs with N = 5 and N = 4; (iii) for case III, the frequency characteristic plots of both the approximants exhibit similarity. While the phase response exhibits proximity with the theory throughout the bandwidth, the magnitude plot starts deviating at the lower frequency end below 0.04 rad/s; (iv) for case IV, the magnitude-frequency profiles for the proposed fourth and fifth-order aproximants are similar, whereas the fifth-order design markedly outperforms the fourth-order model in the range [0.001, 0.026] rad/s regarding the phase response accuracy. Table 5 presents the performance metrics for all the designed CFOPIDs, which show that: (i) for case I, maximum APE for N = 5 (49.841 dB) and N = 6 (49.153 dB) is similar. The fifth-order design outperforms the sixth-order controller about the maximum and mean APE, whereas the sixth-order design achieves smaller mean AME; (ii) for case II, the fifth-order design is inferior to its lower-order counterpart regarding the maximum and mean APE metrics, but significantly outperforms the fourth-order approximant about both the magnitude response indices; (iii) for case III, the performances of both the designed controllers are similar; (iv) substantial improvement in maximum APE is achieved by the fifth-order approximant as compared to the fourth-order one for case IV. However, it is remarked that the high values of maximum absolute errors are typically obtained at the extremum low and high frequencies of the considered intervals, where the approximations are obviously worst. As a consequence, mean absolute errors are affected by such values. On the contrary, the approximations give suitable results inside the frequency interval [ω min , ω max ].

Comparison with the literature
The CFOPID controller can be approximated using Oustaloup's CFOD approximant by a simple substitution method. For instance, the normalized CFOD approximant for s 0.8+0.1 j based on the Oustaloup's method with M = 1, ω b = 0.001 rad/s, and ω h = 1000 rad/s, is given by  The magnitude and phase responses comparison plots of the controllers obtained by using the Oustaloup's CFOD-based substitution and the proposed optimization with larger N are illustrated in Figs. 5(a)-(d) for cases I-IV, respectively. It is observed that: 1. for case I, the proposed and Oustaloup's controller achieve superior accuracy in the magnitude and phase responses, respectively, at the lower frequency values, whereas their performance is similar in the mid-band region. At the higher frequencies, the proposed controller attains better phase accuracy, whereas the Oustaloup's approximant yields superior accuracy in magnitude; 2. for case II, the proposed controller's phase behavior is inferior to Oustaloup's.
However, the proposed approximant attains better proximity to the theoretical magnitude plot over the entire bandwidth; 3. for case III, the magnitude response accuracy of the Oustaloup's CFOPID approximant is superior to that of the designed one for frequencies in [0.001, 0.019] rad/s. The phase response of the proposed approximation suffers lower deviation in the interval [0.001, 0.64] rad/s, while the Oustaloup's model achieves better conformity with the theoretical magnitude-frequency behavior. At the higher frequencies, the behavior is similar; 4. for case IV, in the low frequency region, the proposed model yields an inferior phase response. However, the proposed approximant achieves a smaller value for the maximum deviation in phase response beyond 335 rad/s.
In Table 6, comparisons about the AME and APE metrics between the proposed and Oustaloup's CFOPID approximants are presented. Results reveal that: (i) for case I, the Oustaloup's model yields better performance for maximum and mean AME and mean APE, whereas the proposed design achieves significantly lower maximum APE; (ii) for case II, the magnitude errors are significantly smaller for the proposed design, although the Oustaloup's model outperforms the proposed one about the phase response errors; (iii) for case III, significant reduction in maximum APE (3.265 • ) is obtained by the proposed approximant as compared to the Oustaloup's approximant (36.745 • ), although the Oustaloup's design yields a better magnitude response accuracy; (iv) for case IV, the proposed model provides markedly improved accuracy over Oustaloup's for the maximum APE (21.782 • versus 44.357 • ) but is outperformed for the other three error indices.
Overall, it may be inferred that the proposed approach does not outperform the Oustaloup's method for both the magnitude and phase responses simultaneously. However, the proposed technique results in real coefficients based real rational transfer functions, which is in contrast to the real coefficients based complex transfer function yielded using the Oustaloup's method. Therefore, hardware implementation of the proposed controllers can be realized using standard circuit design techniques readily available in the literature. Furthermore, unlike the Oustaloup's method, there is no restriction on the approximation order (N : even or odd) of the proposed controllers.
Among the various techniques available for continuous-time rational approximation of FOTFs, it is not possible to establish which one is the best [46]. While some methods provide better accuracy regarding the frequency response, others may achieve improved time response performance. The approximation performances can also depend on the non-integer order. The real and imaginary part of the complexvalued impulse response of the proposed sixth-order approximant and the Oustaloup's

Conclusions
An optimal technique that guarantees stable poles, minimum-phase zeros, and realvalued coefficients of rational approximants for the complex fractional-order PID controllers is presented in this paper. The drawbacks of the reported curve-fitting techniques are eliminated through (i) incorporation of constraints, (ii) appropriate selection of the lower bound of decision variables pertaining to the coefficients of the numerator and denominator polynomials of the proposed model, and (iii) utilizing a real-parameter constrained optimization algorithm. The effectiveness of the suggested technique is verified on the Podlubny's, IMC, and El-Khazali's forms of CFOPID controllers. Performance comparisons for different orders (odd and even) of approximation are demonstrated for the design examples. Comparisons with the Oustaloup's CFOPID approximant shows that the proposed approach may achieve better magnitude or phase response fitting, though simultaneous improvements in both are Fig. 7 Impulse response (imaginary part) comparison plots between the proposed (N = 6) and Oustaloup's (N = 12) approximants for the CFOPID controller of case I not possible. The proposed approach may be considered as an effective alternative to the Oustaloup's method that results in complex rational approximants with real coefficients; especially, when the bandwidth considerations are limited to 3∼4 decades. For the considered non-linear optimization problem, two limitations of the proposed approach are (i) an inability to guarantee the generation of a global optimal solution, which is in general true when employing a metaheuristic method, and (ii) large dispersion of coefficients (especially for the denominator polynomial) when N is increased. While this work solely concentrated on the approximation of CFOPID controllers, it will be interesting to consider functions such as the unstable, non-minimum phase, and conjugated-order CFO systems [4,68] in the future. Solving such problems may necessitate incorporation of additional design constraints that may involve employing the multi-or many-objective constrained evolutionary methods.
Funding Open access publishing supported by the National Technical Library in Prague.

Conflict of interest
The authors declare that they have no conflict of interest.
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/.