*Department of Mechanical Engineering, Southern Taiwan University of Science and Technology, Tainan, Taiwan, E-mail*: huangt@stust.edu.tw

## 10.1 INTRODUCTION

The governing equations of dynamic systems or processes can be generalized by arbitrary or fractional-order differential equations including integer-order ones, therefore it would be suitable to apply fractional calculus to controller design. The concept of fractional calculus is first mentioned in a letter Leibniz wrote to de L’Hopital in 1695. And its applications may be dated back to the early 19^{th} century, Abel uses fractional calculus to formulate the “tautochrone” problem and find its solution [1]. Later on, researchers such as Liouville and Heaviside also use fractional calculus to derive mathematical models for systems including potential theory, heat equation, and notch design of a dam, and solve those problems therewith. Then the related research seemed to have stalled inexplicably until the late 20^{th} century. Researches on miscellaneous fractional-order systems, including rheology, viscoelasticity, diffusion, electromagnetic theory, electrochemistry of corrosion, and statistics have been made during the past several decades. Inevitably, researchers would apply the technique to control systems.

The well-known PID controller is typical of controllers used in the industry. It uses the linear combination of the input error value, and it’s integral and derivative as the output. One advantage of PID control is that no complex algorithm or calculation is involved. In addition, there are empirical parameters tuning methods designed for it, such as the Ziegler-Nichols method [2] which requires only the system’s time response without knowing its parametric model. The simplicity and effectiveness of these tuning methods make the PID controller even more appealing. It is desired to introduce the fractional calculus technique to such a controller and find a suitable tuning rule for the adapted controller, namely the fraction order PID controller.

Most researchers propose to optimize fractional-order PID controllers’ parameters by using gain crossover frequency and phase margin [3]. Yet, this method requires fairly good initial guess so as to have the parameters converge to optimum. Valério and Sâ da Costa [4] propose a scheme for tuning fractional-orders PID controllers’ parameters adhering to the concept of the Ziegler-Nichols method. However, it’s useful only when step response with time delay in sigmoidal shape, not suitable for system response without time delay.

This research studied the feasibility of applying the fractional-order PID controller (a.k.a. PI^{λ}Dμ controller, *A, μ* > 0) to a DC motor-driven two-link rotary pendulum system, an underactuated system. Some researchers also apply the fractional-order PID controller to underactuated systems such as an overhead 2D crane and an inverted pendulum on a cart [5, 6] without any specific approach. In view of the simplicity and effectiveness of the Ziegler-Nichols method, a variant of the Ziegler-Nichols method is proposed to determine the parameters of the fractional-order PID controller in this research. Different order PI^{λ}Dμ controllers are tried for comparison in order to find a PI^{λ}Dμ controller of appropriate order for the specified system.

This chapter consists of the following sections: Section 10.1 is the introduction, Section 10.2 briefly introduces the fractional calculus, fractional-order PID controllers, and the proposed tuning method, Section 10.3 presents the modeling of a DC motor-driven two-link rotary pendulum system, Section 10.4 estimates the controller parameters by using local linearization and the Ruth-Hurwitz stability criterion, Section 10.5 gives the simulation results, and Section 10.6 is the conclusion.

## 10.2 FRACTIONAL-ORDER PID CONTROLLERS

*10.2.1 FRACTIONAL CALCULUS*

Fractional calculus operator is defined as follows [7]:

(10.1)

where *a* and *t* are the lower limit and upper limit of the integral, respectively, and *a* is the fractional order, which will be taken as a real number afterwards. If α > 0, * _{a} D^{α}_{t}* is a fractional-order differential, and if α < 0,

*is a fractional-order integral. Note that when the infinitesimal increment*

_{a}D^{α}_{t}*h*in time equals

*t–a*, both the left and right subscripts of the

*D*operator can be omitted, and when

*a*is an integer, the operator functions as the commonly seen integer

*D*operator, that is,

*D*/

^{k}= d^{k}*dt*(

^{k}*k*∈

*N*) as a

*k*-th order differentiator, and

*D*

^{-1}= ∫dτ for the first order integral and

*D*

^{–}

*=∫∫…∫*

^{k}*dτ*..

*dτdτ*(

*k*times

*k*∈

*N*) for the

*k*-th order integral.

