Motion Control of a Quadrotor Aircraft via Singular Perturbations

1. Introduction

The study of the design and applications of unmanned aerial vehicles (UAVs) has been an important topic over the last few years. The current and potential applications of these vehicles are extensive, such as surveillance, rescue, espionage and entertainment, to mention a few.

UAVs can be controlled automatically to track a flight plan. From the theoretical point of view, UAVs are quite challenging since in most cases they are underactuated systems, which means that these systems have more degrees-of-freedom than control inputs.

The quadrotor is an UAV with six degrees-of-freedom and only four control inputs (four rotors). The control of the quadrotor is difficult because the nonlinearity of its dynamic model [1]. It is not possible to control all of the states at the same time. However, one can select some states to be controlled. A possible combination of controlled outputs can be the translational position and the yaw angle. The roll and the pitch angles can be controlled indirectly by introducing stable zero dynamics [2], [3].

Many methods to achieve attitude stabilization and trajectory tracking control have been proposed in recent years. In [4], the authors proposed a partial state backstepping sliding mode controller using a simplified model of the quadrotor. In [5] and [6], the application of the backstepping methodology to the control of quadrotors was revisited, showing that motion control can be achieve within the Lyapunov function framework. In [2], [3] and [7], a visual feedback is added to achieve stabilization of the system at the desired configuration. An experimental evaluation of a control approach was discussed in [8]. The work introduced in [9] proposed a neural network controller for the robust trajectory tracking control of the quadrotor posture. It is assumed that the model is partially known and that the system is subject to wind disturbances. By using the Newton-Euler formalism, a backstepping-based controller was introduced in [10]. [11], the feedback linearization technique is used to derive a tracking controller. An analysis of the internal dynamics is included in that work. More recently, in

[12], a modified sliding mode term is designed with the aim of controlling the attitude of a quadrotor aircraft subject to a class of disturbances. A motion controller composed of six nonlinear PD loops is introduced in

[13]. The controller guarantees that the error trajectories are uniformly bounded. The design of a controller based on the block control technique combined with the super twisting control algorithm for trajectory tracking of a quadrotor aircraft is introduced in [14].

The recent textbooks [15] and [16] can be consulted for further references on the history, modelling, motion planning and control of quadrotors.

Control approaches, such as feedback linearization, backstepping, sliding modes, adaptive control, to mention a few, match the idea that the resulting closed-loop system should satisfy Lyapunov's theory. As an alternative, the theory of singularly perturbed systems has also been recognized as a powerful tool in the analysis and design of controllers for different kinds of systems such as electromechanical, biological and electrical, see for example [17] and [22]. Essentially, this technique is based on analysing the convergence of the solution of a set of differential equations in two (or more) time scales.

Recent applications of the theory of singularly perturbed systems in the derivation of controllers for mechatronic systems can be found in [18], [19], [20] and [21]. The derivation of the controllers only considered the existence of gains assuring the time scale separation while numerical simulations and real-time experiments confirmed their practical viability.

The theory of singularly perturbed systems is well-adapted to the design of trajectory tracking controllers for quadrotor systems, since control laws can be easily synthesized by using different time scales. The error dynamics obtained by using a controller inspired from time scale separation should satisfy a set of conditions to guarantee that the error trajectories really converge to zero as time increases. The theory of singularly perturbed systems established in [17] provides a rigorous framework to prove the convergence of systems with different time scales.

In the next section, references about control designs for quadrotors inspired in time scale separation are provided.

A motion controller was proposed by Zuo [23], which is based on using a new command-filtered backstepping technique to stabilize the attitude and a linear tracking differentiator to eliminate the classical inner/outer-loop structure. The work by Bertrand et al. [24] deals with the regulation problem and the time scale separation by introducing an angular velocity controller in terms of a perturbing parameter together with the angular velocity commands. Esteban et al. [25] propose a control scheme for a radio-controlled helicopter. Stability proof is based on using a three time scale separation.

The design of control algorithms for quadrotor aircraft which ensures tracking of a desired time-varying trajectory via the theory of singularly perturbed systems still requires study since literature related to these subjects is limited.

In this document, a new control algorithm is introduced. The proposed controller has been derived by using the following time scale separation ideas:

An advantage of the proposed scheme is the addition of integral action, which is helpful in compensating model uncertainties and external disturbances. The overall closed-loop system stability is studied rigorously by using the Lyapunov theory and stability of singularly perturbed systems [17].

The present document is organized as follows. Section 2 is devoted to the quadrotor dynamics and the control goal. The proposed motion control algorithm is described in Section 3. Section 4 provides the proof that the state space origin of the resulting closed-loop system is exponentially stable. Section 5 describes the simulation results, which have the purpose of confirming that by using the proposed scheme the quadrotor matches exponentially a desired translation chaotic trajectory and a desired periodic yaw angle.

Simulations considering that the system is affected by perturbations and nonlinearities that appear in a practical implementation are provided in Section 6. There, the time scale separation introduced by the proposed controller is also confirmed. Additionally, Section 6 presents the results by using a desired translational helicoidal trajectory.

Finally, some concluding remarks and discussion on further research are provided in Section 7.

