Robust controller design: Recent emerging concepts for control of mechatronic systems

The recent industrial revolution puts competitive requirements on most manufacturing and mechatronic processes. Some of these are economic driven, but most of them have an intrinsic projection on the loop performance achieved in most of closed loops across the various process layers. It turns out that successful operation in a globalization context can only be ensured by robust tuning of controller parameter as an effective way to deal with continuously changing end-user specs and raw product properties. Still, ease of communication in non-specialised process engineering vocabulary must be ensured at all times and ease of implementation on already existing platforms is preferred. Speciﬁcations as settling time, overshoot and robustness have a direct meaning in terms of process output and remain most popular amongst process engineers. An intuitive tuning procedure for robustness is based on linear system tools such as frequency response and bandlimited speciﬁcations thereof. Loop shaping remains a mature and easy to use methodology, although its tools such as Hinf remain in the shadow of classical PID control for industrial applications. Recently, next to these popular loop shaping methods, new tools have emerged, i.e. fractional order controller tuning rules. The key feature of the latter group is an intrinsic robustness to variations in the gain, time delay and time constant values, hence ideally suited for loop shaping purpose. In this paper, both methods are sketched and discussed in terms of their advantages and disadvantages. A real life control application used in mechatronic applications illustrates the proposed claims. The results support the claim that fractional order controllers outperform in terms of versatility the Hinf control, without losing the generality of conclusions. The paper pleads towards the use of the emerging tools as they are now ready for broader use, while providing the reader with a good perspective of their potential.


Introduction
With the recent industrial revolution at hand, a great burden is placed on the performance of the global manufacturing process, due to continuously varying user-defined specification of the end-product. The implications are mainly economic, but the driving force consists of better process operation in the lower level loops [1] . That is, a good lower loop performance facilitates high-end performance of supervisory loops where the economic cost is most relevant. Recent surveys have outlined emerging tools in the field of process control with the major reason for their success being the intrinsic robustness they offer at a higher degree of freedom to operate and tune controller parameters [2][3][4] .
In a recent summary on the industrial relevance of various control algorithms [5] , it is overwhelming the general agreement on top-three most prevalent in practice: PID (proportionalintegral-derivative) control, MPC (model predictive control) and system identification. At the bottom of the same list are robust control mature algorithms such as Hinf and LQR (linear quadratic regulator), which are based on loop shaping tuning methods.
At the core of the robust design in control systems are the frequency based loop shaping methods [6] . This is of particular interest in smart manufacturing, whereas mechatronic applications of human-robot collaborative systems play a key role [7,8] . Hinf control is a mature methodology and has an established tradition in control theory [9] . Robustness is a key feature in the optimization procedure used to find the controller parameters. Given user defined specifications, translated in frequency response domain, the controller can be optimally tuned in frequency domain with resulting poles and zeros. Other notable applications in mechatronic system controls are based on backstepping strategies [10][11][12][13][14] , sliding mode [15] or combinations thereof with information technology tools [16] . In miscellaneous applications of such basic elements of control systems as motors, the essential features are robustness and real time implementability [17][18][19][20] .
To fulfil a high degree of robustness over a specified bandwidth in single controllers designed for process operation, the control design uses specifications in terms of open loop Table 1 Real response performances comparison between PID, FO (fractional order) PID and Hinf controllers.