One most used fractional-order differential definition is the Grünwald-Letnikov definition [7], which is expressed as follows:

(10.2)

where

(10.3)

with

Meanwhile a famous fractional-order integral definition is the Riemann-Liouville definition [7], which is expressed as follows:

(10.4)

where

(10.5)

is the well-known Euler’s gamma function.

As engineering problems described by ordinary differential equations in integer order can be solved using the Laplace transform, so are those systems in fractional-order differo-integral equations.

According to the Riemann-Liouville definition, the Laplace transform for fractional-order differo-integrals is:

(10.6)

when (*n–1*) ≤ α < *n* If _{0}*D _{t} ^{α–j–1} f* (0) = 0,

*j=*0, 1, 2,…,

*n–1*, then

*ℒ{*

_{0}

*D*(

^{α}_{t}f*t*)}=

*s*(

^{a}F*s*)

Now the order of the Laplace operator *s* is a fraction in fractional-order differo-integral.

*10.2.2 PID CONTROLLERS AND ZIEGLER-NICHOLS METHOD*

*10.2.2.1 Pid Controllers*

PID controllers have been popularly used after its debut. Its control effort *u*(*t*) in time domain can be expressed as:

(10.7)

where *K _{p}, K_{l}* and

*K*are the proportional gain, the integral gain, and the derivative gain, respectively. And the corresponding transfer function obtained by taking the Laplace transform of Eq. (10.7) is:

_{D}(10.8)

Alternatively, the parameter set is defined somehow in a different way as follows:

(10.9)

And the corresponding transfer function of Eq. (10.9) is:

(10.10)

By comparison of these two expressions, it is easy to conclude that *K _{1} = K_{p}/T_{1}* , and

*K*. Furthermore, it’s noteworthy that the 3 terms inside the parenthesis of Eq. (10.9) have the same dimension of

_{D}= K_{P}T_{D}*e*(

*t*), where the unit of

*T*and

_{1}*T*(both of physical quantity “time”) is “second” used to cancel the unit of “

_{D}*dt*” in the integral and derivative, respectively.

*10.2.2.2 Ziegler-Nichols Method*

The Ziegler-Nichols method [2] is a heuristic method for tuning PID controllers. It is meant to make the closed-loop control system’s response acceptable though not optimal. The procedure of the Ziegler-Nichols method is first to set the integral gain and derivative gain to zero which results in a P controller. Then vary the proportional gain to yield a sustained periodic oscillation in the output (or close to it) by trial and error. Mark the critical gain as the ultimate gain *K _{u}* and the corresponding period as the ultimate period

*T*. Then, according to the Ziegler-Nichols method, the empirical PID parameters are assigned to be

_{u}*K*0.6

_{p}=*K*,

_{u}*T*, and

_{I}= T_{u}/2*T*8, which is equivalent to have the alternative parameter set:

