Motion Control of a Quadrotor Aircraft via Singular Perturbations

Abstract

In this paper, a new motion controller for a quadrotor aircraft is introduced. A reformulation of the control inputs of the dynamic model is discussed and then the control algorithm is given in a constructive form. The stability proof of the state space origin of the overall closed-loop system relies on the theory of singularly perturbed systems. Numerical simulations corroborate the viability of the proposed control scheme and the conclusions concerning stability. A set of simulations under practical conditions is also presented, where the system is affected by different types of disturbances and nonlinearities such as noise and actuator saturation.

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:

• Differential equations with fast time scale are generated by using an angular proportional-integral velocity controller.

• Differential equations with medium time scale are generated by a kinematic-like scheme, which generates angular position commands.

• Differential equations with slow time scale are generated by the translational position controller. The translational dynamics is controlled indirectly by designing proper orientation commands. As seen later, those commands are designed under the assumption that the angular position is null.

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. $\Vert x\text{⇊}\Vert =\sqrt{x_1^2+\cdots +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

$\frac{d}{dt}p=\dot{p},$

(1)

$\frac{d}{dt}\dot{p}=-gz{\hspace{0.17em}}_e+\frac{T}{m}R(\mathit{\eta })z_e\mathrm{,}$

(2)

$\frac{d}{dt}\mathit{\eta }=W(\mathit{\eta })\mathbf{\Omega }\mathrm{,}$

(3)

$I_f\frac{d}{dt}\mathbf{\Omega }=-\mathbf{\Omega }\times I_f\mathbf{\Omega }+{\mathit{\tau }}_a\mathrm{,}$

(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, $\text{⇊}\mathit{\eta }\text{⇊}=[\varphi \text{⇈}\theta \text{⇈}\psi {]}^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, $\Omega ={\left[{\Omega }_x{\Omega }_y{\Omega }_z\right]}^T$ is the vector of angular velocity with respect to the body reference frame,

$z_e=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]$

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

$R(\mathit{\eta })=\left[\begin{array}{ccc}{c}_{\theta }{c}_{\psi }& s_{\theta }{c}_{\psi }{s}_{\varphi }-s_{\psi }{c}_{\varphi }& s_{\theta }{c}_{\psi }{c}_{\varphi }+s_{\psi }{s}_{\varphi }\\ c_{\theta }{s}_{\psi }& s_{\theta }{s}_{\psi }{s}_{\varphi }+c_{\psi }{c}_{\varphi }& s_{\theta }{s}_{\psi }{c}_{\varphi }-c_{\psi }{s}_{\varphi }\\ -s_{\theta }& c_{\theta }{s}_{\varphi }& c_{\theta }{c}_{\varphi }\end{array}\right]$

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

Figure 1.

The quadrotor aircraft and reference frames

$f=\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\\ f_4\end{array}\right]\mathrm{,}$

which will be considered the control input,

$T=f_1+f_2+f_3+f_4=\left[\begin{array}{cccc}1& 1& 1& 1\end{array}\right]f\mathrm{,}$

(5)

is the total thrust and the vector

$\mathit{\tau }{\hspace{0.17em}}_a=\left[\begin{array}{cccc}0& l& 0& -l\\ -l& 0& l& 0\\ -\lambda & \lambda & -\lambda & \lambda \end{array}\right]f\mathrm{,}$

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

Additionally, concerning the equation (3),

$W(\mathit{\eta })=\left[\begin{array}{ccc}1& s_{\varphi }{t}_{\theta }& c_{\varphi }{t}_{\theta }\\ 0& c_{\varphi }& -s_{\varphi }\\ 0& s_{\varphi }/c_{\theta }& c_{\varphi }/c_{\theta }\end{array}\right]\mathrm{.}$

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]

$\frac{d}{dt}R(\mathit{\eta })=R(\mathit{\eta })S(\mathbf{\Omega })\mathrm{,}$

where

$S(\mathbf{\Omega })=\left[\begin{array}{ccc}0& -{\mathrm{\Omega }}_z& {\mathrm{\Omega }}_y\\ {\mathrm{\Omega }}_z& 0& -{\mathrm{\Omega }}_x\\ -{\mathrm{\Omega }}_y& {\mathrm{\Omega }}_x& 0\end{array}\right]\mathrm{.}$

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

$\tilde{\eta }=\mathit{\eta }{\hspace{0.17em}}_d(t)-\mathit{\eta }(t)$

(6)

is the orientation error and