Brake
Rise  frequency response. Some of these specifications are illustrated in Fig. 1 , introducing the principle of band-limited robustness (Ro) properties in the Nyquist plane. The controller parameters have to be designed in such a way as to present iso-properties within a specified frequency range, either for phase, either for magnitude intervals. In terms of fractional calculus, this is achieved by continuous fraction expansion, or, in control engineering terms, by pole-zero interlacing. The spacial distribution of the pole-zero interlacing will dictate the slope of the magnitude/phase in frequency domain [21][22][23] .
The generic property of the resulted process and controller closed loop system is robustness to gain and phase variations in the process dynamics. Various methods for pole-zero interlacing exist, all based on frequency domain specifications. A comprehensive review is given in [21] , with manifold applications in [24][25][26] , including mechatronic systems [27][28][29] . This nonrational controller needs to be deployed in real-time systems and discrete-time approximations are described in [21,[24][25][26] while an efficient, low-order, stable solution for discretizing any type of non-rational function is given in [30] .
Originally emerged from fractional calculus, fractional order control (FOC) is increasingly used in control engineering since the generalization of the PID controller using any real order was introduced by Podlubny [Podlubny1999] . Podlubly named this generalized PID controller as the fractional order controller, PI λ D μ where λ ∈ (0, 2) is the order of integration and μ ∈ (0, 2) is the order of differentiation. The superiority of this generalized PID controller in comparison with its classical integer order version has also been demonstrated [21, Podlubny1999, 32] . The design methods range from frequency domain tuning techniques [21,33,34] , to methods using advanced optimization techniques [35][36][37][38] and to methods with time domain specifications [39] . Multiple applications of fractional order controllers have been described while presenting the fractional calculus as a robust control design tool [29,[40][41][42] . Most of the fractional order controller parameter tuning algorithms use an imposed performance criteria for the robustness of the designed control system [24] .
In this paper, we provide an answer for the reader in terms of the following questions: • How does a robustly tuned PID control, known as the 'working horse for industry', perform against emerging fractional order control?. • How does Hinf control, as an example of 'classical' robust control method, perform against emerging fractional order control?.
The present paper pleads to the community for becoming aware of the powerful emerging tools for robust control system design, while summarizing the main tuning methods used in practice for basis loops in mechatronic field. Emerging control design methods such as fractional order systems and controls are introduced from a practical perspective. Next, we provide an analysis and comparison by means of an ubiquitous industrial element, i.e. a DC motor with varying load conditions. This case study has been selected being one of the core problems in control applications [43] , but without losing the generality of conclusions. All designed controllers are experimentally tested on the modular servo system designed by Inteco [44,45] .
The paper is structured as follows. The next section presents the commonly used tuning principles for PID control in industry. The third section provides the principles of robust control by means of Hinf control design method, followed in the fourth section by the principles underlying fractional order control design. A fifth section provides the result for the real life setup, followed by a section containing recommendations for the user. A conclusion summarizes the main idea of the paper. The Appendix contains essential information for those who aim embracing the novel concepts.

Common PID tuning rules in industry
These methods presented hereafter have been selected based on their industrial relevance. They have been extensively published in a manifold of works, see for instance [21,23,25,[46][47][48] . However, they have all been used in dynamic systems control with linear parameter varying conditions and thus have been tested against a great deal of robustness. These works are suitable for the context of the paper, and they are summarized here for the sake of completeness. A comprehensive method overview and comparison can be found in [49] .

Autotuning methods
The indirect tuning methods are those which prior to control parameter tuning require identification of basic step response data to first order plus dead time (FOPDT) or second order plus dead time (SOPDT) [50,51] . By contrast, direct methods skip this identification step. Fig. 2 provides this concept in graphical form.
Identification for the purpose of control is a demanding step in the process of model development. Nevertheless, latest reviews on their industrial relevance indicate that system identification plays an important role in practice [5,52] . Identification methods vary in terms of complexity, depending on the scope of the target model usage. Methods for system identification are available in both time-and frequency-domains. Most commonly used in industry are the step test response data, relay test data and sinusoid test data [1,53,54] . Suitable methods for industrial process control are event-based algorithms, where basic identification plays an important role before tuning controller parameters [55][56][57] .
Relay based methods have been some of the first used to automatic tune PID-type controllers as they were possible to perform using hydraulic and pneumatic instrumentation com- monmly available in industry [47,48] . Recent revisions thereof allow a more robustness to short data series [58,59] .
The AMIGO method uses a FOPDT process approximation of the form: with K the gain, τ the time delay and T the time constant of the approximation to step response data of the process. It then uses a parallel PID configuration with the following tuning rules: The SIMC method uses a SOPDT process approximation of the form: with K the gain, τ the time delay and T 1 > T 2 the time constants of the approximation to step response data of the process. It then uses a series PID configuration with the following tuning rules:

Frequency response tuning methods
The frequency response function (FRF) of a dynamical system is a measure of the modulus and phase of the output signal as a function of an input frequency, relative to the input signal applied to the system. To appropriately characterize the process dynamics in a given frequency interval, the gathered information must cover the modulus, phase and their corresponding slopes with respect to frequency. Classical methods for estimating FRF are based on available input and output data and application of the Fast Fourier Transform (FFT). These procedures require multiple or persistent exciting tests with input signals of various frequencies, such that the frequency response can be estimated over the required frequency range [62] .
When the FRF is required around certain frequency, it is pragmatic to reduce the number of necessary experiments, complexity and time-to-deliver by using an efficient and reliable algorithm. The minimal required information is the modulus, the phase and the corresponding frequency response slope in/around a specified frequency [43,49,63] .
A large number of applications require the frequency response slope and several detection algorithms are available. Relay-based methods are used in [48] , with the identification method being automated and thus useful for autotuning applications. The slope of the frequency response modulus is estimated at the gain crossover frequency. Computation of the frequency response slope has also found applications in the estimation of non-stationary sinusoidal parameters for sinusoids with linear amplitude/frequency modulation. An enhanced algorithm for frequency domain demodulation of spectral peaks is proposed in [64] , which delivers an approximate maximum likelihood estimate of the frequency slope.
In [65] , Bode's integrals are used to approximate frequency response slope of a system at a given frequency, without any model of the process. This information is then used to design a PID controller for slope adjustment of the Nyquist diagram and improve the overall closed-loop performance. The frequency response slope is also employed in the estimation of the gradient and the Hessian of a frequency criterion in an iterative PID controller tuning method.
Emerging industrial controllers of higher degree of freedom are fractional order controllers [3,4,24] . The phase slope of the FRF has been used in the design of fractional order PD/PI controller based on an auto-tuning method that requires knowledge of the process modulus, phase and phase slope at an imposed gain crossover frequency [32] .
A popular direct PID tuner is based on the common relay feedback test, with amplitude d . This test identifies the modulus at the critical frequency, i.e. the point of intersection with the negative real axis in the Nyquist plane. The period of oscillation T c , and the amplitude of the oscillation a are then used to tune the PID-type control parameters. The most known method is the modified Ziegler Nichols with tuning rules [47] : This commonly used in industry tuner provides a standard robustness of 0.5 (on a range from 0 to 1) in the Nyquist plane.
A method based on the same test was proposed in [66] , with a degree of freedom to specify the desired phase margin of the loop. In this way, the user may alter the robustness provided by the tuner, and vary its value according to the process dynamics. Using the same relay feedback test and specified phase margin PM value (commonly selected between 40 • and 75 • ), the tuning rules are: Notice the relation T i = 4T d ; this is a commonly used choice to simplify the tuning procedure; it assumes two identical real zeros in the PID controller. Some other choices of this ratio and their analysis can be found in [67] .

Classical robust control: Hinf controller design
A manifold of mechanical actuators and other mechatronic systems in general are controlled by robust control techniques, with prevalence of Hinf control algorithm. The reason is due to the nonlinear characteristic (static or/and dynamic) that mechanical actuators exhibit and temperature dependent dynamics. The instrumentation is often described as being part of a linear parameter varying (LPV) system, requiring high degrees of robustness to maintain performance over wide operation range. An important aspect is then played by the model mismatch or misalignment. The control objective is to design robust controllers to ensure both stability and good performances of the system despite disturbances and model uncertainty. A popular solution is Hinf-optimal control, where the model used for design incorporates the possible uncertainties, illustrated hereafter.
Consider a simplified process transfer function H (s) = K T s+1 . Representing the gain uncertainty as an upper linear fractional transformation [6] , results in the block scheme from The corresponding state space equations are: resulting the packed matrix [6] : The controller C is designed to achieve robust stability, meaning that the closed-loop system remains internally stable for all possible plant models P and satisfy the performance criterion [6] : where W p , W u are weighting functions chosen to represent the frequency characteristics of performance and control-effort constraints. To design robust controllers is possible using computer aided design tools and dedicated control systems software such as Matlab [6] . Dedicated automatic tuning functions are available within the Robust Control Toolbox in Matlab and examples are manifold. Despite this, industry relevance of robust control and Hinf control in particular is lower than of PID control, as reported in [5] .

Design principles
The fractional order PI λ D μ controller may be defined as a generalized classical PID, having an integrator and a differentiator of order λ and μ, respectively, where generally λ, μ ∈ (0, 2). The transfer function of this fractional order controller is given as: The tuning of the fractional order PID controller in Eq. (10) implies determining the five controller parameters, the proportional, integrative and derivative gains, k p , k i and k d , as well as the fractional orders of integration and differentiation, λ and μ. To determine these parameters, five nonlinear equations are used resulting from five performance criteria [21,68] . The performance criteria usually refer to the following: • a gain crossover frequency, ω gc , that ensures a certain closed loop settling time. The gain crossover frequency requirement is translated into the modulus condition as indicated next: where H ol ( s ) is the open loop transfer function.
• A phase margin PM that limits the closed loop overshoot. This is translated into the phase condition as given below: • The iso-damping property, which ensures a constant overshoot value of the closed loop system despite gain variations. This is specified in terms of the open loop phase derivative: • Good output disturbance rejection. This is specified as a constraint on the sensitivity function as follows: for all frequencies ω ≤ ω S rad/s. • High frequency noise rejection. This is specified as a constraint on the complementary sensitivity function as: for all frequencies ω ≥ ω T rad/s.
The complex frequency domain form of the fractional order PI λ D μ controller is given as: The tuning of the fractional order PI λ D μ controller is now defined as the solution of the system consisting of complex nonlinear equations given by Eqs. (11) -(15) .

Optimization
To solve Eqs. (11) - (15) and determine the controller parameter values, the Matlab optimization toolbox is considered, using the fmincon() function. This function attempts to solve problems of the form: min X F ( X ) subject to some linear constraints specified as inequalities A * X ≤ B and equalities A eq * X = B eq , nonlinear constraints specified also as inequalities C ( X ) ≤ 0 and equalities C eq (X ) = 0. The solution is bounded in the interval LB ≤ X ≤ UB , where X = [ k p k i k d λ μ] is a vector containing the five controller parameters and LB and UB are the minimum and maximum admissible values for the controller parameters.
The modulus condition in Eq. (11) is further considered as the main function to be minimized, whereas the phase margin condition in Eq. (12) and the robustness condition in Eq. (13) are defined as nonlinear equality constraints. The complementary sensitivity and sensitivity functions in Eqs. (14) and (15) , respectively, are defined as nonlinear inequality constraints. Since there are no linear constraints, the A , B , A eq and B eq parameters in the fmincon() function are not defined.
For the fmincon() function to be successful, a proper set of initial values X 0 = [0. 01 ; 3 ; 0. 5 ; 0. 3 ; 0. 5] needs to be specified. The fmincon() function uses the interior-point algorithm, which can handle large, sparse problems, as well as small, dense problems. The advantage of this algorithm is that it satisfies bounds at all iterations, and can also recover from Not a Number or Infinity results.

A basic element in mechatronic industry: the DC motor
LPV mechatronic systems are a class of (non)linear systems whose properties depend on a set of time-varying parameters that are not known in advance. A way to deal with model uncertainty is to probe the system dynamics in various conditions to extract a set of models characterizing these possible situations. This is a tedious but necessary task when accurate performance is required [69] . If performance is not the primary goal, but rather compensation of unknown disturbances, or unknown dynamics, then the control design parameters are tuned for slower control actions, with high robustness.
LPV mechatronic systems are encountered in heavy duty machines (e.g. mining, agriculture, road engineering) and in automotive industry, mainly due to changing environmental conditions [70] . Space engineering applications have increased degrees of LPV dynamics, due to their size, relative position, remote location and complexity. The mechanical design of such systems may not be optimal from control point of view, but it is from a practical point of view (e.g. remote access, testing, validation protocols, etc) and it poses additional challenges for control [71] . Moreover, many types of motors such as brushless, induction and DC motors are part of other complex systems, e.g. windmill parks, unmanned systems and parts of manufacturing and production systems [72] .
The DC motor is a versatile execution element which requires a certain degree of robustness due to varying operation conditions, load changes and other varying variables linked to it, making it a linear parameter varying system. We have explicitly chosen this example due to its simple dynamic operation. Consequently, differences among the Hinf and FOC methods will be solely due to the control strategies and their intrinsic properties, and not due to some exotic process dynamics.
The experimental unit consists in the modular servo system designed by Inteco [44] used in the particular configuration indicated in Fig. 4 . The plant is composed of a tachogenera-tor (used to measure the rotational speed), inertia load, backlash, incremental encoder, and gearbox with output disk.
The mathematical model of the modular servo system has been determined using classical experimental identification method as being: where ω is the angular velocity of the rotor, and k = 158 and T = 0. 95 s are the motor nominal gain and time constant, respectively. This model is used to design the robust controller using Hinf rules and the fractional order PID controller. For sake of comparison, a classical integer order PID controller has also been tuned using computer aided tools from Matlab, providing the parameters: K p = 0. 02, T i = 2 and T d = 0. 5 in the textbook structure. All controllers have a specified phase margin of P M = 60 • . The next step consists in the tuning of the Hinf robust controller in terms of Eq. (9) for the process described in Eq. (8) using the tuning procedure presented in Section 3 . The considered gain uncertainty is k = k N (1 + p k δ k ) , with k N = 158 , p k = 0. 25 and −1 ≤ δ k ≤ 1 .
To implement the fractional order integrator and differentiator, we use the efficient method for approximation described in Appendix [30] . The approximation parameters have a low frequency bound of 10 −2 and a high frequency bound of 10 2 rad/s, as well as an order of approximation N = 5 . To test the last two remaining performance specifications, the Bode diagrams for the sensitivity function and complementary sensitivity function need to be analyzed. The sensitivity and complementary functions are given in Fig. 6 , confirming that all tuning conditions are met.
For the sake of comparison, a classical integer order PID controller has also been designed using the same phase margin specification. The results of the comparison between the classical PID control and fractional order control are depicted in the set of Figs. 7-10 . The test has been performed over a wide range of the operating interval of the DC motor, for various conditions (i.e. changing the magnetic brake at 50% and at 25%, affecting both gain and time constant values of the system). The FO-PID outperforms the PID controller in all cases, while exhibiting similar or less control effort. The comparison of the experimental results with the two designed robust Hinf and fractional order -controllers reveals the same output performances, as indicated in Fig. 11 and with similar control effort. Robustness tests have been performed at different speed values for setpoint, while maintaining the same brake position (i.e. at 25%), illustrated in Figs. 12 and 13 .
From the experimental validation it follows that FO-PID control outperforms classical PID control and has a similar performance as the Hinf control. The details are summarized in the following table.

Recommendations
As a suggested list of recommendations for industrial use of the robust methods presented in this paper, we propose the following: • If automatic tuning from industrial devices is available, use the standard PID tuning rules available at hand -however, take notice to retune the PID parameters periodically for avoiding performance deterioration; this has been addressed in [74] . • If robustness is not an issue, but rather fast control (not accurate, but in large tolerance intervals) then classical PID control may be used;  • If robustness is important, if available, use standard tools from control-related instrumentation industrial architecture. • If the above is not available, consider using control specific analytical tools, such as gain margin, phase margin, specifications based tuning rules. • For accurate robustness specification in frequency domain (e.g. bandwidth, etc) consider loop shaping tools, e.g. classical robust control design.   • If the system exhibits significant LPV dynamics, consider fractional order control for both optimization in terms of output performance and control effort (e.g. in loops where energy is costly, or the fuel consumption is related to varying operating conditions). • If the LPV dynamics are extreme, consider scheduling a set of fractional order controllers for stringent performance requirements.
These items are by no means exhaustive, but they offer a good support for the beginner user or plant operator in terms of tuning control loops for industry. Fig. 12. Step response of the modular servo system considering 100 rad reference angle. Fig. 13. Step response of the modular servo system considering 40 rad reference angle.

Conclusions
Essentially, this paper pleads for new emerging concepts in robust control design. We start by providing a revisited summary of the robust control methods at hand and most prevalent design methods for industrial process control. A mature robust control design method -Hinf, and a newly emerging one -fractional order PID control, are described.
An example of great simplicity and prevalence in mechatronic applications and basic loop control in manifold of production systems is used: the DC motor. A comparison is performed, illustrating the great potential of the fractional order control in robust design. The experimental results suggest that one can choose the fractional order PI λ D μ controller to maintain the same robustness of the control system as with the classical robust controllers, such as the Hinf controller. Both methods are frequency based design tools and their intrinsic complexity is similar. However, the fractional order control is a direct generalization of the classical PID control, a very popular control in industry, hence it holds the potential to become easier accepted by end-users than Hinf control.

Acknowledgments
This research was supported by the Janos Bolyai research fellowship of the Hungarian Academy of Sciences (Dulf).
This work was in part supported by a STSM Grant from COST action 15225.

Appendix A. On Fractional Order Systems and Controls
In this section we have selected a basic set of information to support the readers who aim to embrace the concepts of fractional order systems and control. The information is a summary from the following textbooks and articles: [21, 24, 25, 30, Podlubny1999] .

A1. Preliminaries
The fractional order system models are a generalization of the LTI system models. For instance, in the case of continuous time models, we have the transfer function given in the form: 0 a n s α n + a n−1 s α n−1 + · · · + a 0 s α 0 (A.1) For commensurate-order system, the continuous-time transfer function is given by:

A2. Stability
Stability analysis of fractional-order systems studies the solutions of the differential equations [75] . Alternatively, one may study the transfer function of the system (A.1) . Consider: a n s α n + a n−1 s α n−1 + · · · + a 0 s α 0 (A.3) with α i ∈ R + a multi-valued function of the complex variable s whose domain can be seen as a Riemann surface of a number of sheets. This is finite only in the case of ∀ i , α i ∈ Q + , where the principal sheet is defined by −π < arg(s) < π . For α i ∈ Q + , i.e. α = 1 /q, q a positive integer, the q sheets of the Riemann surface are determined by: s = | s| e jφ , (2k + 1) π < φ < (2k + 3) π, K = −1 , 0, . . . , q − 2 (A.4)  The case of k = −1 represents the principal sheet . For the mapping ω = s α , these sheets become regions in the plane ω defined by: The complete mapping is illustrated in Figs. A.1 and A.2 for the case of ω = s 1 / 3 . These three sheets correspond to: for π < arg(s) < 3 π second sheet 1 for 3 π < arg(s) < 5 π third sheet Hence, an equation of the type: a n s α n + a n−1 s α n−1 + · · · + a 0 s α 0 = 0 (A.6) which in general is not a polynomial, will have an infinite number of roots, among which only a finite number of roots will be on the principal sheet of the Riemann surface. It can be said that the roots which are in the secondary sheets of the Riemann surface are related to solutions that are always monotonically decreasing functions and only the roots that are in the principal sheet of the Riemann surface are responsible for different dynamics; e.g. i) damped oscillation, ii) oscillation of constant amplitude, or iii) oscillation of increasing amplitude with monotonic growth. A fractional-order system with an irrational-order transfer function G (s) = P(s) Q(s) is boundedinput bounded-output (BIBO) stable if and only if the following condition is fulfilled: The previous condition is satisfied if all the roots of Q(s) = 0 in the principal Riemann sheet have negative real parts. For the case of commensurate-order systems, whose characteristic equation is a polynomial of the complex variable λ = s α , the stability condition is expressed as: where λ i are the roots of the characteristic polynomial in λ. For the particular case of α = 1 , the well known stability condition for linear time-invariant systems of integer order is recovered: Consider a system defined by the transfer function: G (s) = 1 a n s nα + a n−1 s (n−1) α + · · · + a 1 s α + a 0 (A.10) where α = 1 q , with q, n ∈ Z + , a k ∈ R . Introducing the mapping λ = s α to obtain the function G ( λ) and applying the condition (A.8) , the stability of the system can be studied by evaluating the function G ( λ) along the curve Γ defined in the λ-plane: and illustrated in Fig. A.3 :