_{D=}T_{u}/*K*0.6

_{p}=*K*,

_{u}*K*1.2

_{I}=*K*, and

_{u}/T_{u}*KD = 3K*40.

_{u}T_{u}/*10.2.3 FRACTIONAL-ORDER PID CONTROLLERS AND TUNING RULE*

*10.2.3.1 FRACTIONAL-ORDER PID CONTROLLERS*

Oustalop [8] proposes to apply fractional-order controllers to dynamic systems, and names such robust controllers as “Commande Robuste d’Ordre Non-Entier” (CRONE). He also introduces the fractional-order PID controller, or the PI^{λ}D^{μ} controller (*λ* and μ are non-negative fractions). Similar to the PID controllers, the control effort of the fractional-order PID controller *u*(*t*) in the time domain can be expressed as:

(10.11)

where *K _{p}, K_{I,}*

_{X}, and

*K*are the proportional gain, the fractional-order integral gain, and the fractional-order derivative gain, respectively.

_{D,μ}And the corresponding transfer function takes the following formml:

(10.12)

when λ= μ = 1 in Eq. (10.12), G_{C}(*s*) is an integer PID controller; *λ* = 1, *μ* = 0, G_{C}(*s*) an integer PD controller; *λ* = 0, *μ* = 1, G_{C}(*s*) an integer PI controller; and *λ* = *μ* = 0, G_{C}(*s*) an integer P controller. Without limiting *λ* and μ to be 0 or 1, the fractional-order differo-integral offers more flexibility than the integer-order differo-integral in the controller design.

*10.2.3.2 TUNING RULE: GENERALIZED ZIEGLER-NICHOLS METHOD*

It is desired to find a general tuning rule, just like the Ziegler-Nichols method, of setting empirical values for fractional-order PID controllers’ parameters. As mentioned in Section 10.2.2.1, the three terms inside the parenthesis of Eq. (10.9) are of the same dimension, i.e., the unit of *T _{I}* and

*T*cancels the unit of “

_{D}*dt*” in the integral and derivative, respectively. Keeping this in mind, now the orders of the fractional integral and derivative are

*λ*and μ in Eq. (10.11), respectively. In order to cancel the unit of “(

*dt*)

*” in the integral and “(*

^{λ}*dt*)

*μ*“ in the derivative, (

*T*) and (

_{I}*T*)

_{D}*can be employed in the control law for the fractional integral and derivative, respectively. Therefore, the control law for the fractional PID controller can be expressed as:*

^{μ}(10.13)

which yields

Note that λ *=* μ *=* 1 is a special case of the fractional-order PID controller, and indeed it becomes the commonly seen PID controller with parameters set by the Ziegler-Nichols method. Hence Eq. (10.13) generalizes the integer-order PID controller where λ *=* μ *=* 1.

*10.2.4 IMPLEMENTATION OF FRACTIONAL-ORDER PID CONTROLLERS*

*10.2.4.1 APPROXIMATION OF FRACTIONAL-ORDER LAPLACE OPERATOR*

The fractional-order differentiator *s ^{r}* (

*r*is a positive fraction) in continuous-time domain can be expressed in terms of

*z*

^{–1}in discrete-time domain by a generating function

*s*=

*w*(

*z*

^{–1}) just like the integer-order differentiator

*s*(

*z*is an operator in Z transform); similarly it applies to the fractional-order integrator

*s*(

^{–r}*r*is a positive fraction). Adopting the Tustin generating function, the fractional differo-integral would be discretized as [9]:

(10.14)

where can be expanded in Taylor series as follows:

(10.15)

Therefore,

(10.16)

and

(10.17)

where *A _{n}* (

*z*

^{–1},

*r*), and likewise

*A*(

_{n}*z*

^{–1},

*–r*), can be found by the following iterations:

(10.18)

with the coefficient:

Without doubt the accuracy of the approximation is determined by the selection of *n*.

*10.2.4.2 IMPLEMENTATION OF FRACTIONAL-ORDER PID CONTROLLERS*

By substituting the above discrete-time approximation of the fractional-order Laplace operator *s ^{r}* and

*s*into Eq. (10.11), the fractional-order PID controller in Eq. (10.10) can be approximated in discrete-time domain as:

^{–r}(10.19)

where the index *k* stands for the *k*-th sampling time, e.g., *e*(*k*)*= e*(*t* = *kT*_{s}) with sampling time *T _{s}*.

## 10.3 AN UNDERACTUATED SYSTEMML: A ROTARY PENDULUM SYSTEM

A rotary pendulum system (Figure 10.1) which uses a permanent magnet DC motor to drive a two-link rotary pendulum is the underactuated system employed to demonstrate the performance of the fractional-order PID controller in this research. The rotary pendulum system is an underactuated system because there is only one actuator to drive link 1 but no actuator to drive link 2. Both links have one rotational degree of freedom (d.o.f.) around different axes. The motion of link 2 is coupled with link 1 so it has to be considered when controlling link 1.

*10.3.1 DC MOTOR MODEL*

A DC motor model is used to drive the two-link rotary pendulum. The armature voltage is *e _{a}*, and the current through the armature is

*i*. The armature current generates a motor torque

*τ*, where

_{M}= K_{T}i*K*is the torque constant of the DC motor. Assume that the angular displacement of the motor shaft is

_{T}*θ*and the load is

*τ*, then the mechanical dynamics of the DC motor is:

_{L}(10.20)

where *J _{M}* and

*B*are the rotational inertia and the rotational damping coefficient of the motor, respectively.

_{M}As the motor shaft rotates at a speed *ω = dθ/dt*, the so-called back e.m.f. is induced by the relationship *e _{b.e.m.f.} = K_{E}*ω where

*K*is the back e.m.f. constant of the DC motor and

_{E}*K*=

_{E}*K*in the SI unit system. Then applying the Kirchhoff’s voltage law to the armature gives the electrical dynamics of the DC motor below:

_{T}(10.21)

where *L _{a}* and

*R*are the inductance and resistance of the armature, respectively.

_{a}*10.3.2 TWO-LINK ROTARY PENDULUM MODEL*

The mechanical part of the rotary pendulum, driven by the torque τ* _{L}*, consists of two rotational links of 1 d.o.f. each: link 1 and link 2 as shown in Figure 10.1. Link 1 is connected to the motor shaft and rotates around the center O, and link 2 is connected to the other end of link 1 and rotates around point P in a plane normal to link 1. The masses of link 1 and 2 are

*m*

_{1}and

*m*

_{2}, respectively. The lengths of link 1 and 2 are

*l*

_{1}and

*l*

_{2}, respectively. And the moment of inertia of link 1 around O is

*J*

_{1}. Assume that θ

_{1}and θ

_{2}are the angular displacements of link 1 and 2, respectively. Then the angular velocities of link 1 and 2 are

*θ̇̇*

_{1}and

*θ̇̇*

_{2}respectively. Since the center O is on the motor shaft, it is straightforward

*θ*

_{1}=

*θ*.

The dynamic model of the two-link rotary pendulum can be derived using the Euler-Lagrange equation [10]:

(10.22)

where the Lagrangian *L* is the difference of the kinetic energy *T* and potential energy *V* of the system (*L = T – V*), *q _{i}* is the

*i*-th generalized coordinate, and

*Q*is the

_{nc,i}*i*-th non-conservative generalized force. The kinetic energy is:

(10.23)

where only the first term accounts for link 1, and the other 4 terms are from link 2. Meanwhile, the potential energy only comes from link 2, which is:

(10.24)

Therefore the Lagrangian is:

(10.25)

where *θ*_{1} and *θ̱̱*_{2}are chosen as the 2 generalized coordinates, then the corresponding Euler-Lagrange equations are:

(10.26)

and

(10.27)

The dynamic equation of the mechanical subsystem then can be expressed as follows:

(10.28)

(10.29)

where

*10.3.3 COMPLETE SYSTEM MODEL*

The complete system consists of a DC motor represented by Eq. (10.20) and (10.21), and a two-link rotary pendulum represented by Eq. (10.28) and (10.29). Since Eq. (10.28) gives the term *τ _{L}*, it can be substituted into Eq. (10.20), and using the fact that θ

_{1}= θ, then the whole system therewith can be represented by the dynamic model below.

(10.30)

(10.31)

(10.32)

## 10.4 ESTIMATION OF CONTROLLER PARAMETERS

The aforementioned fractional-order Laplace operator approximation is used to design the fractional-order PID controller for the rotary pendulum system, and the proposed fractional-order Ziegler-Nichols method is employed to tune parameters. Then the results will be compared with the result of using the traditional Ziegler-Nichols method tuned integer-order PID controller. One important step of the Ziegler-Nichols method is to find the ultimate gain *K _{u}* and the corresponding period

*T*for consistent oscillation. However, the trial-and-error approach would be time-consuming; it is desired to have the theoretical groundwork to make a reasonable estimation of

_{u}*K*to save trouble.

_{u}The Ruth-Hurwitz stability criterion [11] can be used to find the ultimate gain *K _{u}* in a linear system. Since the rotary pendulum system is nonlinear, it helps if the local linearization method is used to linearize the system around an operating point before using the Ruth-Hurwitz stability criterion to find the ultimate gain

*K*. The approach will be detailed in this section.

_{u}*10.4.1 LOCAL LINEARIZATION OF THE TWO-LINK ROTARY PENDULUM SYSTEM*

The equilibrium point of the rotary pendulum, i.e., *θ*_{2} = 0, is chosen to be the operating point for local linearization. First assume that *θ*_{2}≈ 0. Next linearize the dynamic Eqs. (10.31)–(10.32) by using the approximation sin *θ*_{2}≈ *θ*_{2}, cos *θ*_{2} ≈ 1, *θ*_{2} ≈ 0, and at equilibrium, then setting the higher order terms to zeros. Hence, the linearized model of the rotary pendulum can be expressed as:

(10.33)

(10.34)

(10.35)

where *Ĵ*_{11}=*J _{M}+J*

_{1+}

*m*

_{2}

*l*

_{1}

^{2},

*Ĵ*

_{12=}1/2

*m*

_{2}

*l*

_{1}

*l*

_{2}

*Ĵ*

_{22=}1/3

*m*

_{2}

*l*

_{2}

^{2}

Taking the Laplace transform of Eqs. (10.33)–(10.35) yields the following system of equations:

(10.36)

(10.37)

(10.38)

Since Eq. (10.38) contains only two unknowns Θ_{1} and Θ_{2}, Θ_{2} can be expressed in terms of Θ_{1} as:

(10.39)

By substituting the result into Eq. (10.37), Θ_{1} can be expressed in terms of *I* as:

(10.40)

and the transfer function of link 1’s angular displacement to armature current is:

(10.41)

Alternatively, Eq. (10.40) can be rewritten in a different way as:

(10.42)

Then by substituting Eq. (10.42) into Eq. (10.36), Θ_{1} can be expressed in terms of *E _{a}* as:

(10.43)

And the transfer function of link 1’s angular displacement to armature voltage can be expressed as:

(10.44)

where

*10.4.2 TRANSFER FUNCTION OF THE CLOSED-LOOP CONTROLLED ROTARY PENDULUM SYSTEM*

Figure 10.2 shows the block diagram of a typical closed-loop control system. The plant is the DC motor-driven rotary pendulum system, the reference input *r* is the desired link 1’s angular displacement input *θ _{1d}*, the system output

*y*is the actual link 1’s angular displacement output

*θ*, and the control effort

_{1a}*u*is the armature voltage

*e*in this research. When using a P controller with gain

_{a}*K*to control the rotary pendulum system, the closed-loop transfer function of the actual angular displacement output

_{p}*θ*to the desired angular displacement input

_{1a}*O*becomes:

_{1d}(10.45)

where *k _{1}= KT J_{22}* and

*k*. And its characteristic equation is:

_{2}= KT m_{2}l_{2}g(10.46)

*10.4.3 ESTIMATION OF THE ULTIMATE GAIN USING THE RUTH TABLE*

The following DC motor parameters and two-link rotary pendulum parameters are used: *L _{a}* = 4.6909 × 10

^{–3},

*R*= 2.5604 Ω,

_{a}*K*= 1.8259 × 10

_{T}^{–2}N·m/A,

*K*= 1.8259 ×10

_{F}^{–2}V/(rad/s),

*J*= 7.2019 × 10

_{M}^{–5}kg·m

^{2},

*B*= 9.1358 × 10

_{M}^{–5}N·m(rad/s),

*m*

_{1}0.056 kg,

*l*

_{1}= 0.16 m,

*J*

_{1}0.001569 kg·m

^{2},

*m*

_{2}= 0.022kg, and

*l*

_{2}= 0.16m. The coefficients of the characteristic equation calculated are

*a*

_{0}≈ 1.5691 × 10

^{–9},

*a*

_{1}≈ 8.5655 × 10

^{–7},

*a*

_{2}≈ 4.6355 × 10

^{–7},

*a*

_{3}≈ 1.949 × 10

^{–4},

*a*

_{4}≈ 1.959 × 10

^{–5},

*k*

_{1}≈ 3.4278 × 10

^{–6}, and

*k*

_{2}≈ 6.3051 × 10

^{–4}. And the Ruth table [11] is shown in Table 10.1.

S^{5} |
1.5691 × 10^{–9} |
4.6355 × 10^{–7} |
1.959 × 10^{–5} |

s^{4} |
8.5655 × 10^{–7} |
(1.949 × 10^{–4} + 3.4278 × 10^{–6}K )_{p} |
6.3051 × 10^{–4}K_{P} |

s^{3} |
A_{1} |
A_{2} |
– |

s^{2} |
B_{1} |
6.3051 × 10^{–4}K_{P} |
– |

s^{1} |
C_{1} |
– | |

s^{0} |
6.3051 × 10^{–4}K_{p} |
– | – |