Notation: The variables and constants are represented by lowercase letters, Roman letters and Greek characters, x, θ ∈ ℝ. Vectors of dimension n are defined by bold lowercase letters and bold italic Greek letters, i.e., x , eta; ∈ ℝ ⁿ . The matrices of size n × m are defined by italic capital letters, i.e., B ∈ ℝ^n×m. I_3×3 denotes the identity matrix of size 3. An over dot placed on a variable, vector or matrix, indicates differentiation with respect to time, i.e., ẋ ∈ ℝ ⁿ . For x ∈ ℝ, we have that c_x = cos(x), s_x = sin(x) and t_x = tan(x). The notation arctan (x) ∈ ℝ means the standard inverse trigonometric function for all x ∈ ℝ. λ_min {A( x )} and λ_Max {A( x )} denote the minimum and maximum eigenvalues of a symmetric positive definite matrix A( x ), for all x ∈ ℝ ⁿ , respectively. ‖ x ⇊ ‖ = x 1 2 + ⋯ + x n 2 stands for the norm of vector x ∈ ℝ ⁿ . For vector x or matrix B with an arbitrary dimension, its transpose is expressed by x ^T or B^T, respectively. The cross product of two vectors a and b ∈ ℝ³ is represented by a × b ∈ ℝ³. A specific element of a vector or matrix is invoked by index notation, for example, x ₂, B₃₄ ∈ ℝ.

2. Quadrotor dynamics and control goal

2.1. Quadrotor dynamics

The equations of motion of the quadrotor are studied in the recent textbooks [27], [16] and [15]. The quadrotor dynamics is represented by a set of six ordinary differential equations of second order. In order to represent the quadrotor dynamics, we have used three differential equations derived in the inertial frame and three differential equations derived in the body frame. The inertial frame is denoted by E and the body frame as B. In Figure 1, a detailed representation of the quadrotor aircraft is given.

However, as discussed in [23], [28], [30], [31], [32] and [33], the quadrotor dynamics can also be represented by using three equations for the translation dynamics obtained in the inertial frame E and three equations for the orientation dynamics obtained in the body frame B. Neglecting the aerodynamic effects and the gyroscopic torque, that motion representation in state variables is given by

d d t p = p ̇ ,

(1)

d d t p ̇ = − g z e + T m R ( η ) z e ,

(2)

d d t η = W ( η ) Ω ,

(3)

I f d d t Ω = − Ω × I f Ω + τ a ,

(4)

where the constant m is the mass of the quadrotor, g is gravity acceleration, p = [x y z] ^T ∈ ℝ³ means the linear positions with respect to inertial reference frame, ⇊ η ⇊ = [ ϕ ⇈ θ ⇈ ψ ] T is the vector of Euler angles, which represent the roll, pitch and yaw angles in the inertial reference frame, I_f = diag{I_xx, I_yy, I_zz} is a matrix that contains the inertia constants, Ω = [ Ω x Ω y Ω z ] T is the vector of angular velocity with respect to the body reference frame,

z e = [ 0 0 1 ]

is a unitary vector in the z direction of the inertial frame which is used for compactness of the notation,

R ( η ) = [ c θ c ψ s θ c ψ s ϕ − s ψ c ϕ s θ c ψ c ϕ + s ψ s ϕ c θ s ψ s θ s ψ s ϕ + c ψ c ϕ s θ s ψ c ϕ − c ψ s ϕ − s θ c θ s ϕ c θ c ϕ ]

is an orthogonal rotation matrix that describes a clockwise/left-handed rotation with x-y-z convention, the lift forces f_i generated by the four propellers can be arranged in the vector

OPEN IN VIEWER

Figure 1.

The quadrotor aircraft and reference frames

f = [ f 1 f 2 f 3 f 4 ] ,

which will be considered the control input,

T = f 1 + f 2 + f 3 + f 4 = [ 1 1 1 1 ] f ,

(5)

is the total thrust and the vector

τ a = [ 0 l 0 − l − l 0 l 0 − λ λ − λ λ ] f ,

represents the torques applied to the quadrotor, which are generated by the lift forces of the rotor propellers.

Additionally, concerning the equation (3),

W ( η ) = [ 1 s ϕ t θ c ϕ t θ 0 c ϕ − s ϕ 0 s ϕ / c θ c ϕ / c θ ] .

In order to guarantee the non singularity of W( η ) for θ = ±π/2, we assume that the quadrotor does not perform acrobatic manoeuvres.

Notice that the relationship (3) from angular velocity in the body frame Ω ∈ ℝ³ to the time derivative of the Euler angles in the inertial frame η̇ ∈ ℝ³ has been obtained from the equation [29]

d d t R ( η ) = R ( η ) S ( Ω ) ,

where

S ( Ω ) = [ 0 − Ω z Ω y Ω z 0 − Ω x − Ω y Ω x 0 ] .

Throughout this paper the vector η ∈ ℝ³ will be denoted either as the orientation or as the angular position of the quadrotor.

2.2. Control problem formulation

Considering the vectors p _d (t) and η _d (t), which are the desired Cartesian position and the desired orientation, respectively, the vector

η ˜ = η d ( t ) − η ( t )

(6)

is the orientation error and

p ˜ = p d ( t ) − p ( t )

(7)

is the Cartesian position error.

Roughly speaking, the control problem consists in designing a control algorithm f (t) ∈ ℝ⁴ such that the limit