A3. Time and frequency domain response
The general equation for the response of a fractional-order system in the time domain can be determined by using the analytical methods. The response depends on the roots of the characteristic equation, having six different cases: 1. There are no roots in the Riemann principal sheet. In this case the response will be a monotonically decreasing function. 2. Roots are located in ( s ) < 0, (s) = 0. In this case the response will be as previous case. 3. Roots are located in ( s ) < 0, ( s ) = 0; the response will be a function with damped oscillations. 4. Roots are located in (s) = 0, ( s ) = 0; the response will be a function with oscillations of constant amplitude. 5. Roots are located in ( s ) > 0, ( s ) = 0; the response will be a function with oscillations of increasing amplitude 6. Roots are located in ( s ) > 0, (s) = 0; the response will be a monotonically increasing function.
In general, the frequency response has to be obtained by the direct evaluation of the irrational-order transfer function of the fractional-order system along the imaginary axis for s = jω, ω ∈ (0, ∞ ). However, for the commensurate order systems we can obtain the frequency response by the addition of the individual contributions of the terms of order α resulting from the factorization of the function: , For each of these terms, referred to as (s α + γ ) ±1 , the magnitude curve will have a slope from 0 to ± α20 dB/dec for higher frequencies, and the phase plot will go from 0 to ±α π 2 .