$\tilde{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

$\underset{t\to \infty }{{\displaystyle \lim }}\left[\begin{array}{c}\tilde{\eta }(t)\\ \tilde{p}(t)\end{array}\right]=\mathrm{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(\mathit{\eta })\mathit{\mu }$

(9)

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

$B(\mathit{\eta })=\left[\begin{array}{cccc}\frac{1}{4c_{\theta }{c}_{\varphi }}& 0& -\frac{1}{2l}& -\frac{1}{4\lambda }\\ \frac{1}{4c_{\theta }{c}_{\varphi }}& \frac{1}{2l}& 0& \frac{1}{4\lambda }\\ \frac{1}{4c_{\theta }{c}_{\varphi }}& 0& \frac{1}{2l}& -\frac{1}{4\lambda }\\ \frac{1}{4c_{\theta }{c}_{\varphi }}& -\frac{1}{2l}& 0& \frac{1}{4\lambda }\end{array}\right]\mathrm{.}$

(10)

Matrix B( η ) is obtained under the construction

$B(\mathit{\eta })={\left[\left[\begin{array}{cccc}{c}_{\theta }{c}_{\varphi }& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 1& 1& 1\\ 0& l& 0& -l\\ -l& 0& l& 0\\ -\lambda & \lambda & -\lambda & \lambda \end{array}\right]\right]}^{-1}\mathrm{.}$

which is singular for ${\phi }^{*}\mathrm{,}{\theta }^{*}\in ℝ$ such that cos(θ*) = 0 or cos(φ*) = 0.

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

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

$\left[\begin{array}{cccc}1& 1& 1& 1\end{array}\right]B(\mathit{\eta })\mathit{\mu }=\frac{1}{c_{\varphi }{c}_{\theta }}{\mu }_1\mathrm{,}$

(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

$\frac{d}{dt}p=\dot{p}\mathrm{,}$

(12)

$\frac{d}{dt}\dot{p}=-gz{\hspace{0.17em}}_e+\frac{1}{m}\left[\begin{array}{c}{R}_{13}\\ R_{23}\\ R_{33}\end{array}\right]\frac{{\mu }_1}{R_{33}}\mathrm{,}$

(13)

$\frac{d}{dt}\mathit{\eta }=W(\mathit{\eta })\mathbf{\Omega }\mathrm{,}$

(14)

$I_f\frac{d}{dt}\mathbf{\Omega }=-\mathbf{\Omega }\times I_f\mathbf{\Omega }+\mathit{\mu }{\hspace{0.17em}}_r\mathrm{,}$

(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

${\mathit{\mu }}_r=I_f\dot{\Omega }{\hspace{0.17em}}_d+\mathbf{\Omega }\times I_f\mathbf{\Omega }{\hspace{0.17em}}_d+K_{pw}\tilde{\Omega }+K_{iw}\mathit{\xi }\mathrm{,}$

(16)

$\frac{d}{dt}\mathit{\xi }=\stackrel{\sim }{\mathbf{\Omega }}=\mathbf{\Omega }{\hspace{0.17em}}_d-\mathbf{\Omega }\mathrm{,}$

(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]

$\mathbf{\Omega }{\hspace{0.17em}}_d=W{(\mathit{\eta })}^{-1}[\dot{\eta }{\hspace{0.17em}}_d+K_{po}\tilde{\eta }]\mathrm{,}$

(18)

with K_po a 3 × 3 symmetric positive definite matrix.

Thus, the angular error dynamics is given by

$\frac{d}{dt}\tilde{\eta }=-K_{po}\tilde{\eta }+W(\mathit{\eta })\tilde{\Omega }\mathrm{,}$

(19)

$I_f\frac{d}{dt}\tilde{\Omega }=-\mathbf{\Omega }\times I_f\tilde{\Omega }-K_{pw}\tilde{\Omega }-K_{iw}\mathit{\xi }\mathrm{,}$

(20)

$\frac{d}{dt}\mathit{\xi }=\stackrel{\sim }{\mathbf{\Omega }}\mathrm{.}$

(21)

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

$\left[\begin{array}{c}\tilde{\eta }\\ \tilde{\Omega }\\ \tilde{\xi }\end{array}\right]=0\in {ℝ}^9\mathrm{.}$

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

$\frac{d}{dt}\tilde{p}=\dot{\tilde{p}}\mathrm{,}$

(22)

$\frac{d}{dt}\dot{\tilde{p}}=\ddot{p}{\hspace{0.17em}}_d+gz{\hspace{0.17em}}_e-\frac{1}{m}\left[\begin{array}{c}{R}_{13}\\ R_{23}\\ R_{33}\end{array}\right]\frac{{\mu }_1}{R_{33}}\mathrm{.}$

(23)

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

$\tilde{\eta }(t)=0\Rightarrow \mathit{\eta }(t)=\mathit{\eta }{\hspace{0.17em}}_d(t)\mathrm{,}\hspace{0.17em}\forall \hspace{0.17em}t\ge t_0\mathrm{.}$

(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

$\frac{d}{dt}\tilde{p}=\dot{\tilde{p}}\mathrm{,}$

(25)

$\frac{d}{dt}\dot{\tilde{p}}=\ddot{p}{\hspace{0.17em}}_d+gz{\hspace{0.17em}}_e-\frac{1}{m}\left[\begin{array}{c}{R}_{13d}\\ R_{23d}\\ R_{33d}\end{array}\right]\frac{{\mu }_1}{R_{33d}}\mathrm{.}$

(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 $\gamma {\hspace{0.17em}}_{p1}\in {ℝ}^3$ we have

${\ddot{p}}_d+gz{\hspace{0.17em}}_e-\frac{1}{m}\left[\begin{array}{c}{R}_{13d}\\ R_{23d}\\ R_{33d}\end{array}\right]\frac{{\mu }_1}{R_{33d}}=\gamma {\hspace{0.17em}}_{p1}\mathrm{.}$

(27)

It is worth noticing that ${\gamma }_{p1}={\left[{\gamma }_{p11}\hspace{0.17em}{\gamma }_{p12}\hspace{0.17em}{\gamma }_{p13}\right]}^T\in {ℝ}^3$ .

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

$\frac{d}{dt}\tilde{p}=\dot{\tilde{p}}\mathrm{,}$

(28)

$\frac{d}{dt}\dot{\tilde{p}}={\mathit{\gamma }}_{p1}\mathrm{.}$

(29)

Therefore, the selection of the PD structure

${\mathit{\gamma }}_{p1}=-K_{pp}\tilde{p}-K_{dp}\dot{\tilde{p}}$

(30)

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

Let us rewrite equation (27) as

$\left[\begin{array}{c}{R}_{13d}\\ R_{23d}\\ R_{33d}\end{array}\right]\frac{{\mu }_1}{m{R}_{33d}}=\ddot{p}{\hspace{0.17em}}_d+gz{\hspace{0.17em}}_e-{\mathit{\gamma }}_{p1}={\mathit{\gamma }}_1\mathrm{,}$

(31)

where ${\gamma }_{p1}=\left[{\gamma }_{p11}\hspace{0.17em}{\gamma }_{p21}\hspace{0.17em}{\gamma }_{p13}\right]\in {ℝ}^3$ . It is possible to look at the equation (31) as a nonlinear algebraic equation system where the unknown variables are ${\phi }_d\left(t\right),\hspace{0.17em}{\phi }_d\left(t\right)$ and ${\mu }_1\left(t\right)\in ℝ$ , while ψ_d(t) is arbitrarily specified. A solution for the nonlinear algebraic equation system (31) is given by

$\begin{array}{l}\mathit{\eta }{\hspace{0.17em}}_d=\left[\begin{array}{c}{\varphi }_d(t)\\ {\theta }_d(t)\\ {\psi }_d(t)\end{array}\right]\\ =\left[\begin{array}{c}\arctan \left(c_{{\theta }_d}\left[s_{{\psi }_d}\frac{m{\gamma }_{11}}{{\mu }_1}-c_{{\psi }_d}\frac{m{\gamma }_{12}}{{\mu }_1}\right]\right)\\ \arctan \left(c_{{\psi }_d}\frac{m{\gamma }_{11}}{{\mu }_1}+s_{{\psi }_d}\frac{m{\gamma }_{12}}{{\mu }_1}\right)\\ {\psi }_d(t)\end{array}\right]\mathrm{,}\end{array}$

(32)

${\mu }_1=m{\gamma }_{13}=m[{\ddot{p}}_{d3}+g-{\gamma }_{p13}]\mathrm{,}$

(33)

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

Notice that the computation of ${\dot{\eta }}_d$ and ${\ddot{\eta }}_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 ${\dot{\Omega }}_d$ which is explicitly given by

$\dot{\Omega }{\hspace{0.17em}}_d=W{(\mathit{\eta })}^{-1}[\ddot{\eta }{\hspace{0.17em}}_d+K_{po}\dot{\tilde{\eta }}]+\left[\frac{d}{dt}\left[W\left(\mathit{\eta }\right){\hspace{0.17em}}^{-1}\right]\right]\left[\dot{\eta }{\hspace{0.17em}}_d+K_{po}\tilde{\eta }\right]\mathrm{.}$

(34)

At the same time, the calculation of ${\dot{\eta }}_d$ and ${\ddot{\eta }}_d$ requires the signals ${\dot{\gamma }}_{p1}$ and ${\ddot{\gamma }}_{p1}$ , which, under definition (30), depend on the unmeasurable signals $\ddot{\stackrel{\sim }{p}}$ and ${\stackrel{\sim }{p}}^{(3)}$ .

In order to respect the assumption that the only measured signals are those of the state ${[p{\hspace{0.17em}}^T\dot{p}{\hspace{0.17em}}^T{\mathit{\eta }}^T\mathbf{\Omega }{\hspace{0.17em}}^T]}^T\in {ℝ}^{12}$ , the following definition of γ _p1 is proposed:

${\dot{\gamma }}_{p1}=-K_{f1}[{\gamma }_{p1}-{\gamma }_{p2}]$

(35)

$\dot{\gamma }{\hspace{0.17em}}_{p2}=-K_{f2}[\gamma {\hspace{0.17em}}_{p2}-[-K_{pp}\tilde{p}-K_{dp}\dot{\tilde{p}}]]\mathrm{,}$

(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_{pp}\stackrel{\sim }{p}(t)-K_{dp}\dot{\stackrel{\sim }{p}}(t)]$ . Therefore, the closed-loop dynamics of the translational error (28)-(29) can be recovered in an exponential way.

The explicit definition of ${\dot{\phi }}_d,{\dot{\theta }}_d,{\ddot{\phi }}_d$ and ${\ddot{\theta }}_d$ are provided in Appendix since they are necessary to compute Ω _d in (18) and ${\dot{\Omega }}_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.

3.3. Summary of the control algorithm

Next, we provide a summary of the equations to implement the proposed trajectory tracking controller. Specifically, we provide the steps to achieve a program for numerical simulation purposes.

• Step 1: Specify the desired signals for translation position p _d (t) and desired yaw angle ψ_d(t), which must be four and two times differentiable, respectively.

• Step 2: Compute the translational errors $\stackrel{\sim }{p}$ and $\dot{\stackrel{\sim }{p}}$

• Step 3: Compute the signal γ_p1 through the linear filter (35)-(36).

• Step 4: Compute µ₁ in (33)

• Step 5: Compute ${\eta }_d,\hspace{0.17em}{\dot{\eta }}_d$ , and ${\ddot{\eta }}_d$ . The explicit expressions for ${\dot{\phi }}_d,\hspace{0.17em}{\dot{\theta }}_d,\hspace{0.17em}{\ddot{\phi }}_d,\hspace{0.17em}{\ddot{\theta }}_d$ are provided in Appendix.

• Step 6: Compute Ω _d and $\dot{\Omega }{\hspace{0.17em}}_d$ in (18) and (34), respectively.

• Step 7: Compute the signal ${\mathit{\mu }}_r=[{\mu }_2{\mu }_3{\mu }_4{]}^T\in {ℝ}^3$ in (16)-(17).

• Step 8: Finally, with µ₁ and µ _r compute the effective control action f in equation (9).

In the next section, a theoretical framework for the exponential stability of the overall closed-loop system is provided. Firstly, the closed-loop system is obtained and then stability will be discussed by using the theory of singularly perturbed systems [17].

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

Consider the following parametrization of the control gains

$K_{pw}=\frac{{\bar{K}}_{pw}}{{\varepsilon }_1}\mathrm{,}\hspace{0.17em}{K}_{iw}=\frac{{\bar{K}}_{iw}}{{\varepsilon }_1}\mathrm{,}$

(37)

$K_{f1}=\frac{{\bar{K}}_{f1}}{{\varepsilon }_1}\mathrm{,}\hspace{0.17em}{K}_{f2}=\frac{{\bar{K}}_{f2}}{{\varepsilon }_1}\mathrm{,}$

(38)

and

$K_{po}=\frac{{\bar{K}}_{po}}{{\varepsilon }_2}\mathrm{,}$

(39)

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

${\varepsilon }_2>>{\varepsilon }_1>0.$

(40)

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

$\frac{d}{dt}\tilde{p}=\dot{\tilde{p}}\mathrm{,}$

(41)

$\frac{d}{dt}\dot{\tilde{p}}=\ddot{p}{\hspace{0.17em}}_d+gz{\hspace{0.17em}}_e-\frac{1}{m}\left[\begin{array}{c}{R}_{13}({\mathit{\eta }}_d\mathrm{,}\tilde{\eta })\\ R_{23}({\mathit{\eta }}_d\mathrm{,}\tilde{\eta })\\ R_{33}({\mathit{\eta }}_d\mathrm{,}\tilde{\eta })\end{array}\right]\frac{{\mu }_1}{R_{33}({\mathit{\eta }}_d\mathrm{,}\tilde{\eta })}\mathrm{,}$

(42)

$\frac{d}{dt}\mathit{\xi }=\tilde{\Omega }\mathrm{,}$

(43)

${\varepsilon }_2\frac{d}{dt}\tilde{\eta }=-{\bar{K}}_{po}\tilde{\eta }+{\varepsilon }_{2}W({\mathit{\eta }}_d\mathrm{,}\tilde{\eta })\stackrel{\sim }{\Omega }\mathrm{,}$

(44)

${\varepsilon }_1\frac{d}{dt}\tilde{\Omega }=-I_f^{-1}[{\varepsilon }_1\mathbf{\Omega }({\varepsilon }_2)\times I_f\tilde{\Omega }+{\bar{K}}_{pw}\tilde{\Omega }+{\bar{K}}_{iw}\mathit{\xi }]\mathrm{,}$

(45)

${\varepsilon }_1\frac{d}{dt}{\gamma }_{p1}=-{\bar{K}}_{f1}[{\gamma }_{p1}-{\gamma }_{p2}]\mathrm{,}$

(46)

${\varepsilon }_1\frac{d}{dt}\gamma {\hspace{0.17em}}_{p2}=-{\bar{K}}_{f2}[\gamma {\hspace{0.17em}}_{p2}-[-K_{pp}\tilde{p}-K_{dp}\dot{\tilde{p}}]]\mathrm{,}$

(47)

where, by virtue of the angular velocity error $\stackrel{\sim }{\Omega }$ in (17), $\Omega \left({\varepsilon }_2\right)={\Omega }_d\left({\varepsilon }_2\right)-\stackrel{\sim }{\Omega }$ , 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 ${\eta }_d\in {ℝ}^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=\left[\begin{array}{c}{x}_s\\ -\\ x_m\end{array}\right]=\left[\begin{array}{c}\tilde{p}\\ \dot{\tilde{p}}\\ \mathit{\xi }\\ -\\ \tilde{\eta }\end{array}\right]$

  are the slow states with respect to

  $z=\left[\begin{array}{c}\tilde{\Omega }\\ {\gamma }_{p1}\\ {\gamma }_{p2}\end{array}\right]\mathrm{,}$

  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{\hspace{0.17em}}_s=\left[\begin{array}{c}\tilde{p}\\ \dot{\tilde{p}}\\ \mathit{\xi }\end{array}\right]$

  are slow states with respect to

$x_m=\tilde{\eta }\mathrm{.}$

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

$\left[\begin{array}{l}{\left[{\stackrel{\sim }{p}}^T{\dot{\stackrel{\sim }{p}}}^T{\xi }^T{\stackrel{\sim }{n}}^T\right]}^T\hfill \\ {\left[{\stackrel{\sim }{\Omega }}^T\hspace{0.17em}{\gamma }_{p1}^T\hspace{0.17em}{\gamma }_{p2}^T\right]}^T\hfill \end{array}\right]=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\mathrm{,}x)\mathrm{,}$

with $x=[{\stackrel{\sim }{p}}^T{\dot{\stackrel{\sim }{p}}}^T\mathit{\xi }{\hspace{0.17em}}^T{\stackrel{\sim }{\eta }}^T{]}^T$ and $z={\left[{\stackrel{\sim }{\Omega }}^T\hspace{0.17em}{\gamma }_{p1}^T\hspace{0.17em}{\gamma }_{p2}^T\right]}^T$ , for ${\varepsilon }_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

$\tilde{\Omega }{\hspace{0.17em}}^{*}=-{\bar{K}}_{pw}^{-1}{\bar{K}}_{iw}\mathit{\xi }\mathrm{,}$

(48)

${\mathit{\gamma }}_{p1}^{*}=-K_{pp}\tilde{p}-K_{dp}\dot{\tilde{p}}\mathrm{,}$

(49)

${\gamma }_{p2}^{*}=-K_{pp}\tilde{p}-K_{dp}\dot{\tilde{p}}\mathrm{.}$

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

$\frac{d}{dt}\tilde{p}=\dot{\tilde{p}}\mathrm{,}$

(51)

$\frac{d}{dt}\dot{\tilde{p}}={\ddot{p}}_d+g{z}_e-\frac{1}{m}\left[\begin{array}{c}{R}_{13}({\mathit{\eta }}_d^{*}\mathrm{,}\tilde{\eta })\\ R_{23}({\mathit{\eta }}_d^{*}\mathrm{,}\tilde{\eta })\\ R_{33}({\mathit{\eta }}_d^{*}\mathrm{,}\tilde{\eta })\end{array}\right]\frac{{\mu }_1}{R_{33}({\mathit{\eta }}_d^{*}\mathrm{,}\tilde{\eta })}\mathrm{,}$

(52)

$\frac{d}{dt}\mathit{\xi }=-{\bar{K}}_{pw}^{-1}{\bar{K}}_{iw}\mathit{\xi }\mathrm{,}$

(53)

${\varepsilon }_2\frac{d}{dt}\tilde{\eta }=-{\bar{K}}_{po}\tilde{\eta }-{\varepsilon }_{2}W({\mathit{\eta }}_d^{*}\mathrm{,}\tilde{\eta }){\bar{K}}_{pw}^{-1}{\bar{K}}_{iw}\mathit{\xi }\mathrm{.}$

(54)

As previously mentioned, the states

$x={\left[x_s^T\hspace{0.17em}|x_m^T\right]}^T=[{\tilde{p}}^T\hspace{0.17em}\dot{\tilde{p}}{\hspace{0.17em}}^T\hspace{0.17em}{\mathit{\xi }}^T|{\tilde{\eta }}^T{]}^T$

are slow with respect to

$z={\left[{\ddot{\Omega }}^T\hspace{0.17em}{\gamma }_{p1}^T\hspace{0.17em}{\gamma }_{p2}^T\right]}^T\mathrm{.}$

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\mathrm{,}x)\mathrm{,}$

is defined. Let us remember that $x=[\stackrel{\sim }{p}{\hspace{0.17em}}^T\dot{\stackrel{\sim }{p}}{\hspace{0.17em}}^T\mathit{\xi }{\hspace{0.17em}}^T\mathit{\eta }{\hspace{0.17em}}^T]$ and $z=[\stackrel{\sim }{\Omega }{\hspace{0.17em}}^T{\gamma }_{p1}^T{\gamma }_{p2}^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

$\frac{d}{d\sigma }y{\hspace{0.17em}}_1=-I_f^{-1}{\bar{K}}_{pw}{y}_1\mathrm{,}$

(55)

$\frac{d}{d\sigma }{y}_2=-{\bar{K}}_{f1}[y_2-y_3]\mathrm{,}$

(56)

$\frac{d}{d\sigma }{y}_3=-{\bar{K}}_{f2}{y}_3\mathrm{,}$

(57)

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

$W(y_1,y_2,y_3)=\frac{1}{2}{y}_1^T{I}_f{y}_1+\beta \frac{1}{2}{y}_2^T{y}_2+\frac{1}{2}{y}_3^T{y}_3,\beta >0,$

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

$\begin{array}{c}\frac{d}{d\sigma }W(y_1\mathrm{,}{y}_2\mathrm{,}{y}_3)\\ =-{\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right]}^T\left[\begin{array}{ccc}{\bar{K}}_{pw}& 0& 0\\ 0& \beta {\bar{K}}_{f1}& -\frac{1}{2}\beta {\bar{K}}_{f1}0\\ 0& -\frac{1}{2}\beta {\bar{K}}_{f1}& {\bar{K}}_{f2}\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right]\mathrm{.}\end{array}$

By invoking Sylvester's theorem, the condition

$\beta <\frac{4{\lambda }_{min}\left\{{\bar{K}}_{f1}\right\}{\lambda }_{min}\left\{{\bar{K}}_{f2}\right\}}{{\lambda }_{Max}{\left\{{\bar{K}}_{f1}\right\}}^2}$

guarantees the inequality

$\frac{d}{d\sigma }W(y_1\mathrm{,}{y}_2\mathrm{,}{y}_3)\le -{\kappa }_{1}W(y_1\mathrm{,}{y}_2\mathrm{,}{y}_3)\mathrm{,}{\kappa }_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\hspace{0.17em}{y}_2{(0)}^T\hspace{0.17em}{y}_3{(0)}^T]}^T\in B_{\rho }$

implies that

${[y_1{(\sigma )}^T\hspace{0.17em}{y}_2{(\sigma )}^T\hspace{0.17em}{y}_3{(\sigma )}^T]}^T\to 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

$\underset{t\to \infty }{{\displaystyle \lim }}\left[\begin{array}{c}x(t)\\ z(t)\end{array}\right]=\underset{t\to \infty }{{\displaystyle \lim }}\left[\begin{array}{c}\tilde{p}(t)\\ \dot{\tilde{p}}(t)\\ \mathit{\xi }(t)\\ \tilde{\eta }(t)\\ \tilde{\Omega }(t)\\ \gamma {\hspace{0.17em}}_{p1}(t)\\ \gamma {\hspace{0.17em}}_{p2}(t)\end{array}\right]=0\mathrm{.}$

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=[\stackrel{\sim }{p}{(t)}^T\dot{\stackrel{\sim }{p}}{(t)}^T\mathit{\xi }{(t)}^T{]}^T$ and the medium time states x _m = ῆ(t) is provided as follows:

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

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 $\left[\begin{array}{c}\hfill {\left[{\stackrel{\sim }{p}}^T\hspace{0.17em}{\dot{\stackrel{\sim }{p}}}^T{\xi }^T\hspace{0.17em}{\stackrel{\sim }{\eta }}^T\right]}^T\hfill \\ \hfill {\left[{\stackrel{\sim }{\Omega }}^T\hspace{0.17em}{\gamma }_{p1}^T\hspace{0.17em}{\gamma }_{p2}^T\right]}^T\hfill \end{array}\right]=0$ .

• 2) The quasi-steady-state solution

  $\bar{z}={\bar{z}}^{*}=x_m=\bar{h}(t\mathrm{,}\hspace{0.17em}{x}_s)$

  is given by the isolated root

  $\tilde{\eta }{\hspace{0.17em}}^{*}=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.,

  $\frac{d}{dt}\tilde{p}=\dot{\tilde{p}}\mathrm{,}$

  (59)
  $\frac{d}{dt}\dot{\tilde{p}}=-K_{pp}\tilde{p}-K_{dp}\dot{\tilde{p}}\mathrm{,}$

  (60)
  $\frac{d}{dt}\mathit{\xi }=-{\bar{K}}_{pw}^{-1}{\bar{K}}_{iw}\mathit{\xi }\mathrm{.}$

  (61)

  Let us consider the Lyapunov function

  $\begin{array}{l}V(\stackrel{\sim }{p},\dot{p},\xi )=\frac{1}{2}{[\alpha \stackrel{\sim }{p}+\dot{\stackrel{\sim }{p}}]}^T[\alpha \stackrel{\sim }{p}+\dot{\stackrel{\sim }{p}}]\\ +\frac{1}{2}\tilde{p}{\hspace{0.17em}}^T[K_{pp}+\alpha K_{dp}-{\alpha }^2{I}_{3\times 3}]\tilde{p}+\frac{1}{2}\mathit{\xi }{\hspace{0.17em}}^T{\bar{K}}_{iw}\mathit{\xi }\mathrm{,}\end{array}$

  (62)

  with α > 0 small enough so that

  $K_{dp}-\alpha I_{3\times 3}>0.$

  The time derivative of $V\left(\stackrel{\sim }{p},\dot{p},\xi \right)$ along the trajectories of the slow dynamics (59)–(61) is given as

  $\begin{array}{l}\dot{V}(\tilde{p}\mathrm{,}\dot{p}\mathrm{,}\mathit{\xi })=-{\dot{\tilde{p}}}^T[K_{dp}-\alpha I_{3\times 3}]\dot{\tilde{p}}\\ -\alpha {\tilde{p}}^T{K}_{pp}\tilde{p}-{\mathit{\xi }}^T{\bar{K}}_{iw}{\bar{K}}_{pw}^{-1}{\bar{K}}_{iw}\mathit{\xi }\mathrm{,}\end{array}$

  (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

  $\dot{V}(\tilde{p}\mathrm{,}\dot{p}\mathrm{,}\mathit{\xi })\le -{\kappa }_{2}V(\tilde{p}\mathrm{,}\dot{p}\mathrm{,}\mathit{\xi })\mathrm{,}$

  with some κ₂ > 0, is satisfied for all ${[{\stackrel{\sim }{p}}^T{\dot{p}}^T\mathit{\xi }{\hspace{0.17em}}^T]}^T\in B_r$ . Therefore, the state-space origin of the slow dynamics (59)–(61) is exponentially stable. In other words,

  $\underset{t\to \infty }{{\displaystyle \lim }}{\left[\tilde{p}\begin{array}{c}{(t)}^T\hspace{0.17em}\dot{\tilde{p}}{(t)}^T\mathit{\xi }{(t)}^T\end{array}\right]}^T=0$

  with exponential convergence rate as time t increases.

• 5) By deriving

  $\bar{y}=\bar{z}-\bar{h}(t\mathrm{,}{x}_s)=x_m-\bar{h}(t\mathrm{,}{x}_s)$

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

$\frac{d}{d\bar{\alpha }}\tilde{\eta }=-{\bar{K}}_{po}\tilde{\eta }$

(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

$\underset{t\to \infty }{{\displaystyle \lim }}{[\tilde{p}{(t)}^T\hspace{0.17em}\dot{\tilde{p}}{(t)}^T\hspace{0.17em}\mathit{\xi }{(t)}^T\hspace{0.17em}\tilde{\eta }(t)]}^T=0\mathrm{.}$

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

${\varepsilon }_2^{*}>{\varepsilon }_2>>{\varepsilon }_1^{*}>{\varepsilon }_1>\mathrm{0,}$

which can be used as a tuning guideline.

Figure 4.

Slow time scale states: the position errors p͂ >(t), translational velocity errors $\dot{\stackrel{\sim }{p}}(t)$ , and integral of the angular velocity error $\mathit{\xi }(t)$ for the different values of the perturbing parameters ε₁ and ε₂. Medium time scale states: orientation error $\stackrel{\sim }{\eta }(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.

Figure 5.

Fast time scale: Time evolution of $y_1(t)=\stackrel{\sim }{\Omega }(t)+{\bar{K}}_{pw}^{-1}{\bar{K}}_{iw}\mathit{\xi }(t){,y}_2(t)={\gamma }_{p1}(t)+K_{pp}\stackrel{\sim }{p}(t)+K_{dp}\dot{\stackrel{\sim }{p}}(t)$ and $y_3(t)={\gamma }_{p2}(t)+K_{pp}\stackrel{\sim }{p}(t)+K_{dp}\dot{\stackrel{\sim }{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)=\left[\begin{array}{c}{x}_d(t)\\ y_d(t)\\ z_d(t)\end{array}\right]={\widehat{a}}_o\left[\begin{array}{c}\widehat{\mathrm{x}}(t)\\ \widehat{\mathrm{y}}(t)\\ \widehat{\mathrm{z}}(t)+1\end{array}\right]\hspace{0.17em}\in {ℝ}^3\mathrm{,}$

(65)

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

$\left[\begin{array}{c}\dot{\widehat{x}}\\ \dot{\widehat{y}}\\ \dot{\widehat{z}}\end{array}\right]={\widehat{\tau }}_o\left[\begin{array}{c}\widehat{y}\left[\widehat{z}-1+{\widehat{x}}^2\right]+\gamma \widehat{x}\\ \widehat{x}\left[3\widehat{z}+1-{\widehat{x}}^2\right]+\gamma \widehat{y}\\ -2\widehat{z}\left[\alpha +\widehat{x}\widehat{y}\right]\end{array}\right]\hspace{0.17em}\in {ℝ}^3\hspace{0.17em},$

(66)

being γ = 0.1 y α = 0.3. Besides, ${\widehat{\tau }}_o=0.2,{\widehat{a}}_o=1$ . Initial conditions are x̂(0) = −0.16, ŷ(0) = 0.1 and ẑ(0) = 0.1. The signals ${\dot{p}}_d,\hspace{0.17em}{\ddot{p}}_d,\hspace{0.17em}{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

${\psi }_d=a_1\sin \left(2\pi {\tau }_1^{-1}t\right)+a_2\sin \left(2\pi {\tau }_2^{-1}t\right)\hspace{0.17em}\in ℝ\hspace{0.17em}\mathrm{,}$

(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,

$\begin{array}{l}{K}_{pw}={\varepsilon }_1^{-1}diag\left\{\mathrm{0.005,0.005,0.075}\right\}\mathrm{,}\\ K_{iw}={\varepsilon }_1^{-1}diag\left\{\mathrm{0.005,0.005,0.075}\right\}\mathrm{,}\\ K_{f1}={\varepsilon }_1^{-1}diag\left\{\mathrm{25,25,25}\right\}\mathrm{,}\\ K_{f2}={\varepsilon }_1^{-1}diag\left\{\mathrm{50,50,50}\right\}\mathrm{,}\hspace{0.17em}\text{and}\hspace{0.17em}\\ K_{po}={\varepsilon }_2^{-1}diag\left\{\mathrm{5,5,5}\right\}\mathrm{,}\hspace{0.17em}\text{where}\hspace{0.17em}\end{array}$

(68)

${\varepsilon }_1=\mathrm{0.1,}\hspace{0.17em}\mathrm{0.4,}\hspace{0.17em}\mathrm{0.7,}\hspace{0.17em}\mathrm{1.0,}$

(69)

and

${\varepsilon }_2=\mathrm{0.5,}\hspace{0.17em}\mathrm{1.5,}\hspace{0.17em}\mathrm{2.5,}\hspace{0.17em}\mathrm{3.5,}\hspace{0.17em}\mathrm{4.5.}$

(70)

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

$\begin{array}{l}{K}_{pp}=diag\left\{\mathrm{2,0.75,10}\right\}\mathrm{,}\\ K_{dp}=\hspace{0.17em}diag\left\{\mathrm{3,1,20}\right\}\mathrm{.}\end{array}$

(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

$\bar{x}(t)=\left[\begin{array}{c}{\bar{x}}_1(t)\\ {\bar{x}}_2(t)\\ {\bar{x}}_3(t)\end{array}\right]=\left[\begin{array}{c}\tilde{p}(t)\\ \dot{\tilde{p}}(t)\\ \mathit{\xi }(t)\end{array}\right]\mathrm{,}$

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

$\bar{z}(t)=\tilde{\eta }(t)$

and the fast time scale states

$y(t)=z(t)-h\left(t\mathrm{,}x(t)\right)\mathrm{,}$

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

$x(t)={\left[\bar{x}{(t)}^T|\bar{z}{(t)}^T\right]}^T:={\left[{\tilde{p}}^T\hspace{0.17em}{\dot{\tilde{p}}}^T\hspace{0.17em}{\xi }^T|{\eta }^T\right]}^T.$

More specifically, the fast states are

$y=\left[\begin{array}{c}{y}_1(t)\\ y_2(t)\\ y_3(t)\end{array}\right]=\left[\begin{array}{c}\tilde{\Omega }(t)+{\bar{K}}_{pw}^{-1}{\bar{K}}_{iw}\mathit{\xi }(t)\\ {\mathit{\gamma }}_{p1}(t)+K_{pp}\tilde{p}(t)+K_{dp}\dot{\tilde{p}}(t)\\ {\mathit{\gamma }}_{p2}(t)+K_{pp}\tilde{p}(t)+K_{dp}\dot{\tilde{p}}(t)\end{array}\right]\mathrm{.}$

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 $\dot{\stackrel{\sim }{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 $\dot{\stackrel{\sim }{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.

Figure 7.

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

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:

• Controller. This block corresponds to the controller law described by equations (9), (16)-(17), (18), (32)-(33) and (35)-(36).

• Saturation. This block limits the forces generated by the propellers. The lower limit is 0 [N] and the upper limit is 1.0 [N]. It is worth noting that this nonlinearity provides only positive control action, which has the physical meaning that the propellers keep the same direction of rotation at all times.

• Low pass filter. This block is represented in the Laplace domain and is referred to as the low pass filter $\frac{a}{s+a}$ , which has the purpose of modelling the actuator dynamics. We have selected a = 20π, which means a response of $\frac{5}{a}$ [sec] for a step input.

• Sensor. Noise is added to all translational and angular measurements. The signals are updated with a rate of 100 Hz.

  The maximum amplitude of the noise in the translational displacement is 0.01 [m]. Similarly, the maximum amplitude of the noise in the angular displacement is 0.05 [rad]. Additionally, the amplitudes of the noise added to the time derivative of the translational displacement and angular displacement are 1.0 [m/sec] and 5.0 [rad/sec], respectively.

• Wind. We have added perturbing forces to the propellers, which emulate the effect of the wind over the quadrotor propellers. The perturbing force applied to each propeller has the direction of the thrust force and is defined as $f_w(t)=\frac{1}{20}\frac{mg}{4}\sin (2\pi t/5)$ [N].

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:

• The first one corresponds to the equations described by the Rabinovich-Fabrikant system in (65)-(66).

• The second one is the helicoidal position trajectory which is usually found in the literature and described by

$p_d(t)=\left[\begin{array}{c}{x}_d(t)\\ y_d(t)\\ z_d(t)\end{array}\right]=\left[\begin{array}{c}1.7\cos (\frac{1}{12}t)\\ 1.7\sin (\frac{1}{12}t)\\ \frac{1}{320}t+1\end{array}\right]\hspace{0.17em}\in {ℝ}^3\mathrm{.}$

(72)

The purpose of using a chaotic-type desired position trajectory as the one in (65)-(66) is to assess the performance for a flight with acrobatic characteristics. The desired yaw angle ψ_d(t) specified for the simulations with disturbances is given in (67).

The gains used in the simulations are in (68)–(71). The different values of the perturbing parameters ε₁ and ε₂ in (69) and (70), respectively, were used to prove numerically the threefold time scale separation.

Figure 9 depicts the actual path p (t) performed by the quadrotor, the norm of the desired translational velocity $\Vert {\dot{p}}_d(t)\Vert$ , and the norm of the actual velocity $\Vert \dot{p}(t)\Vert$ for both simulations. Specifically, the left hand-side plots correspond to the chaotic-type desired position in (65)-(66), while the right hand-side plots stand for the helicoidal position in (72). The perturbing parameters ε₁ = 1 and ε₂ = 1 were used.

The purpose of showing the velocity signals ‖ ṗ (t)‖ and ‖ ṗ _d (t)‖ in the lower plots of Figure 9 is to appreciate the complexity of the manoeuvres that are encoded for each type of desired position p _d (t) and to assess the tracking accuracy of the controller in both simulations. We can see that an acceptable performance is achieved.

By using both types of desired position trajectory p _d (t) and the different values of the perturbing parameters ε₁ and ε₂, Figure 10 describes the performance of the position error x̄ ₁(t) = p͂ (t) (upper plots), the orientation error x̄ (t) = ῆ (t) (middle plots) and the signals ‖ x̄ (t)‖ = ‖ x _s (t)‖ [blue], ‖ z̄ (t)‖ = ‖ x _m (t)‖ [black], and ‖ y (t)‖ [red] (lower plots). Some of the signals have been scaled in magnitude in order to facilitate comparison.

From Figure 10 it can be appreciated that for ε₁ = 0.25 and any value of ε₂ the signal ‖ y (t)‖ converges to zero much faster than x̄(t) = x_s(t) and z̄(t) = x_m(t). For any fixed value of ε₁, it is observed that the states z̄(t) = x_m(t) are the most sensitive to the variations of the parameter ε₂.

The left hand-side plots of Figure 11 describe the forces f_i(t) generated by the propellers for the chaotic-type desired position trajectory p _d (t) in (65)-(66) and all the combinations of perturbing parameters ε₁ and ε₂. Likewise, in the left hand-side plots of Figure 11, the respective results are shown for the helicoidal desired position trajectory p _d (t) in (72). The dashed line in each plot indicates the value $\frac{1}{4}mg$ . Notice that the forces f_i(t) show slight changes from one value of ε₂ to the other. An increase of high frequency contents in the forces f_i(t) is observed for small values of the parameters ε₁. The reason for the high frequency components is the amplification of the noise present in the measurements of position and velocity. Another important observation about the actuator response is that the presence of saturation seems to have little effect on the stability of the closed-loop system. However, for very small values of the perturbing parameters the system may become unstable.

Figure 10.

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). The position error x̄ ₁ (x˜ [blue], y˜ [red] and z˜ [black]) are observed in the upper plots, the orientation error z̄ (φ˜ [blue], θ˜ [red] and ψ˜ [black]) are given in the middle plots, and the norm of the slow time scale states ‖x̄(t)‖ = ‖x_m(t)‖ [blue], the medium time scale states ‖z̄(t) = x_m(t)‖ [black], and the fast time scale states ‖ y (t)‖ [red], are depicted in the lower plots. A darker tone in each colour means a bigger value in the parameter ε₂. The changes in the parameter ε₁are indicated in the third axis of each plot.

The numerical results shown in Figures 9, 10 and 11 allowed verifying that the threefold time scaling is still present in spite of the fact that practical disturbances are affecting the overall closed-loop system.

7. Concluding remarks and further research

The proposed trajectory tracking controller allows choosing any desired translational position p _d (t) at least four times differentiable, and any desired yaw angle φ_d(t) at least two times differentiable. The desired roll ψ_d(t) and the desired pitch θ_d(t) are designed as functions of the quadrotor states and the output of a second order linear filter.

The overall closed-loop system stability is proven by using the framework of the theory of singularly perturbed systems. This analysis suggested that by using high-gain angular velocity control, the coupling between the translational error dynamics and the angular error dynamics is diminished.

Numerical simulations confirmed that the motion control objective is satisfied with the proposed scheme. In addition, the different time scales introduced by the controller were numerically verified under ideal conditions and under the situation that the system is affected by different types of perturbations that appear in a practical implementation.

A possible disadvantage of the proposed approach is the use of Euler angles to represent the orientation dynamics in equations (14)-(15), which implies that the matrix W( η ) loses rank for θ = ±π/2(2n + 1), being n = 1, 2, · · ·,. A way to avoid the problem of singularities is to chose a non-minimal representation of the orientation such as the unit quaternion. See references [37], [38] and [39] for examples of the unit quaternion in robot control. Application of the unit quaternion to the proposed two-loop control scheme requires a full redesign of the controller and therefore is left for further research.

Figure 11.

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). Using different combinations of the parameters ε₁ and ε₂, the plots show the time evolution of the forces f (t) generated by the propellers. The vertical axis of each plot shows the applied force. Notice that the forces are limited as 0 ≤ f_i(t) ≤ 1.0 [N]. The increase in the colour tone corresponds to the increase in the parameter ε₂. The changes in the parameter ε₁ are indicated in the third axis of each plot.

Another problem that is left for further research is the estimation of ε₁ and ε₂ in Proposition 1 and 2, respectively, which can achieved by using a Lyapunov function.