lim t → ∞ [ η ˜ ( t ) p ˜ ( t ) ] = 0,

(8)

is guaranteed for a compact set of initial conditions of the quadrotor state variables.

As will be seen later, the specification of the desired position p_d(t) should be differentiable at least up to order four and the desired angular position η_d(t) should be differentiable at least up to order two, while the desired roll angle φ_d(t) and desired pitch angle θ_d(t) should introduce stable position error dynamics, with the angle specification of the desired yaw ψ_d(t) assumed to be at least twice differentiable.

3. Development of the proposed algorithm

The proposed controller f ∈ ℝ⁴ has a complex structure. It will be described in a constructive form. First, consider that the control law f has the structure

f = B ( η ) μ

(9)

where µ ∈ ℝ⁴ is a signal to be defined and

B ( η ) = [ 1 4 c θ c ϕ 0 − 1 2 l − 1 4 λ 1 4 c θ c ϕ 1 2 l 0 1 4 λ 1 4 c θ c ϕ 0 1 2 l − 1 4 λ 1 4 c θ c ϕ − 1 2 l 0 1 4 λ ] .

(10)

Matrix B( η ) is obtained under the construction

B ( η ) = [ [ c θ c ϕ 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 ] [ 1 1 1 1 0 l 0 − l − l 0 l 0 − λ λ − λ λ ] ] − 1 .

which is singular for φ * , θ * ∈ ℝ such that cos(θ*) = 0 or cos(φ*) = 0.

The matrix B( η ) in (10) satisfies the following:

Property: Consider any vector µ ∈ ℝ⁴. Then,

[ 1 1 1 1 ] B ( η ) μ = 1 c ϕ c θ μ 1 ,

(11)

where µ₁ is the first element of the vector µ ∈ ℝ⁴.

By developing the calculations, the proof of property (11) is obtained.

By using property (11) it is possible to show that under the control law (9) the quadrotor dynamics (1)–(4) becomes

d d t p = p ̇ ,

(12)

d d t p ̇ = − g z e + 1 m [ R 13 R 23 R 33 ] μ 1 R 33 ,

(13)

d d t η = W ( η ) Ω ,

(14)

I f d d t Ω = − Ω × I f Ω + μ r ,

(15)

where µ₁ ∈ ℝ is the first component of µ ∈ ℝ⁴ and µ _r = [µ₂ µ₃ µ₄] ^T ∈ ℝ³ contains the other ones.

The block diagram implementation of the proposed controller can be appreciated in the Figure 2. As will be seen in the next section, the proposed scheme has an inner loop of angular velocity control and an outer loop of translational and angular position control.

3.1. Control of the angular dynamics

We solve the problem by using the following controller

μ r = I f Ω ̇ d + Ω × I f Ω d + K p w Ω ˜ + K i w ξ ,

(16)

d d t ξ = Ω ∼ = Ω d − Ω ,

(17)

where the K_pw and K_iw are 3×3 positive definite matrices. It should be noted that this control structure matches the PI velocity controller for the robot manipulators reported in [35], which is adapted from the structure of the well-known PD+ motion controller [34]. The addition of integral action has the advantage of improving the rejection of time-varying disturbances and model uncertainties.

To guarantee that the actual orientation η (t) converges to the desired orientation η _d (t) (to be explicitly computed later) the following kinematic-like controller is adopted [26]

Ω d = W ( η ) − 1 [ η ̇ d + K p o η ˜ ] ,

(18)

with K_po a 3 × 3 symmetric positive definite matrix.

Thus, the angular error dynamics is given by

d d t η ˜ = − K p o η ˜ + W ( η ) Ω ˜ ,

(19)

I f d d t Ω ˜ = − Ω × I f Ω ˜ − K p w Ω ˜ − K i w ξ ,

(20)

d d t ξ = Ω ∼ .

(21)

The only equilibrium point of system (19)–(21) is given at the state space origin, that is,

[ η ˜ Ω ˜ ξ ˜ ] = 0 ∈ ℝ 9 .

OPEN IN VIEWER

Figure 2.

Block diagram of the proposed motion control algorithm which is based on a primary model-based angular velocity tracking controller and secondary kinematic-like translational and angular position controller.

3.2. Control of the translational dynamics

By using the definition of p͂ (t) in (7), the translational error dynamics can be written as

d d t p ˜ = p ˜ ̇ ,

(22)

d d t p ˜ ̇ = p ̈ d + g z e − 1 m [ R 13 R 23 R 33 ] μ 1 R 33 .

(23)

Now, in order to design a way to stabilize exponentially the system (22)-(23), let us assume that

η ˜ ( t ) = 0 ⇒ η ( t ) = η d ( t ) , ∀ t ≥ t 0 .

(24)

Assumption (24) will be relaxed later when the stability analysis is presented. It is clear that in practice, either asymptotical or exponential tracking of the desired angular position η _d (t) can be guaranteed.

Invoking assumption (24), the translational error dynamics (22)-(23) becomes

d d t p ˜ = p ˜ ̇ ,

(25)

d d t p ˜ ̇ = p ̈ d + g z e − 1 m [ R 13 d R 23 d R 33 d ] μ 1 R 33 d .

(26)