A4. Implementation aspects
A stable and efficient method for direct approximation of fractional order systems in the form of discrete-time rational transfer functions or non-rational transfer functions (NRTF), consists of four steps, summarized from [30] . The corresponding Matlab implementation is also given in the same reference.
Step 1) Discretize the fractional-order Laplace operator using a suitable generating function , that is equal to an interpolation between Euler and Tustin discretization rules: with α ∈ [0 ÷1] and T the sampling period. When α = 0, we obtain the classical Euler discretization rules, while when α = 1 we obtain the Tustin rules. The choice of the parameter α has a weighting effect on the frequency response, penalizing the errors on the magnitude or on the phase: larger value decreases the phase error near the Nyquist frequency, while a lower value ensures a lower magnitude error. Hence, this parameter may be tuned for a trade-off in performance at high frequencies among the magnitude and the phase of the process. The first step produces a discrete time fractional order system, G (z −1 ) , replacing s by the form in Eq. (A.19) , with given α weighting parameter and maximum frequency ω h values. The maximum sampling period is selected according to the Nyquist sampling theorem. The rational discrete-time approximation of the general fractional order system is determined in the frequency range ω ∈ (0, ω h ).
Step 2) Calculate the frequency response of the discrete-time fractional order system. The frequency response is computed based on the classical relation between the continuous time and discrete time domains, with z = e sT , with s = jω, where ω is a vector spaced equally in the frequency interval: with N s the total number of samples (the higher this value, the better is the approximation at lower frequencies). We obtain a vector of frequency response values of the fractional order discrete-time transfer function: Step 3) Calculate the impulse response of the discrete-time fractional order system. This step employs the inverse Fast Fourier Transform (FFT) algorithm, which converts the previously computed frequency domain response into a time domain response, at discrete instants [0, T , 2T , . . . , (N S − 1) T ] . The result of this step consists of a vector with N S impulse response value: N s nk n = 0, 1 , · · · , N S − 1 (A. 22) where G [ k ] is the frequency response of the original G NRTF ( s ).
Step 4) Determine a rational discrete-time transfer function that produces a similar impulse response as obtained from the inverse FFT. To determine the rational discrete-time transfer  where c i and d i are the coefficient determined on the desired order N of the resulted approximation, with i = 1 , 2, · · · , N . The accuracy of the approximation increases with N. Lower values for α improve the approximation of the magnitude curve, while higher values improve the approximation of the phase curve.

