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, _aD^α_t is a fractional-order integral. Note that when the infinitesimal increment 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 ₀D_t ^α–j–1 f (0) = 0, j= 0, 1, 2,…, n–1, then ℒ{₀D^α_tf(t)}= s^aF(s)

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

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_D 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:

(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₁ = K_p/T₁ , and K_D= K_PT_D. Furthermore, it’s noteworthy that the 3 terms inside the parenthesis of Eq. (10.9) have the same dimension of e(t), where the unit of T₁ and T_D (both of physical quantity “time”) is “second” used to cancel the unit of “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_u. Then, according to the Ziegler-Nichols method, the empirical PID parameters are assigned to be K_p= 0.6K_u , T_I= T_u/2 , and T_D=T_u/ 8, which is equivalent to have the alternative parameter set: K_p= 0.6K_u , K_I= 1.2K_u/T_u , and KD = 3K_uT_u/40.

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_D,μ are the proportional gain, the fractional-order integral gain, and the fractional-order derivative gain, respectively.

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_D cancels the unit of “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_I) and (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.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 τ_M= K_Ti, where K_T is the torque constant of the DC motor. Assume that the angular displacement of the motor shaft is θ and the load is τ_L, then the mechanical dynamics of the DC motor is:

(10.20)

where J_M and B_M are the rotational inertia and the rotational damping coefficient of the motor, respectively.

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_E is the back e.m.f. constant of the DC motor and K_E = K_T in the SI unit system. Then applying the Kirchhoff’s voltage law to the armature gives the electrical dynamics of the DC motor below:

(10.21)

where L_a and R_a are the inductance and resistance of the armature, respectively.

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₁ and m₂, respectively. The lengths of link 1 and 2 are l₁ and l₂, respectively. And the moment of inertia of link 1 around O is J₁. Assume that θ₁ and θ₂ are the angular displacements of link 1 and 2, respectively. Then the angular velocities of link 1 and 2 are θ̇̇₁and θ̇̇₂ respectively. Since the center O is on the motor shaft, it is straightforward θ₁ = θ.

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_nc,i is the 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 θ₁ and θ̱̱₂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 θ₁ = θ, then the whole system therewith can be represented by the dynamic model below.

(10.30)

(10.31)

(10.32)

10.4.1 LOCAL LINEARIZATION OF THE TWO-LINK ROTARY PENDULUM SYSTEM

The equilibrium point of the rotary pendulum, i.e., θ₂ = 0, is chosen to be the operating point for local linearization. First assume that θ₂≈ 0. Next linearize the dynamic Eqs. (10.31)–(10.32) by using the approximation sin θ₂≈ θ₂, cos θ₂ ≈ 1, θ₂ ≈ 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 Ĵ₁₁=J_M+J₁₊m₂l₁², Ĵ₁₂₌1/2m₂l₁l₂ Ĵ₂₂₌1/3m₂l₂²

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 Θ₁ and Θ₂, Θ₂ can be expressed in terms of Θ₁ as:

(10.39)

By substituting the result into Eq. (10.37), Θ₁ 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), Θ₁ 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 θ_1a, and the control effort u is the armature voltage e_a in this research. When using a P controller with gain K_p to control the rotary pendulum system, the closed-loop transfer function of the actual angular displacement output θ_1a to the desired angular displacement input O_1d becomes:

(10.45)

where k₁= KT J₂₂ and k₂= KT m₂l₂g. And its characteristic equation is:

(10.46)

FIGURE 10.2 Block diagram of a closed-loop control system.

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_a = 2.5604 Ω, K_T = 1.8259 × 10^–2 N·m/A, K_F = 1.8259 ×10^–2 V/(rad/s), J_M = 7.2019 × 10^–5 kg·m², B_M = 9.1358 × 10^–5 N·m(rad/s), m₁ 0.056 kg, l₁ = 0.16 m, J₁ 0.001569 kg·m², m₂ = 0.022kg, and l₂ = 0.16m. The coefficients of the characteristic equation calculated are a₀ ≈ 1.5691 × 10^–9, a₁ ≈ 8.5655 × 10^–7, a₂ ≈ 4.6355 × 10^–7, a₃ ≈ 1.949 × 10^–4, a₄ ≈ 1.959 × 10^–5, k₁ ≈ 3.4278 × 10^–6, and k₂ ≈ 6.3051 × 10^–4. And the Ruth table [11] is shown in Table 10.1.