where the d subscript in R_ijd means that the signal R_ij depends on η _d Equaling the right-hand side of (26) to some control vector γ p 1 ∈ ℝ 3 we have

p ̈ d + g z e − 1 m [ R 13 d R 23 d R 33 d ] μ 1 R 33 d = γ p 1 .

(27)

It is worth noticing that γ p 1 = [ γ p 11 γ p 12 γ p 13 ] T ∈ ℝ 3 .

If the desired Euler angles are specified so that equation (27) is satisfied, the translational error dynamics (25)-(26) is then given by

d d t p ˜ = p ˜ ̇ ,

(28)

d d t p ˜ ̇ = γ p 1 .

(29)

Therefore, the selection of the PD structure

γ p 1 = − K p p p ˜ − K d p p ˜ ̇

(30)

stabilizes exponentially the translational error dynamics (28)-(29).

Let us rewrite equation (27) as

[ R 13 d R 23 d R 33 d ] μ 1 m R 33 d = p ̈ d + g z e − γ p 1 = γ 1 ,

(31)

where γ p 1 = [ γ p 11 γ p 21 γ p 13 ] ∈ ℝ 3 . It is possible to look at the equation (31) as a nonlinear algebraic equation system where the unknown variables are φ d ( t ) , φ d ( t ) and μ 1 ( t ) ∈ ℝ , while ψ_d(t) is arbitrarily specified. A solution for the nonlinear algebraic equation system (31) is given by

η d = [ ϕ d ( t ) θ d ( t ) ψ d ( t ) ] = [ arctan ( c θ d [ s ψ d m γ 11 μ 1 − c ψ d m γ 12 μ 1 ] ) arctan ( c ψ d m γ 11 μ 1 + s ψ d m γ 12 μ 1 ) ψ d ( t ) ] ,

(32)

μ 1 = m γ 13 = m [ p ̈ d 3 + g − γ p 13 ] ,

(33)

where γ₁₁, γ₁₂ and γ₁₃ are obtained from (31).

Notice that the computation of η ̇ d and η ̈ d is necessary to implement the control law µ _r in (16), since they are required in the desired angular velocity Ω _d in (18) and in the desired angular acceleration Ω ̇ d which is explicitly given by

Ω ̇ d = W ( η ) − 1 [ η ̈ d + K p o η ˜ ̇ ] + [ d d t [ W ( η ) − 1 ] ] [ η ̇ d + K p o η ˜ ] .

(34)

At the same time, the calculation of η ̇ d and η ̈ d requires the signals γ ̇ p 1 and γ ̈ p 1 , which, under definition (30), depend on the unmeasurable signals p ∼ ̈ and p ∼ ( 3 ) .

In order to respect the assumption that the only measured signals are those of the state [ p T p ̇ T η T Ω T ] T ∈ ℝ 12 , the following definition of γ _p1 is proposed:

γ ̇ p 1 = − K f 1 [ γ p 1 − γ p 2 ]

(35)

γ ̇ p 2 = − K f 2 [ γ p 2 − [ − K p p p ˜ − K d p p ˜ ̇ ] ] ,

(36)

where K_f1 and K_f2 are 3 × 3 symmetric positive definite matrices. Thus, the larger the numerical value of K_f1 and K_f2, the faster the convergence of the signal γ _p1(t) to [ − K p p p ∼ ( t ) − K d p p ∼ ̇ ( t ) ] . Therefore, the closed-loop dynamics of the translational error (28)-(29) can be recovered in an exponential way.

The explicit definition of φ ̇ d , θ ̇ d , φ ̈ d and θ ̈ d are provided in Appendix since they are necessary to compute Ω _d in (18) and Ω ̇ d in (34).

The rest of this section is devoted to explaining the steps to be followed in order to implement the proposed motion controller.

4. Overall closed-loop system stability via three time scale

Consider the following parametrization of the control gains

K p w = K ̄ p w ɛ 1 , K i w = K ̄ i w ɛ 1 ,

(37)

K f 1 = K ̄ f 1 ɛ 1 , K f 2 = K ̄ f 2 ɛ 1 ,

(38)

and

K p o = K ̄ p o ɛ 2 ,

(39)

where ε₁ and ε₂ are positive constants denoted as perturbing parameters, which satisfy

ɛ 2 > > ɛ 1 > 0.

(40)

Therefore, it is possible to show that the overall closed-loop system can be written as

d d t p ˜ = p ˜ ̇ ,

(41)

d d t p ˜ ̇ = p ̈ d + g z e − 1 m [ R 13 ( η d , η ˜ ) R 23 ( η d , η ˜ ) R 33 ( η d , η ˜ ) ] μ 1 R 33 ( η d , η ˜ ) ,

(42)

d d t ξ = Ω ˜ ,

(43)

ɛ 2 d d t η ˜ = − K ̄ p o η ˜ + ɛ 2 W ( η d , η ˜ ) Ω ∼ ,

(44)

ɛ 1 d d t Ω ˜ = − I f − 1 [ ɛ 1 Ω ( ɛ 2 ) × I f Ω ˜ + K ̄ p w Ω ˜ + K ̄ i w ξ ] ,

(45)

ɛ 1 d d t γ p 1 = − K ̄ f 1 [ γ p 1 − γ p 2 ] ,

(46)