A5. Anti-Reset windup
Integrator windup or reset windup, refers to the closed loop situation where a large change in setpoint occurs (say, e.g. a positive change) and the integral term accumulates a significant error during the rise (windup). Often, this results in overshooting and monotonically increasing dynamics as this accumulated error is unbound (offset by errors in the other direction). The concept is illustrated in Fig. A.4 .
Integral windup occurs in presence of physical system limitations, i.e. process input saturation: the input of the process being limited at the top or bottom of its scale, making the error bound to a non-zero value. In this case, a common anti-windup solution is the integrator being turned off for periods of time until the response falls naturally back into an acceptable range [1] .  Among the additional functionalities that a FO-PID controller should possess, the antiwindup plays a major role. In fact, it is well known that the performance of any PID-type controller can be severely limited in practical cases by the presence of actuator saturation as depicted in Fig. A.5 .
Consider the control scheme of Fig. A.5 , where u is the controller output, u is the actual control effort, y is the process output, r is the set-point reference value, and e = r − y is the control error. The integrator windup mainly occurs when a step is applied to the reference set-point signal rather than to the manipulated variable (i.e. for a load disturbance). There are different anti-windup techniques that are typically employed for integer-order PID controllers, which can be also applied to fractional-order controllers. In particular, the best performance is obtained with the back-calculation technique.
Back-calculation consists of recomputing the integral term once the controller saturates. In particular, the integral value is reduced by feeding back the difference of the saturated and unsaturated control signal, as shown in Fig. A.6 where T t is an additional parameter called tracking time constant. The value of T t determines the rate at which the integral term is reset and its choice determines the performances of the total closed loop. The following definitions are used for the terms T t = T i for PI and T t = √ T i T d for PID. The same rules to calculate the tracking time are used for FO-PID. The following figure shows the back-calculation antiwindup scheme.
The ideal result is showed in Fig. A.7 , where the integral windup effect is adequately eliminated from the closed loop response.