A Robust Tuning Method for Fractional Order Pi Controllers

Proceedings of the 2nd IFAC Workshop on Fractional Differentiation and its Applications Porto, Portugal, July 19-21, 2006 A ROBUST TUNING METHOD FOR FRACTIONAL ORDER PI CONTROLLERS YangQuan Chen ∗,1 Huifang Dou ∗ Blas M. Vinagre ∗∗ Concha A. Monje ∗∗ ∗ Center for Self-Organizing and Intelligent Systems (CSOIS), Department of Electrical and Computer Engineering, Utah State University, Logan, UT 84322-4160, USA. ∗∗ Department of Electronic & Electromechanical Engineering, Industrial Engineering School, University of Extremadura, Avda. De, Elvas s/n, 06071-Badajoz, Spain Abstract: The application of fractional controller attracts more attention in the recent years. In this paper, a new tuning method for PIα controller design is proposed for a class of unknown, stable, and minimum phase plants. We are able to design a PIα controller to ensure that the phase Bode plot is flat, i.e., the phase derivative w.r.t. the frequency is zero, at a given gain crossover frequency so that the closed-loop system is robust to gain variations and the step responses exhibit an iso-damping property. Several relay feedback tests can be used to identify the plant gain and phase at the given frequency in an iterative way. The identified plant gain and phase at the desired tangent frequency are used to estimate the derivatives of amplitude and phase of the plant with respect to frequency at the same frequency point by Bode’s integral relationship. Then, these derivatives are used to design a PIα controller for slope adjustment of the Nyquist plot to achieve the robustness of the system to gain variations. No plant model is assumed during the PIα controller design. Only several relay tests are needed. Keywords: Fractional order controller, PIα controller tuning, relay feedback test, Bode’s integral, flat phase condition, iso-damping property. 1. INTRODUCTION In recent years, an increasing number of studies can be found related to the application of fractional controllers in many areas of science and engineering (Manabe, 1961; Oustaloup et al., 1995; Oustaloup et al., 1996; Raynaud and 1 Corresponding author: Dr. YangQuan Chen. E-mail: [email protected]; Tel. 01435-7970148; Fax: 01-435-7973054. URL: http://www.csois.usu.edu/people/yqchen. YangQuan Chen is supported in part by the TCO Bridging Fund of Utah State University (2005-2006). Blas M. Vinagre is partially supported by the Research Grant 2PR02A024 (Junta de Extremadura and FEDER). Zergaı̈noh, 2000; Podlubny, 1999; Vinagre and Chen, 2002). This is due to a better understanding of the fractional order calculus potentials revealed by many phenomena such as viscoelasticity and damping, chaos, diffusion and wave propagation. In theory, the control systems can include both the fractional order dynamic system to be controlled and the fractional-order controller. However, in control practice, it is more common to consider the fractional-order controller. This is due to the fact that the plant model may have already been obtained as an integer order model in the classical sense. In most cases, our objective is to apply the fractional-order control (FOC) to enhance the system control performance (Vinagre and Chen, 2002). For example, a generalization of the PID controller, namely the PIλ Dµ controller, involving an integrator of order λ and a differentiator of order µ where λ and µ can be real numbers, was proposed in (Podlubny, 1999), where the better response of this type of controller was demonstrated in comparison with the classical PID controller, when used for the control of fractional order systems. However, in general, there is no systematic way for setting the fractional orders λ and µ. α is derived to ensure that the slope of the Nyquist curve is equal to the phase of the open loop system at a given frequency. Section 3 presents an approximation for solving the fractional order α. The controller design procedure are summarized in Sec. 4. In Sec. 5, extensive illustrative simulations are given to demonstrate the effectiveness of the proposed design method. Finally, Sec. 6 concludes this paper with some remarks on further investigations. In this paper, we will concentrate on the fractional order PI controller, i.e., PIα controller 2. SLOPE ADJUSTMENT OF PHASE BODE PLOT 1 ), sα (1) In this section, we will show how Ki and α are related under the new condition (2). where α is a real number and α ∈ (0, 2). For the systematic design of α, a new tuning condition, called flat phase condition, first proposed in (Chen et al., 2003), will be used which can give a relationship between Ki and α. Specifically, in addition to the gain and phase margin specifications, we propose to add an additional condition that the phase Bode plot at a specified frequency wc where the sensitivity circle tangentially touches the Nyquist curve is locally flat. When achieved, this new condition will make the system more robust to gain variations. This additional condition can be 6 G(s) |s=jwc = 0 with its equivalent expressed as d ds expression given as following: Substitute s by jw in close loop system (3) so that the close loop system can be written as G(jw) = C(jw)P (jw), where C(s) = Kp (1 + Ki 6 dG(s) |s=jwc = 6 G(s)|s=jwc , ds (2) where wc is the frequency at the point of tangency where the sensitivity circle tangentially touches the Nyquist curve. In (2), G(s) = C(s)P (s) (3) is the transfer function of the open loop system including the controller C(s) and the plant P (s). At the first look of (2), it seems complicated since the derivative of the phase of the system at wc has to be known. Fortunately, Bode’s integrals (Bode, 1945) can be used to approximate the derivatives of the amplitude and the phase of a system with respect to frequency at a given frequency. To obtain the approximate derivatives, the knowledge of the amplitude and the phase of the system at the given frequency together with the static gain of the system. In practice, wc can be set as the gain crossover frequency. Our objective in this paper is to devise a way to retrieve the parameters Kp , Ki and α of the controller C(s) to ensure the flat phase condition (2). Then, we can adjust Kp to make the sensitivity circle exactly tangentially touches the Nyquist curve on the flat phase. The remaining part of this paper is organized as follows. In Sec. 2, the relationship between Ki and C(jw) = Kp (1 + Ki = Kp [(1 + 1 ) (jw)α Ki απ απ Ki cos ) − j α sin ] α w 2 w 2 (4) is the PIα controller obtained from (1). The phase of the closed-loop system is given by G(jw) = 6 C(jw) + 6 P (jw) 6 = φ0 + tan−1 [ wα+1 sin (α+1)π + Ki w 2 wα+1 cos (α+1)π 2 ]− (α + 1)π .(5) 2 where φ0 = 6 P (jw). Then, the derivative of the closed-loop system G(jw) with respect to w can be written as follows: dG(jw) dC(jw) dP (jw) = P (jw) + C(jw) . (6) dw dw dw From (2), the phase of the derivative of the open loop system should be known in advance which obviously can not be obtained directly from (6). So, we need to simplify (6). Consider (6). The PIα controller C(jw) is given by (4) whose derivative with respect to w is that jαKp Ki dC(jw) . =− dw (jw)α+1 For calculation of dP (jw) dw , (7) we have lnP (jw) = ln|P (jw)| + j 6 P (jw). (8) Differentiating (8) with respect to w gives dlnP (jw) 1 dP (jw) = dw P (jw) dw = d6 P (jw) dln|P (jw)| +j . dw dw (9) Straightforwardly, we arrive at 4 dP (jw) dln|P (jw)| d6 P (jw) = P (jw)[ +j ]. (10) dw dw dw 2 0 −2 Substituting (4), (7) and (10) into (6) gives −4 −6 sa Ki sp dG(jw) )( + j ) = Kp P (jw)[(1 + dw (jw)α w w jαKi ], (11) − (jw)α+1 where sa (w) and sp (w), first introduced in (Karimi et al., 2002), are defined as follows: dln|P (jw)| sa (w) = w , dw sp (w) = w d6 P (jw) . dw (12) −12 −14 −16 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Fig. 1. ∆ vs. α. The illustration for αzero with different sp : -0.5, -1, -1.5, -2 origin approximation (13) + sa Ki + sp wα cos (α+1)π − αKi sa wα sin (α+1)π 2 2 1.6 1.4 dG(jw) (α + 1)π = φ0 − + tan−1 [ dw 2 − sp wα sin (α+1)π − sp K i sa wα cos (α+1)π 2 2 −10 1.8 Hence, the slope of the Nyquist curve at any specific frequency w is given by 6 −8 1.2 1 0.8 0.6 0.4 ].(14) Following the condition (2), one obtains the relationship between Ki and α as follows: 0.2 0 −3 −2.5 −2 −1.5 −1 −0.5 0 Fig. 2. Comparison of the original αzero and sp relationship and αzero (sp ) wα (α + 1)π (α + 1)π √ [αcos − 2sp sin + ∆],(15) are performed. However, from (15), we notice that 2sp 2 2 ∆ can not be negative so as to ensure that Ki is a real number. Therefore, for each sp , α must be where ∆ = α2 cos2 (α+1)π −4αsp cos (α+1)π sin (α+1)π + limited within a more restricted interval. In Fig. 1, 2 2 2 (α+1)π it is shown that for different sp , there exists an 2 2 4s2p sin − 4s . It should be mentioned that p 2 α zero , such that, when α > αzero , ∆ > 0. Clearly, due to the nature of the quadratic equation, an alα there exist a relationship between αzero and sp . (α+1)π w [αcos 2 − ternative relationship, i.e., Ki = 2s p Using least squares fitting techniques, αzero (sp ) √ can be approximately expressed by 2sp sin (α+1)π − ∆], has been discarded. Also 2 noted is that in (15) only sp presents. So, there 2.0093s2p − 0.5211sp + 0.0035 is no need to compute sa . αzero = . (17) s2p − 0.9359sp + 0.0474 The approximation of sp can be given as follows (Karimi et al., 2002): Figure 4 shows that the accuracy of the approxd6 P (jw) imate function (17) to the actual αzero and sp |w0 sp (w0 ) = w0 relationship is practically acceptable. dw 2 ≈ 6 P (jw0 ) + [ln|Kg | − ln|P (jw0 )|] (16) π 3. PHASE MARGIN ADJUSTMENT Ki = where |Kg | = P (0) is the static gain of the plant, 6 P (jw0 ) is the phase and |P (jw0 )| is the gain of the plant at the specific frequency w0 . For the systems containing an integrator, because of the phase of the integrator is constant and its derivative is zero, sp should be estimated by using the partial model of the system without the integrator. Note that, the pure time delay has no effect on the estimation of sp . For most of the plants, sp can be selected between -3 and 0. In general, sp depends on the system dynamics and the frequency at which the simulations or experiments To determine all the three parameters for PIα controller, we have already established the relationship (15) in the previous section. However, we still need two other equations. Assume that the phase of the open loop system at the gain crossover frequency wc is 6 G(s)|s=jwc = φ0 + tan−1 [ − (α + 1)π . 2 (α+1)π + Ki wc 2 ] (α+1)π α+1 wc cos 2 wcα+1 sin (18) The corresponding gain is r • i) Given wc , the gain crossover frequency; • ii) Given Φm , the desired phase margin; Ki απ 2 Ki απ 2 ) + ( α sin ) = 1.(19)• iii) From the real plant, obtain the measure|G(jwc )| = Kp |P (jwc )| (1 + α cos wc 2 wc 2 ments of 6 P (jwc ) and |P (jwc )| using the iterative relay tests proposed in (Chen et Denote Φm the desired phase margin, i.e., 6 G(s)|s=jwc = al., 2003); Φm − π. Straightforwardly, we have • iv) Calculate an estimation of sp (wc ) according to (16); Φcontroller = Φm − π − φ0 • v) Compute α and Ki from (22) and (15), respectively; wα+1 sin (α+1)π + Ki wc (α + 1)π ]− .(20) = tan−1 [ c α+1 2 (α+1)π • vi) Obtain Kp from (19). 2 wc cos 2 However, it is very complex to solve (20) together with (15) to get α, Kp and Ki . However, from an important observation that by substituting (15) into (20), Φcontroller is the function only of sp and α, not explicitly of w any more, we can proceed to use the LS fitting again to approximate the function in (20). We propose to use the following form of approximation: Φcontroller ≈ A(sp )α2 + B(sp )α + C(sp ) , α2 + D(sp )α + E(sp ) (21) Remark 4.1. Due to the constraint in αzero (sp ), wc should not be chosen too aggressively. As usual, Φm should be selected from 30◦ to 60◦ . 5. ILLUSTRATIVE EXAMPLES The PIα controller design method presented above will be illustrated via simulation examples. In the simulation, the following plants, studied in (Wallén et al., 2002), will be used. α ∈ (αzero (sp ), 2), 1 , n = 1, 2, 3, 4; (s + 1)(n+3) (23) 1 ; s(s + 1)3 (24) P6 (s) = 1 e−s ; (s + 1)3 (25) P7 (s) = 1 e−s ; s(s + 1)3 (26) 1 e−s ; (s + 1) (27) Pn (s) = where A(sp ), B(sp ), C(sp ), D(sp ) and E(sp ) are polynomial functions of sp . Our fitting results are summarized below for completeness: A(sp ) = −0.00652s7p − 0.07259s6p − 0.32682s5p − P5 (s) = 0.7568s4p −0.92446s3p − 0.44551s2p + 0.19469sp + 0.00283, B(sp ) = 0.0273s7p + +4.57371s3p 0.30814s6p + + 3.04877s2p 1.41817s5p + 3.42016s4p P8 (s) = + 0.30284sp − 0.01085, C(sp ) = −0.02871s7p − 0.32823s6p − 1.54191s5p − 3.85236s4p −5.52107s3p − 4.39267s2p − 1.42674sp + 0.01003, D(sp ) = 0.02154s7p + 0.2571s6p + 1.26183s5p + 3.3037s4p +5.04888s3p + 4.74463s2p + 3.03777sp − 2.09475, E(sp ) = −0.02433s7p − 0.29619s6p − 1.49144s5p − 4.05076s4p −6.55861s3p − 6.81121s2p − 5.17001sp + 0.10642. So, α can be sovled from the following approximate relationship: A(sp )α2 + B(sp )α + C(sp ) = Φm − π − φ0 (22) α2 + D(sp )α + E(sp ) 5.1 General plants Pn (s) Let us consider the following fifth order plant first, i.e., P2 (s), which is also discussed in (Karimi et al., 2002). The specifications are set as wc =0.295 rad./sec. and Φm = 45◦ . The PIα controller designed by using the proposed tuning formulae is C2α (s) = 1.378(1 + s0.168 1.383 ). The Bode and the Nyquist diagrams are compared in Fig. 3. 1 200 0.5 0 −200 0 −400 Clearly, given sp , it is much easier to obtain α by solving (22) than by solving (20). −0.5 −600 −3 10 −2 10 −1 10 0 10 1 10 2 10 3 10 −1 −100 −1.5 −200 −2 Remark 3.1. For α ∈ (αzero (sp ), 2) and sp ∈ (−3, 0), the precision of the estimation is found to be acceptable via our extensive numerical experiments. −300 −2.5 −400 −500 −3 10 −2 10 −1 10 0 10 1 10 (a) Bode plot 2 10 3 10 −3 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 (b) Nyquist plot Fig. 3. Bode and Nyquist plots for C2α (s)P2 (s). 4. THE PIα CONTROLLER DESIGN The procedures to determine the PIα controller parameters are briefly summarized in the following: From the Bode diagram in Fig. 3(a), it is seen that the phase curve near the gain crossover frequency is flat due to the proposed design method. The phase margin exactly equals 45◦ . That means the controller moves the point P (0.295j) of the Step Response 1.6 200 0 1 For the fourth order plant: P1 (s) = (s+1) 4 , the 0.2512 proposed controller is 0.695(1 + s1.369 ) with respect to β=0.5, wc =0.374 rad./sec. and Φm =45◦ . The controller designed by the modified Ziegler1 ). The results are Nichols method is 0.062(1 + 0.22s summarized in Fig. 5. Nyquist Diagram Step Response 1 1.6 0.5 1.4 0 1.2 1 Amplitude −0.5 −1 0.8 −1.5 0.6 −2 0.4 −2.5 −3 −3 0.2 −2.5 −2 −1.5 −1 −0.5 0 0.5 0 1 0 20 40 60 Real Axis 80 100 120 140 160 Time (sec) (a) Nyquist plots (b) Step responses Fig. 5. Comparisons of Bode plots and step responses for P1 (s) (Bode plots - Dashed line: The modified Ziegler-Nichols, Solid line: The proposed; Step responses - Solid line: The modified proposed controller with gain variations 1, 1.1, 1.3; Dotted line: The modified Ziegler-Nichols controller with the same gain variations 1, 1.1, 1.3) 1.4 −200 1.2 −400 1 −3 10 −2 10 −1 10 0 10 1 10 2 10 3 10 Amplitude −600 0.8 0 0.6 −100 0.4 −200 −300 0.2 −400 0 −500 The other plants shown in (23) have similar simulation results. We briefly summarized the results as follows for further illustrations: Imaginary Axis Nyquist curve to a point of C(jw)P (jw) on the unit circle having a phase of −135◦ and at the same time makes the Nyquist curve match the constraint of (2). A bad new is, from Fig. 3(b), the Nyquist curve of the open loop system is not tangential to the sensitivity circle at the flat phase region. But if we are allowed to adjust the open loop gain, we can shift the gain Bode plot get a different gain crossover frequency. Define the frequency interval corresponding to the flat phase is [wl , wh ]. So, the gain crossover frequency wc can be moved in [wl , wh ] by adjusting Kp by wl wh , wc ]. In this Kp′ = βKp where β belongs to [ w c case, setting β = 0.5 gives the modified proposed ′ controller C2α (s) = 0.689(1 + s0.168 1.383 ). For comparison, the PI controller designed by the modified Ziegler-Nichols method is C2mZN (s) = 0.344(1 + 1 1.237s ). The Bode plots are compared in Fig. 4(a). The step responses of the close loop system are compared in Fig. 4(b). Comparing the closed-loop system with the proposed modified controller to that with the modified Ziegler-Nichols controller, the overshoots of the step response from the proposed scheme remain invariant under gain variations. However, the overshoots of the modified Ziegler-Nichols controller change remarkably. −3 10 −2 10 −1 10 0 10 1 10 (a) Bode plots 2 10 3 10 0 20 40 60 80 100 120 140 160 Time (sec) 1 For the sixth order plant: P3 (s) = (s+1) 6 , the 0.132 proposed controller is 0.526(1 + s1.385 ) with respect to β=0.4, wc =0.242 rad./sec. and Φm =45◦ . The controller designed by the modified Ziegler1 Nichols method is 0.289(1 + 1.327s ). The results are summarized in Fig. 6. (b) Step responses Nyquist Diagram In practice, the fractional order integrator in the proposed PIα controller can not be exactly achieved since it is an infinite dimensional filter. A band-limit implementation of the fractional order integrator is important in practice, i.e., the finite-dimensional approximation of the fractional order system should be done in a proper range of frequencies of practical interest (Chen and Moore, 2002). The approximation method we use in this paper is the Oustaloup Recursive Algorithm (Oustaloup et al., 2000). In our simulations, for approximation of the fractional order integrator, the frequency range of practical interest is selected to be from 0.001Hz to 1000Hz. The sampling time and the number of the recursive zero-pole pairs are assigned as 0.001 sec and 13, respectively. 1.6 1.4 0 1.2 1 Amplitude −0.5 Imaginary Axis Fig. 4. Comparisons of Bode plots and step responses (Bode plots - Dashed line: The modified Ziegler-Nichols C2mZN (s)P2 (s), Solid ′ line: The proposed C2α (s))P2 (s); Step responses - Solid line: The modified proposed controller with gain variations 1, 1.1, 1.3; Dotted line: The modified Ziegler-Nichols controller with the same gain variations 1, 1.1, 1.3) Step Response 1 0.5 −1 0.8 −1.5 0.6 −2 0.4 −2.5 −3 −3 0.2 −2.5 −2 −1.5 −1 −0.5 0 Real Axis (a) Nyquist plots 0.5 1 0 0 20 40 60 80 100 120 140 160 Time (sec) (b) Step responses Fig. 6. Comparisons of Bode plots and step responses for P3 (s) (Bode plots - Dashed line: The modified Ziegler-Nichols, Solid line: The proposed; Step responses - Solid line: The modified proposed controller with gain variations 1, 1.1, 1.3; Dotted line: The modified Ziegler-Nichols controller with the same gain variations 1, 1.1, 1.3) 1 For the seventh order plant: P4 (s) = (s+1) 7 , the 0.105 proposed controller is 0.516(1 + s1.389 ) with respect to β=0.4, wc =0.206 rad./sec. and Φm =45◦ . The controller designed by the modified Ziegler1 ). The results Nichols method is 0.164(1 + 0.949s are summarized in Fig. 7. From these general plant class Pn (s), the effectiveness of the proposed PIα controller is clearly demonstrated. Nyquist Diagram 1.6 1.4 0 1.2 1 Amplitude Imaginary Axis −0.5 −1 0.8 −1.5 0.6 −2 0.4 −2.5 0.2 −3 −3 phase, open loop unstable systems. We are particularly interested in fractional-order PI control of biomimetic systems. Step Response 1 0.5 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 0 7. ACKNOWLEDGMENTS 0 20 40 Real Axis 60 80 100 120 140 160 Time (sec) (a) Nyquist plots (b) Step responses Fig. 7. Comparisons of Bode plots and step responses for P4 (s) (Bode plots - Dashed line: The modified Ziegler-Nichols, Solid line: The proposed; Step responses - Solid line: The modified proposed controller with gain variations 1, 1.1, 1.3; Dotted line: The modified Ziegler-Nichols controller with the same gain variations 1, 1.1, 1.3) 5.2 Plant with an integrator P5 (s) (Omitted due to space limit) 5.3 Plant with a time delay P6 (s) (Omitted due to space limit) 5.4 Plant with an integrator and a time delay P7 (s) (Omitted due to space limit) 5.5 First order plus time-delay plant P8 (s) (Omitted due to space limit) 6. CONCLUSION A new PIα tuning method is proposed for a class of unknown plants in this paper. Given the gain crossover frequency, the phase margin and with an additional condition that the phase Bode plot at the specified frequency is locally flat, we can design the PIα controller to ensure that the closedloop system is robust to gain variations and to ensure that the step responses exhibit an isodamping property. Comparing with the “flat phase” PID controller proposed in (Chen et al., 2003), PIα , although also only having three parameters, can achieve very good performance. Most importantly, PIα can be easily applied for the first order system while in (Chen et al., 2003), the FOPTD plants cannot be well handled. This makes PIα more advantages in practice because every system can be approximated by the first order plus a time delay model. Further investigations include the experimental verification and exploration of nonminimum The authors are grateful to Professor Li-Chen Fu, Editorin-Chief of Asian Journal of Control for providing a complimentary copy of the “Special Issue on Advances in PID Control”, Asian J. of Control (vol. 4, no. 4). The simulation study was helped by C. H. Hu. REFERENCES Bode, H. W. (1945). Network Analysis and Feedback Amplifier Design. Van Nostrand. New York. Chen, Y. Q. and K. L. Moore (2002). Discretization schemes for fractional-order differentiators and integrators. IEEE Trans. Circuits Syst. I 49, 363–367. Chen, Y. Q., C. H. Hu and K. L. Moore (2003). Relay feedback tuning of robust PID controllers with iso-damping property. In: Proceedings of The 42nd IEEE Conference on Decision and Control. Hawaii. Karimi, A., D. Garcia and R. Longchamp (2002). PID controller design using Bode’s integrals. In: Proceedings of the American Control Conference. Anchorage, AK. pp. 5007–5012. Manabe, S. (1961). The non-integer integral and its application to control systems. ETJ of Japan 6(3-4), 83–87. Oustaloup, A., B. Mathieu and P. Lanusse (1995). The crone control of resonant plants: Application to a flexible transmission. Eur. J. Contr. Oustaloup, A., F. Levron, F. Nanot and B. Mathieu (2000). Frequency band complex non integer differentiator: Characterization and synthesis. IEEE Trans. Circuits Syst. I 47, 25– 40. Oustaloup, A., X. Moreau and M. Nouillant (1996). The crone suspension. Control Eng. Pract. 4(8), 1101–1108. Podlubny, Igor (1999). Fractional-order systems and PIλ Dµ -controllers. IEEE Trans. Automatic Control 44(1), 208–214. Raynaud, H. F. and A. Zergaı̈noh (2000). Statespace representation for fractional order controllers. Automatica 36(7), 1017–1021. Vinagre, Blas M. and YangQuan Chen (2002). Lecture note on fractional calculus applications in automatic control and robotics. In: The 41st IEEE CDC2002 Tutorial Workshop # 2 (Blas M. Vinagre and YangQuan Chen, Eds.). Las Vegas, Nevada, USA. pp. 1–310. [Online] http://mechatronics.ece.usu.edu /foc/cdc02 tw2 ln.pdf. Wallén, A., K. J. Åström and T. Hägglund (2002). Loop-shaping design of PID controllers with constant ti /td ratio. Asian Journal of Control 4(4), 403–409.