ɛ 1 d d t γ p 2 = − K ̄ f 2 [ γ p 2 − [ − K p p p ˜ − K d p p ˜ ̇ ] ] ,

(47)

where, by virtue of the angular velocity error Ω ∼ in (17), Ω ( ɛ 2 ) = Ω d ( ɛ 2 ) - Ω ∼ , with Ω_d defined in (18). Notice that we have written Ω (ε₂), since Ω _d (ε₂), in order to emphasize that the perturbing parameter ε₂ is present into the definition of Ω _d . It is worth noticing that the vector η d ∈ ℝ 3 is obtained from (32) while the scalar µ₁ is defined in (33).

The overall closed-loop system (41)–(47) has the standard form of a singularly perturbed system with a three time scale, which has been conveniently derived through the parametrization of the controller gains (37)–(39).

The stability proof of the system (41)–(47) and the obtention of the time scales are based on the following steps:

The first step in the stability proof, by assumption (40), is the separation of the system in a twofold time scale. Thus,

x = [ x s − x m ] = [ p ˜ p ˜ ̇ ξ - η ˜ ]

are the slow states with respect to

z = [ Ω ˜ γ p 1 γ p 2 ] ,

which at the same time represents the states with fast time scale.

For the second step, exponential convergence of the solutions x (t) is proven by considering that

x s = [ p ˜ p ˜ ̇ ξ ]

are slow states with respect to

x m = η ˜ .

Thus, we have that x_s represents the states with slow time scale, x _m are the states with medium time scale, while z denotes the states with fast time scale. Notice, however, that the slow time and medium time states in [x^T_s x ^T_m] ^T are both slow in comparison to the fast states in z. Our main theoretical result is established as follows:

Proposition 1: Consider the singularly perturbed system (41)–(47). Then, there always exists ε*₁ > 0 such that the state space origin of the closed-loop system (41)–(47) is exponentially stable with ε*₁ > ε₁ > 0.

Proof: The proof of Proposition 1 is achieved by verifying the conditions of the theorem 9.3 in [17]. The proof of Proposition 1 is set down in the following five items:

1) The closed-loop system (41)–(47) has a unique equilibrium point at

[ [ p ∼ T p ∼ ̇ T ξ T n ∼ T ] T [ Ω ∼ T γ p 1 T γ p 2 T ] T ] = 0.

2) One important step in the application of theorem 9.3 [17] is to verify that the subsystem (45)–(47) has isolated roots of the form

z = z * = h ( t , x ) ,

with x = [ p ∼ T p ∼ ̇ T ξ T η ∼ T ] T and z = [ Ω ∼ T γ p 1 T γ p 2 T ] T , for ɛ 1 = 0 . The above expression is known as quasi-steady-state solution [17]. Hence by substituting ε₁ = 0 into (45)–(47) we have that any isolated root z * ∈ ℝ⁹ satisfies

Ω ˜ * = − K ̄ p w − 1 K ̄ i w ξ ,

(48)

γ p 1 * = − K p p p ˜ − K d p p ˜ ̇ ,

(49)

γ p 2 * = − K p p p ˜ − K d p p ˜ ̇ .

(50)

It is noteworthy that the isolated root (49) evaluated in (31)–(33) gives the specific value η * _d for the desired orientation.

3) The system (41)–(47) and the isolated roots (48)–(50) have bounded partial derivatives in compact sets.

4) By considering ε₁ = 0 and using the isolated roots (48)–(50) in equations (41)–(44), the following slow dynamics system is obtained:

d d t p ˜ = p ˜ ̇ ,

(51)

d d t p ˜ ̇ = p ̈ d + g z e − 1 m [ R 13 ( η d * , η ˜ ) R 23 ( η d * , η ˜ ) R 33 ( η d * , η ˜ ) ] μ 1 R 33 ( η d * , η ˜ ) ,

(52)

d d t ξ = − K ̄ p w − 1 K ̄ i w ξ ,

(53)

ɛ 2 d d t η ˜ = − K ̄ p o η ˜ − ɛ 2 W ( η d * , η ˜ ) K ̄ p w − 1 K ̄ i w ξ .

(54)

As previously mentioned, the states

x = [ x s T | x m T ] T = [ p ˜ T p ˜ ̇ T ξ T | η ˜ T ] T

are slow with respect to

z = [ Ω ̈ T γ p 1 T γ p 2 T ] T .

The local exponential stability of the state space origin of the slow dynamics (41)–(44) can be derived by using again a time scale approach, interpreting x _s as slow with respect to x _m . For the sake of simplicity, we have stated such a result in Proposition 2, which will be given later.

5) To obtain the boundary layer system, the change of variable

y = z − h ( t , x ) ,

is defined. Let us remember that x = [ p ∼ T p ∼ ̇ T ξ T η T ] and z = [ Ω ∼ T γ p 1 T γ p 2 T ] , with z* = h (t,x) given in (48)–(50).

By deriving y with respect to σ = t/ε₁ and setting ε₁ = 0, we obtain the boundary layer system

d d σ y 1 = − I f − 1 K ̄ p w y 1 ,

(55)

d d σ y 2 = − K ̄ f 1 [ y 2 − y 3 ] ,

(56)

d d σ y 3 = − K ̄ f 2 y 3 ,

(57)

which is a linear system. The Lyapunov function for system (55)–(57) is defined as

W ( y 1 , y 2 , y 3 ) = 1 2 y 1 T I f y 1 + β 1 2 y 2 T y 2 + 1 2 y 3 T y 3 , β > 0 ,

whose derivative with respect to the fast time scale σ is given by

d d σ W ( y 1 , y 2 , y 3 ) = − [ y 1 y 2 y 3 ] T [ K ̄ p w 0 0 0 β K ̄ f 1 − 1 2 β K ̄ f 1 0 0 − 1 2 β K ̄ f 1 K ̄ f 2 ] [ y 1 y 2 y 3 ] .

By invoking Sylvester's theorem, the condition

β < 4 λ m i n { K ̄ f 1 } λ m i n { K ̄ f 2 } λ M a x { K ̄ f 1 } 2

guarantees the inequality

d d σ W ( y 1 , y 2 , y 3 ) ≤ − κ 1 W ( y 1 , y 2 , y 3 ) , κ 1 > 0.

Therefore, because (55)–(57) is a linear system, we have sufficient conditions to claim that there is an arbitrarily large compact set B_ρ such that

[ y 1 ( 0 ) T y 2 ( 0 ) T y 3 ( 0 ) T ] T ∈ B ρ

implies that

[ y 1 ( σ ) T y 2 ( σ ) T y 3 ( σ ) T ] T → 0

as the scaled time σ → ∞.

By theorem 9.3 [17], there are sufficient conditions to claim that there exists ε*₁ such that ε*₁ > ε₁ > 0 guarantees that there is an arbitrarily large compact set R_A of initial conditions where the trajectories of the overall closed-loop system (41)–(47) achieves the following limit

lim t → ∞ [ x ( t ) z ( t ) ] = lim t → ∞ [ p ˜ ( t ) p ˜ ̇ ( t ) ξ ( t ) η ˜ ( t ) Ω ˜ ( t ) γ p 1 ( t ) γ p 2 ( t ) ] = 0 .

The result stated in Proposition 1 means that the goal (8) holds by using the control algorithms in Section 3.3.

The proof of Proposition 1 is based on the assumption of the local exponential convergence of the slow dynamics (51)–(54). In order to relax such an assumption, the explicit proof of exponential convergence of the slow time states x s = [ p ∼ ( t ) T p ∼ ̇ ( t ) T ξ ( t ) T ] T and the medium time states x _m = ῆ(t) is provided as follows:

OPEN IN VIEWER

Table 1.

Description of the quadrotor parameters.

Notation	Description	Unit	Value
g	gravity acceleration	[m/s2]	9.80665
m	Quadrotor mass	[kg]	0.141
l	Distance between the quadrotor centre of mass and the rotation axis of propeller	[m]	0.175
Ixx, Iyy	Moment of inertia along φ and θ	[kg m2]	2.1 × 10−3
Izz	Moment of inertia along ψ	[kg m2]	4.2 × 10−3
λ	Propeller parameter related to the reactive torque	[m]	10.2 × 10−3

OPEN IN VIEWER

Figure 3.

Actual quadrotor path p (t) using ε₁ = 1.0, ε₂ = 0.5 and simulation time of 300.0 [sec]. The grey point represents the initial position.

Proposition 2: Consider the singularly perturbed system (51)–(54), which represents the slow dynamics in the proof of Proposition 1. Then, there always exists ε*₂ > 0 such that the state space origin of the closed-loop system (51)–(54) is exponentially stable with ε*₂ > ε₂ > 0.

Proof: The proof of Proposition 2 is achieved by verifying the conditions of the theorem 9.3 in [17]. The proof of Proposition 2 is set down in the following five items:

1) The closed-loop system (51)–(54) has a unique equilibrium point at [ [ p ∼ T p ∼ ̇ T ξ T η ∼ T ] T [ Ω ∼ T γ p 1 T γ p 2 T ] T ] = 0 .

2) The quasi-steady-state solution

z ̄ = z ̄ * = x m = h ̄ ( t , x s )

is given by the isolated root

η ˜ * = 0.

(58)

3) The system (51)–(54) and the isolated root (58) have bounded partial derivatives in compact sets.

4) The slow dynamics is obtained by substituting the isolated root (58) into the system (51)–(53), and invoking the equation (27), i.e.,

d d t p ˜ = p ˜ ̇ ,

(59)

d d t p ˜ ̇ = − K p p p ˜ − K d p p ˜ ̇ ,

(60)

d d t ξ = − K ̄ p w − 1 K ̄ i w ξ .

(61)

Let us consider the Lyapunov function

V ( p ∼ , p ̇ , ξ ) = 1 2 [ α p ∼ + p ∼ ̇ ] T [ α p ∼ + p ∼ ̇ ] + 1 2 p ˜ T [ K p p + α K d p − α 2 I 3 × 3 ] p ˜ + 1 2 ξ T K ̄ i w ξ ,

(62)

with α > 0 small enough so that

K d p − α I 3 × 3 > 0.

The time derivative of V ( p ∼ , p ̇ , ξ ) along the trajectories of the slow dynamics (59)–(61) is given as

V ̇ ( p ˜ , p ̇ , ξ ) = − p ˜ ̇ T [ K d p − α I 3 × 3 ] p ˜ ̇ − α p ˜ T K p p p ˜ − ξ T K ̄ i w K ̄ p w − 1 K ̄ i w ξ ,

(63)

The right-hand side of equation (63) is globally negative definite for a small enough α > 0. Finally, it is possible to prove that the inequality

V ̇ ( p ˜ , p ̇ , ξ ) ≤ − κ 2 V ( p ˜ , p ̇ , ξ ) ,

with some κ₂ > 0, is satisfied for all [ p ∼ T p ̇ T ξ T ] T ∈ B r . Therefore, the state-space origin of the slow dynamics (59)–(61) is exponentially stable. In other words,

lim t → ∞ [ p ˜ ( t ) T p ˜ ̇ ( t ) T ξ ( t ) T ] T = 0

with exponential convergence rate as time t increases.

5) By deriving

y ̄ = z ̄ − h ̄ ( t , x s ) = x m − h ̄ ( t , x s )

with respect to scaled time ᾱ = t/ε₂ and setting ε₂ = 0 in (54), the following boundary layer system is obtained:

d d α ̄ η ˜ = − K ̄ p o η ˜

(64)

which is exponentially stable for all K̄_po > 0.

By theorem 9.3 [17], there are sufficient conditions to claim that there exists ε*₂ such that ε*₂ > ε₂ > 0 guarantees that there is an arbitrarily large compact set R_B of initial conditions where

lim t → ∞ [ p ˜ ( t ) T p ˜ ̇ ( t ) T ξ ( t ) T η ˜ ( t ) ] T = 0 .

An interpretation of the results established in Proposition 1 and 2 is the existence of control gains (37)–(39) so as the local exponential stability of the closed-loop (41)–(47) is assured. Additionally, Proposition 1 and 2 guarantee the existence of perturbing parameters ε₁ and ε₂ such that

ɛ 2 * > ɛ 2 > > ɛ 1 * > ɛ 1 > 0,

which can be used as a tuning guideline.

OPEN IN VIEWER

Figure 4.

Slow time scale states: the position errors p͂ >(t), translational velocity errors p ∼ ̇ ( t ) , and integral of the angular velocity error ξ ( t ) for the different values of the perturbing parameters ε₁ and ε₂. Medium time scale states: orientation error η ∼ ( t ) for the different values of the perturbing parameters ε₁ and ε₂. Some of the signals are scaled in magnitude in order to improve their visualization with respect to the other ones.

OPEN IN VIEWER

Figure 5.

Fast time scale: Time evolution of y 1 ( t ) = Ω ∼ ( t ) + K ̄ p w − 1 K ̄ i w ξ ( t ) , y 2 ( t ) = γ p 1 ( t ) + K p p p ∼ ( t ) + K d p p ∼ ̇ ( t ) and y 3 ( t ) = γ p 2 ( t ) + K p p p ∼ ( t ) + K d p p ∼ ̇ ( t ) for the different values of the perturbing parameters ε₁ and ε₂. Some of the signals are scaled in magnitude in order to improve their visualization with respect to the other ones.

5. Simulation results

We have used the Matlab and Simulink toolboxes to implement the system described in Figure 1. In the simulations we considered that the initial conditions of the quadrotor and the controller are null. The parameters of the quadrotor are given in Table 1.

The desired position is defined as

p d ( t ) = [ x d ( t ) y d ( t ) z d ( t ) ] = a ^ o [ x ^ ( t ) y ^ ( t ) z ^ ( t ) + 1 ] ∈ ℝ 3 ,

(65)

where x̂(t), ŷ(t), ẑ(t), are obtained by means of the Rabinovich-Fabrikant equations [36], that is,

[ x ^ ̇ y ^ ̇ z ^ ̇ ] = τ ^ o [ y ^ [ z ^ − 1 + x ^ 2 ] + γ x ^ x ^ [ 3 z ^ + 1 − x ^ 2 ] + γ y ^ − 2 z ^ [ α + x ^ y ^ ] ] ∈ ℝ 3 ,

(66)

being γ = 0.1 y α = 0.3. Besides, τ ^ o = 0.2 , a ^ o = 1 . Initial conditions are x̂(0) = −0.16, ŷ(0) = 0.1 and ẑ(0) = 0.1. The signals p ̇ d , p ̈ d , p d ( 3 ) and p (⁽⁴⁾ _d are obtained by a successive time derivative of p _d (t) in (66) and substituting in these the right-hand side of equation (66).

The desired yaw angle ψ_d(t) is defined as the time-periodic function

ψ d = a 1 sin ( 2 π τ 1 − 1 t ) + a 2 sin ( 2 π τ 2 − 1 t ) ∈ ℝ ,

(67)

with τ₁ = 30 [sec], τ₂ = 4 [sec], a₁ = 0.5 [rad] and a₂ = 0.1 [rad].

Control gains have been selected using the parametrization (37)-(38) and (39). In particular,

K p w = ɛ 1 − 1 d i a g { 0.005,0.005,0.075 } , K i w = ɛ 1 − 1 d i a g { 0.005,0.005,0.075 } , K f 1 = ɛ 1 − 1 d i a g { 25,25,25 } , K f 2 = ɛ 1 − 1 d i a g { 50,50,50 } , and K p o = ɛ 2 − 1 d i a g { 5,5,5 } , where

(68)

ɛ 1 = 0.1, 0.4, 0.7, 1.0,

(69)

and

ɛ 2 = 0.5, 1.5, 2.5, 3.5, 4.5.

(70)

On the other hand, the gains of the translational controller are given by

K p p = d i a g { 2,0.75,10 } , K d p = d i a g { 3,1,20 } .

(71)

The purpose of the numerical simulations is to confirm that the actual quadrotor trajectories p (t) and ψ(t) converge to the desired position trajectory p _d (t) in (65) and to the desired yaw angle ψ_d(t) in (67). In addition, with the presented simulations the time scale separation introduced by the controller is verified. The simulations show the time evolution of the slow time scale states

x ̄ ( t ) = [ x ̄ 1 ( t ) x ̄ 2 ( t ) x ̄ 3 ( t ) ] = [ p ˜ ( t ) p ˜ ̇ ( t ) ξ ( t ) ] ,

OPEN IN VIEWER

Figure 6.

Time evolution of the norm of the slow time scale states ‖x̄(t)‖, the norm of the medium time scale states ‖z̄(t)‖ and the norm of the fast time scale states ‖y(t)‖ for the different values of the perturbing parameters ε₁ and ε₂.

the medium time scale states

z ̄ ( t ) = η ˜ ( t )

and the fast time scale states

y ( t ) = z ( t ) − h ( t , x ( t ) ) ,

where the quasi-steady-state solution h (t, x (t)) is given explicitly in (48)–(50), being

x ( t ) = [ x ̄ ( t ) T | z ̄ ( t ) T ] T : = [ p ˜ T p ˜ ̇ T ξ T | η T ] T .

More specifically, the fast states are

y = [ y 1 ( t ) y 2 ( t ) y 3 ( t ) ] = [ Ω ˜ ( t ) + K ̄ p w − 1 K ̄ i w ξ ( t ) γ p 1 ( t ) + K p p p ˜ ( t ) + K d p p ˜ ̇ ( t ) γ p 2 ( t ) + K p p p ˜ ( t ) + K d p p ˜ ̇ ( t ) ] .

Figure 3 shows the actual quadrotor path p (t) for a time interval of 300.0 [sec] and the control gains previously defined with ε₁ = 1 and ε₁ = 0.5. The desired trajectory corresponds to the attractor system given by (66).

The rest of the Figures show the response of the system for the different values of the perturbation parameters ε₁ and ε₂. In particular, a simulation has been performed for each combination of the perturbing parameters (ε₁, ε₂), which gives a total of 20 simulations.

These simulations have been carried out in two nested “for” loops. Each iteration corresponds to a simulation with a combination of parameters (ε₁, ε₂). In the outer “for” loop, the parameter ε₁ is increased and in the inner “for” loop, the parameter ε₂ is increased as (69) and (70), respectively. The results are shown in different coloured lines. For each colour, lighter lines mean the former simulations and darker lines mean the latter simulations.

Figure 4 depicts the slow time scale states, given by position errors p͂ (t), velocity errors p ∼ ̇ ( t ) and the integral of the angular velocity error ξ (t). Due to the negligible differences between the compared signals, the plots of p͂ (t) and p ∼ ̇ ( t ) are overlapped for the different values of ε₁ and ε₂. In addition, the convergence of the signal ξ (t) for the different values of ε₁ is barely affected. Figure 4 also shows the medium time scale states, which are given by the orientation errors φ(t), θ͂(t) and ψ͂. The main observation is that the convergence rate is increased if ε₂ is decreased, which is predicted by theory in Proposition 1 and 2.

OPEN IN VIEWER

Figure 7.

Total applied force T(t) and applied torque τ (t) for the different values of ε₁ and ε₂.

OPEN IN VIEWER

Figure 8.

Block diagram of the simulations with practical disturbances.

In Figure 5, the time evolution of the states with fast time scale y (t) can be seen. Notice that the smaller the value of the parameter ε₁, the faster the convergence of the signals y ₁(t), y ₂(t) and y ₃(t). Another important observation is that the parameter ε₂ has no effect on the convergence of y (t), which is consistent with Proposition 2.

Figure 6 shows the time evolution of the norm of the slow time scale states ‖ x̄ (t)‖, the norm of the medium time scale states ‖ z̄ (t)‖ and the norm of the fast time scale states ‖ ȳ (t)‖.

The results in Figures 4, 5 and 6 confirm that the perturbation parameter ε₂ has more effect on the convergence of the medium time scale states, while the perturbing parameter ε₁ causes the signal y (t) to converge to zero much faster than x (t) (the slow and medium time scale states).

Finally, the total applied force T(t) and the applied torque τ_a(t) for the different values of the parameters ε₁ and ε₂ are depicted in Figure 7. It is appreciated that both signals are well-behaved for all time.

6. Simulations with practical disturbances

Simulations have been carried out incorporating the effect of practical disturbances. In Figure 8, a block diagram showing the incorporation of disturbances in the simulations is provided. The blocks are described as follows:

OPEN IN VIEWER

Figure 9.

Left hand-side plots correspond to the results for the chaotic-type trajectory in (65)-(66). Right hand-side plots show the results for the helicoidal trajectory in (72).

Under the conditions described above, two different translational position trajectories have been numerically simulated: