Recently,vertical takeoff and landing (VTOL) unmanned aerial vehicles (UAV) have attracted increasing interest in researches and applications in both military and civil society, such as rescue in disasters,unmanned inspection,and road traffic supervision. The motivation also comes from academic research institutes,since it can be used as low cost testbeds for robotics studies. However,the VTOL UAV model,which has been widely investigated in most works,is known as a class of underactuated system with nonholonomic constraints of second order. According to the necessity of Brockett,there is no gloss or time invariant controller that can stabilize the underactuated system to the equilibrium point^{[1]}. Hence,new methodology should be investigated for this kind of systems.
Trajectory tracking control of VTOL UAV is a challenging work due to its coupling property,external disturbances,system uncertainties, etc. Several inspired approaches have been investigated,such as backstepping control^{[2, 3, 4]},sliding mode control^{[2, 5, 6, 7]}, feedback linearization^{[5]},model predictive control^{[8]}, neural networks^{[9]},fuzzy control^{[10]},observer based control^{[11, 12]},etc. However,there are still some prominent problems to be considered and resolved.
1) The attitude representation and desired attitude extraction. In most previous works ^{[2, 3, 4, 5, 6, 7, 8, 9, 10, 13]},Euler angles are used to represent the attitude of rigid body. However,the simplified kinematics is often used as
$ \dot \phi = \omega _x ,\dot \theta = \omega _y ,\dot \psi = \omega _z,$  (1) 
where $\phi ,\theta ,\psi$ denote the roll,pitch,and yaw angles. $\omega _x ,\omega _y ,\omega _z$ are the angular velocity of the rigid body,respectively. It is pointed out that the original kinematics of Euler angles is described as^{[14]}
$ \left[{\begin{array}{*{20}c} {\dot \phi } \\ {\dot \theta } \\ {\dot \psi } \\ \end{array}} \right] = \left[{\begin{array}{*{20}c} 1 & {\sin \phi \tan \theta } & {\cos \phi \tan \theta } \\ 0 & {\cos \phi } & {  \sin \phi } \\ 0 & {\sin \phi \sec \theta } & {\cos \phi \sec \theta } \\ \end{array}} \right]\left[{\begin{array}{*{20}c} {\omega _x } \\ {\omega _y } \\ {\omega _z } \\ \end{array}} \right]. $  (2) 
The existence of transcendental functions in (2) makes it difficult to design a control strategy. Noticing that (1) is a simplified form of (2) with the assumption that the rigid body rotates only in one direction at a time,and the roll/pitch angle changes when the pitch/roll angle equals to $0^{\circ}$. However,the assumption above is an ideal instance,which is infeasible. Meanwhile,the simplification will decrease the control accuracy. Furthermore,we find that the system model based on Euler angles is not available when pitch angle $\theta = \pm \pi / 2$. Especially,considering the error of sensors and calculation,both attitude estimation and control algorithm based on Euler angles cannot work near the state of $\theta = \pm \pi / 2$. Moreover,most works provided the desired attitude,angular velocity,and angular acceleration directly by the position controller. Only ^{[12, 15]} present the analytical solution of the desired attitude information based on quaternion.
2) Stability of hierarchical control structure. Concerning with the hierarchial control strategy,the position and attitude controllers can be designed separately for translational and rotational subsystems,respectively. Although the above strategy can be introduced for controller design,the stability should be analyzed based on the overall closedloop system,since the attitude's tracking is an asymptotical procedure,which makes the attitude error between the actual and desired one a necessary concern in the analysis. However,in ^{[2, 3, 4, 5, 7, 8, 9, 10, 13]},stability is only analyzed for each subsystem.
3) Dynamics of actuators and its influence on the closedloop stability. In ^{[12, 15, 16]} and ^{[11, 17]},cascade theory and singular perturbation theory are used to acquire the stability of the closedloop system. In ^{[18]},the relationship is given between actuators and control input of a VTOL UAV. However,dynamics of actuators is never considered.
4) Controller design with internal uncertainties and external disturbances. Several researches related to controller design against uncertainties have been studied based on sliding mode control^{[2, 5, 6, 7]},neutral networks^{[9]},fuzzy systems^{[10]},adaptive backstepping control^{[19]},disturbance observer^{[16]},etc. However,the chartering of SMC,convergence rate of weights in neutral networks and fuzzy systems will limit the applications of these methods in practical. The adaptive algorithm can only deal with the external disturbances^{[19]}. The disturbance observer is adopted to deal with the uncertainties^{[20]},however,the saturation of the actuators is not considered.
Based on the above review,trajectory tracking control of a VTOL UAV is explored is this paper,taking both uncertainties and actuators$'$ dynamics into account. The error model of VTOL UAV based on trajectory tracking task is established based on modified Rodrigues parameters (MRPs),based on which analytical expression of desired attitude information is given. Then,considering the dynamics of actuators,the overall system is divided into three subsystems according to their timescale properties,based on which a hierarchical control structure is presented. Thereafter, antiwindup controllers against system uncertainties are proposed based on translational and rotational subsystems,respectively. Stability of the overall closedloop system is analyzed based on singular perturbation theory. In summary,the main contributions of the proposed control strategy are presented as follows:
1) A hierarchical control structure is proposed due to the cascade property between translational and rotational subsystems. Meanwhile,the analytical solution of desired attitude information based on MRPs is given.
2) A modified disturbance rejection controller is proposed for disturbance rejection performance as well as the input saturation of actuators.
3) The singular perturbation theory is employed with consideration of the actuators$'$ dynamics,based on which the strictly Lyapunov stability conclusion is achieved.
The rest of this paper is organized as follows. In Section Ⅱ,the trajectory tracking error model is established based on MRPs,and the analytical solution of desired attitude information is given. In Section Ⅲ,the overall system is divided into three subsystems, based on which a hierarchical strategy is introduced. In Section Ⅳ, antiwindup controllers are proposed based on translational and rotational subsystems,respectively. In Section Ⅴ,singular perturbation theory is introduced to analyzed the stability of the overall closedloop system. Simulations are presented in Section Ⅴ to verify the effectiveness of the proposed control strategy, followed by the conclusions in Section Ⅵ.
Ⅱ SYSTEM MODEL AND PROBLEM FORMULATION A System ModelThere are totally three coordinates used in this paper,earth frame ${\mathcal F}_e$,bodyfixed frame ${\mathcal F}_b$,and orientated frame ${\mathcal F}_d$. We choose MRPs to represent the attitude. MRPs are described as a threedimension vector without restrictions,which is defined as ${\pmb \sigma} = {\pmb r}\tan (\alpha/4)$,where ${\pmb r}$ and $\alpha$ represent the unit vector of rotational axis and rotation angle of the rigid body, respectively. Due to the definition of MRPs,its kinematics is given as:
$ \dot {\pmb \sigma} = G(\pmb \sigma){\pmb \omega},$  (3) 
where ${\pmb \omega} \in {\bf R}^3$ denotes the angular velocity of the VTOL UAV. The matrix $G({\pmb \sigma})$ is given as
$ G({\pmb\sigma}) = \frac{1}{2}\left( {\frac{{1 {\pmb\sigma}^{\rm T}{\pmb \sigma}}}{2} I_3 + [{{\pmb \sigma} \times }] + {\pmb \sigma} {\pmb \sigma} ^{\rm T} } \right),$  (4) 
We consider the VTOL UAV as a rigid body without deformation,and the system model is described as
$ \left\{\begin{aligned} &\dot {\pmb \xi} = {\pmb v},\\ &m\dot {\pmb v} = {\pmb e}_3 mg  R{\pmb e}_3 T + {\pmb d}_1,\\ &\dot {\pmb \sigma} = G({\pmb \sigma}){\pmb \omega} ,\\ &J\dot {\pmb \omega} = {\pmb \omega} \times J{\pmb \omega} + {\pmb \tau} + {\pmb d}_2,\\ \end{aligned} \right. $  (5) 
$ R = I_3  \frac{{4( {1  {\pmb \sigma} ^{\rm T} {\pmb \sigma} } )}}{{( {1 + {\pmb \sigma} ^{\rm T} {\pmb \sigma} } )^2 }}[{{\pmb \sigma} \times }] + \frac{8}{{( {1 + {\pmb \sigma} ^{\rm T} {\pmb \sigma} } )^2 }}[{{\pmb \sigma} \times }]^2. $  (6) 
Remark 1. $T$ and $\pmb \tau$ are the resulted aerodynamic force and moment described in the bodyfixed coordinate,which lead to motion of the VTOL UAV. Different UAVs have different types of actuators,whose aerodynamic characteristics are also different. Without loss of generality,we consider the aerodynamic force and moment as the input of the flight control system. In the next subsection,the dynamics of the thrust and torque caused by the actuators are also taken into account.
B Problem FormulationThe trajectory tracking problem of a VTOL UAV is investigated,and the objective in this work is to design a control thrust $T_{\rm d}$ and control torque ${\pmb \tau}_{\rm d}$,which enable the VTOL UAV to track a desired trajectory quickly and accurately. Define the system errors as: position error $\tilde {\pmb \xi} = {\pmb \xi}  {\pmb \xi}_{\rm d}$,and velocity error $\tilde {\pmb v}={\pmb v}  \dot {\pmb \xi}_{\rm d}$. $\tilde {\pmb \sigma}$, $\tilde {\pmb \omega}$ in (7) are errors of MRPs and angular velocity given as:
$ \tilde {\pmb \sigma} = {\pmb \sigma} \oplus {\pmb \sigma}_{\rm d}^{  1} ,\tilde {\pmb \omega} = {\pmb \omega}  \tilde R{\pmb \omega}_{\rm d},$  (7) 
where ${\pmb \sigma}_{\rm d}^{1}$ is the inverse of ${\pmb \sigma}_{\rm d}$,known as ${\pmb \sigma}_{\rm d}^{  1} =  {\pmb \sigma}_{\rm d}$,and $\tilde R = R R_{\rm d}^{\rm T}$ is known as the error of attitude transition matrix. Operator $ \oplus$ denotes the production of MRPs,which is described in (8) with two MRPs variables of ${\pmb \sigma}_1$ and ${\pmb \sigma}_2$
$ {\pmb \sigma}_1 \oplus {\pmb \sigma}_2 = \frac{{( {1  \ {{\pmb \sigma}_2 } \^2 } ){\pmb \sigma}_1 + ( {1  \ {{\pmb \sigma} _1 } \^2 } ){\pmb \sigma}_2  2{\pmb \sigma}_1 \times {\pmb \sigma}_2 }}{{1 + \ {{\pmb \sigma}_2 } \^2 \ {{\pmb \sigma}_1 } \^2  2{\pmb \sigma}_2^{\rm T}{\pmb \sigma}_1 }},$  (8) 
The analytical solution of ${\pmb \sigma}_{\rm d}$,${\pmb \omega}_{\rm d}$,$\dot {\pmb \omega}_{\rm d}$ is shown in Section III. And Lemma 1 holds for MRPs error.
Lemma 1 ^{[22]}.If the attitude variable pairs $( {{\pmb \sigma},{\pmb \omega}})$ and $( {{\pmb \sigma}_{\rm d},{\pmb \omega}_{\rm d} } )$ both satisfy MRPs kinematics in (3), their relative attitude variable pair also satisfies (3).
Denoting the nominal values of mass and rotational inertia as $m_0$ and $J_0$,then their errors are given as
$ \left\{ \begin{aligned} & \Delta m = m  m_0,\\ & \Delta J = J  J_0. \\ \end{aligned} \right. $  (9) 
The thrust and torque inputs of VTOL UAV are obtained by the servo systems,such as motors,flapping angles,control surfaces,etc., which may affect the stability of the overall closedloop system. Assuming the controller of the actuators make the error dynamics of thrust and torque satisfy:
$ \dot {\tilde T} = \frac{K_T}{t_T}\tilde T,\quad \dot {\tilde {\pmb \tau}} = \frac{K_\tau}{t_{\pmb \tau}}\tilde {\pmb \tau},$  (10) 
Based on descriptions above,the system error model is represented as
$ \left\{ \begin{aligned} &\dot{\tilde {\pmb \xi}} = \tilde {\pmb v},\\ &\dot{\tilde {\pmb v}} = g{\pmb e}_3  \frac{R {\pmb e}_3 T}{m_0} + {\pmb d}'_1  \ddot {\pmb \xi}_{\rm d},\\ &\dot{\tilde {\pmb \sigma}} = G( {\tilde {\pmb \sigma}} )\tilde {\pmb \omega},\\ & \dot{\tilde {\pmb \omega}} = J_0^{  1} ( L({\pmb \omega}) J_0^{*} + {\pmb \tau}) + {\pmb d'_2}  ( \tilde R\dot {\pmb \omega}_{\rm d}  [{\tilde {\pmb \omega} \times }]\tilde R{\pmb \omega}_{\rm d} ),\\ & \dot {\tilde T} = \frac{K_T}{t_T}\tilde T,~~\dot {\tilde {\pmb \tau}} = \frac{K_{\pmb \tau}}{t_{\pmb \tau}}\tilde {\pmb \tau}. \\ \end{aligned} \right. $  (11) 
The compound disturbances on the system dynamics are given as
$ \left\{ \begin{aligned} & {\pmb d}'_1 \!=\!\frac{{\pmb d}_1}{m_0}  \frac{\Delta m}{m_0} (\dot{\tilde {\pmb v}}  g {\pmb e}_3 + \ddot {\pmb \xi}_{\rm d}),\\ & {\pmb d}'_2 \!=\! J_0^{  1}\left[{\pmb d}_2 \!\!\!\! \Delta J \dot{\tilde {\pmb \omega}} \!\!\!\! L({\pmb \omega}) \Delta J^{*}  \Delta J( \tilde R\dot {\pmb \omega}_{\rm d}\!\!\!\! [{\tilde {\pmb \omega} \times }] \tilde R{\pmb \omega}_{\rm d} ) \right],\\ \end{aligned} \right. $  (12) 
where $J^*$ denotes the vector form of the diagonal elements of $J$ and,for ${\pmb \omega} = [{\begin{array}{*{20}c} \omega_1 & \omega_2 & \omega_3 \end{array}}]^{\rm T}$,we have:
$ L({\pmb \omega}) = \left[{\begin{array}{*{20}c} 0 & \omega_2\omega_3 & \omega_2\omega_3 \\ \omega_3\omega_1 & 0 & \omega_3\omega_1 \\ \omega_1\omega_2 & \omega_1\omega_2 & 0 \\ \end{array}} \right]. $  (13) 
A traditional method to design guidance and control strategies in aeronautics is assuming that the controller will enable the rotational dynamics to converge faster than the translational dynamics by using an attitude controller with higher gains.
From the VTOL UAV model shown above,we know that the attitude error will converge asymptotically after the convergence of actuators. The position error will converge asymptotically after the convergence of both attitude error and actuators. According to the convergent speed of the different parts of the overall system, we regard the translational subsystem as slow subsystem, rotational subsystem as fast subsystem,and actuators$'$ dynamics as ultrafast subsystem. The timescale property of each subsystem is shown in Fig. 1.
Download:


Fig. 1. Timescale property of each system. 
The singular perturbation theory can be used in the controller design and stability analysis based on the multitimescale properties of VTOL UAV system. Two timescale factors $\varepsilon_1$ and $\varepsilon_2$ are introduced to formalize the timescale separation. Introducing the following notations:
$ \tilde {\pmb \omega}' = \varepsilon_1 \tilde {\pmb \omega},{\pmb \tau}'_{\rm d} = \varepsilon_1{\pmb \tau}_{\rm d},K'_T = \varepsilon_2 K_T,K'_{\pmb \tau} = \varepsilon_2K_{\pmb \tau},$ 
we finally get the error model of VTOL UAV as
$ \Sigma_1 \left\{ \begin{aligned} &\dot{\tilde {\pmb \xi}} = \tilde {\pmb v},\\ &\dot{\tilde {\pmb v}} = g{\pmb e}_3  \frac{1}{m_0}(R_{\rm d} {\pmb e}_3 T_{\rm d}  {\pmb f}_1) + {\pmb d}'_1  \ddot {\pmb \xi}_{\rm d},\\ \end{aligned} \right. $  (14) 
$ \Sigma_2 \left\{ \begin{aligned} &\varepsilon_1 \dot{\tilde {\pmb \sigma}} = G( {\tilde {\pmb \sigma}} )\tilde {\pmb \omega},\\ &\varepsilon_1 \dot{\tilde {\pmb \omega}} = \varepsilon_1 J_0^{  1} (  L({\pmb \omega}) \Delta J^{*} + {\pmb f}_2) + \varepsilon_1{\pmb d}'_2 + J^{1}_0{\pmb \tau}'_{\rm d}\\ & \varepsilon_1( {\tilde R\dot {\pmb \omega}_{\rm d}  \tilde R[{\tilde {\pmb \omega} \times }] {\pmb \omega}_{\rm d} } ),\\ \end{aligned} \right. $  (15) 
$ \Sigma_3 \left\{ \begin{aligned} & \varepsilon_2 \dot {\tilde T} = \frac{K_T}{t_T}\tilde T ,\\ & \varepsilon_2 \dot {\tilde {\pmb \tau}} = \frac{K_{\pmb \tau}}{t_{\pmb \tau}}\tilde {\pmb \tau},\\ \end{aligned} \right. $  (16) 
where the coupling terms are defined as ${\pmb f}_1 = ( { R_{\rm d} {\pmb e}_3 \tilde T}) + ( {I_3  \tilde R} )( {R_{\rm d} {\pmb e}_3 T_{\rm d}})  ( {I_3  \tilde R} )( {R_{\rm d} {\pmb e}_3 \tilde T}) $,${\pmb f}_2 = {\tilde {\pmb \tau}}$.
The purpose of these two timescale factors is to adjust the gain of the controller for each subsystem,whose convergence speed will be changed correspondingly. From the system model in (5),we find that the transition matrix $R$ and control thrust $T$ will affect the translational motion of VTOL UAV. Notice from (6) that the transition matrix $R$ in terms of MRPs can also be regarded as the output of the rotational subsystem. This can be considered as the cascade property of VTOL UAV. A practical hierarchical strategy is introduced to implement the control system. Consequently, translational and rotational controllers can be designed separately. The translational controller is firstly designed to extract the desired thrust $T_{\rm d}$ and attitude matrix $R_{\rm d}$. The desired attitude information can enable a VTOL UAV to track the desired trajectory. Thereafter,the desired torque vector ${\pmb \tau}_{\rm d}$ is determined by the rotational controller with the desired attitude matrix. At last,$T_{\rm d}$ and ${\pmb \tau}_{\rm d}$ can be treated as the input for the actuators to implement the whole control system. The following Condition is assumed in the controller design. In the procedure of controller design,the subsystems/subsystem whose convergent speeds/speed are/is higher than the corresponding subsystem to be controlled are/is already converged. However,the stability should also be analyzed based on the original error model,and this assumption is only used in the procedure of controller design. The diagram of hierarchical control strategy is shown in Fig. 2.
Download:


Fig. 2. Control structure of the system. 
Since the desired attitude information ${\pmb \sigma}_{\rm d}$, ${\pmb \omega}_{\rm d}$,$\dot {\pmb \omega}_{\rm d}$ is determined by the virtual controller $R_{\rm d}$,we present the analytical solution of these information.
Theorem 1. By introducing the notation ${\pmb \delta} = [ {\begin{array}{*{20}c} \delta_1 & \delta_2 & \delta_3 \end{array}}]^{\rm{T}} \triangleq R_{\rm d}{\pmb e}_3 T_{\rm d}$, it is always possible to extract the desired attitude as:
$ \left\{ \begin{aligned} \eta_{\rm d} &= \sqrt{\frac{\delta_3}{2T_{\rm d}} + \frac{1}{2}},\\ {\pmb \sigma}_{\rm d} &= \frac{1}{2T_{\rm d}\eta_{\rm d}(1+\eta_{\rm d})} \left[{\begin{array}{*{20}c} \delta_2 \\ \delta_1 \\ 0 \end{array}} \right]. \end{aligned} \right. $  (17) 
We assume the virtual controller $R_{\rm d}{\pmb e}_3 T_{\rm d}$ is differentiable,the desired angular velocity is given as
${\omega _{\rm{d}}} = \Gamma (\delta )\mathop \delta \limits^. ,$  (18) 
$ \Gamma (\delta ) = \frac{1}{{T_{\rm{d}}^2\gamma }}\left[{\begin{array}{*{20}{c}} {  {\delta _1}{\delta _2}}&{  \delta _2^2 + {T_{\rm{d}}}\gamma }&{{\delta _2}\gamma }\\ {\delta _1^2  {T_{\rm{d}}}\gamma }&{  {\delta _1}{\delta _2}}&{  {\delta _1}\gamma }\\ {{\delta _2}{T_{\rm{d}}}}&{  {\delta _1}{T_{\rm{d}}}}&0 \end{array}} \right],$  (19) 
Then,the desired angular acceleration is described as
$ \dot {\pmb \omega}_{\rm d} = {\bar \Gamma}({\pmb \delta},\dot {\pmb \delta})\dot {\pmb \delta} + \Gamma({\pmb \delta}) \ddot {\pmb \delta}. $  (20) 
$ \begin{aligned} {\bar \Gamma}({\pmb \delta},\dot {\pmb \delta}) = \frac{(2{\pmb \delta}^{\rm T}\dot {\pmb \delta}\gamma + T_{\rm d}{\pmb \delta}^{\rm T}\dot {\pmb \delta} + T^2_{\rm d}\dot \delta_3)}{T^2_{\rm d}\gamma}\Gamma({\pmb \delta}) + \frac{1}{T^2_{\rm d}\gamma} \Gamma', \end{aligned} $  (21) 
where $\Gamma'$ is a matrix with the size of 3$\times$3.
Proof. We notate ${\pmb \sigma}_{\rm d} = [ {\begin{array}{*{20}c} \sigma_{\rm d1} & \sigma_{\rm d2} & \sigma_{\rm d3} \end{array}}]^{\rm{T}}$. The vector ${\pmb \delta}$ is described from the definition of MRPs as
${\pmb \delta} = \frac{T_{\rm d}}{(1+{\pmb \sigma}^{\rm T}_{\rm d}{\pmb \sigma}_{\rm d})^2} \left[{\begin{array}{*{20}c} 8\sigma_{\rm d1}\sigma_{\rm d3}  4\sigma_{\rm d2} (1{\pmb \sigma}^{\rm T}_{\rm d}{\pmb \sigma}_{\rm d}) \\ 8\sigma_{\rm d2}\sigma_{\rm d3} + 4\sigma_{\rm d1}(1{\pmb \sigma}^{\rm T}_{\rm d}{\pmb \sigma}_{\rm d}) \\ 4(\sigma^2_{\rm d1}\sigma^2_{\rm d2}+\sigma^2_{\rm d3}) + (1{\pmb \sigma}^{\rm T}_{\rm d}{\pmb \sigma}_{\rm d})^2 \end{array}} \right]. $  (22) 
Notice that there is a constraint $\\pmb \delta\ = 1$,hence, only two degreesoffreedom of rotation can be determined by this vector. We assume $\sigma_{\rm d3} = 0$ to calculate $\sigma_{\rm d1}$ and $\sigma_{\rm d2}$. Since $\{\pmb \sigma}_{\rm d}\ \leq 1$,we have:
$ \begin{aligned} \eta_{\rm d} = \sqrt{\frac{\delta_3}{2T_{\rm d}} + \frac{1}{2}} = \frac{1\{\pmb \sigma}^2_{\rm d}\}{1+\{\pmb \sigma}^2_{\rm d}\} = \frac{1{\pmb \sigma}^{\rm T}_{\rm d}{\pmb \sigma}_{\rm d}}{1+{\pmb \sigma}^{\rm T}_{\rm d}{\pmb \sigma}_{\rm d}}, \end{aligned} $  (23) 
and
$ \left \{ \begin{aligned} \delta_1 &= \frac{  4T_{\rm d}\sigma_{\rm d2}(1{\pmb \sigma}^{\rm T}_{\rm d} {\pmb \sigma}_{\rm d})}{(1+{\pmb \sigma}^{\rm T}_{\rm d}{\pmb \sigma}_{\rm d})^2},\\ \delta_2 &= \frac{ 4T_{\rm d}\sigma_{\rm d1}(1{\pmb \sigma}^{\rm T}_{\rm d} {\pmb \sigma}_{\rm d})}{(1+{\pmb \sigma}^{\rm T}_{\rm d}{\pmb \sigma}_{\rm d})^2}. \end{aligned} \right. $  (24) 
Consequently,we can easily prove that (17) holds. From (4) and the following equations
$ \left \{ \begin{aligned} &{\dot T}_{\rm d} = \frac{1}{T_{\rm d}}{\pmb \delta}^{\rm T}\dot {\pmb \delta},\\ &\gamma^2 + \delta^2_1 + \delta^2_2 = 2T_{\rm d}\gamma,\\ \end{aligned} \right. $  (25) 
we can obtain (18)$\sim$(21).
Ⅳ.CONTROLLER DESIGNBased on the above error model,a hierarchical control scheme is present to exploit the cascade property. The control design steps can be summarized as follows:
Step 1. The translational controller is designed based on subsystem $\Sigma_1$ under the assumption that ${\pmb f}_1 = 0$. A linear extended state observer (ESO)^{[23]} is introduced to estimate and compensate the compound disturbances. The desired attitude information is extracted by Theorem 1 as the input of subsystem $\Sigma_2$.
Step 2. The rotational controller is designed based on subsystem $\Sigma_2$ under the assumption that ${\pmb f}_2 = 0$. A linear ESO is also used to suppress the attitude error caused by the compound disturbances.
Step 3. Stability of the overall system is analyzed based on Lyapunov analysis,taking the dynamics of actuators into account.
A Translational Controller DesignWe notate the derivative of ${\pmb d}'_1$ as ${\pmb h}_1(t)$. Then,the secondorder linear ESO for translational subsystem $\Sigma_1$ is
$ \left\{ \begin{aligned} & \dot {\hat {\pmb z}}_1 = g{\pmb e}_3  \frac{R {\pmb e}_3 T}{m_0} + {\hat {\pmb z}}_2  \ddot {\pmb \xi}_{\rm d} + g_1 ({\tilde {\pmb v}}  {\hat {\pmb z}}_1),\\ & \dot {\hat {\pmb z}}_2 = g_2 ({\tilde {\pmb v}}  {\hat {\pmb z}}_1),\\ \end{aligned} \right. $  (26) 
where $g_1$ and $g_2$ are positive constants to be selected.
Transforming (26) to the frequencydomain using the Laplace transform,and substitute (14) into (26),we get:
$ \left\{ \begin{aligned} & s{\hat {\pmb z}}_1 = ({\hat {\pmb z}}_2  {\pmb d}'_1) + s\tilde {\pmb v} + g_1 ({\tilde {\pmb v}}  {\hat {\pmb z}}_1),\\ & s{\hat {\pmb z}}_2 = g_2 ({\tilde {\pmb v}}  {\hat {\pmb z}}_1),\\ \end{aligned} \right. $  (27) 
where $s$ is the Laplace operator.
From (27),we finally get
$ \hat {\pmb z}_2 = \frac{g_2}{s^2 + g_1 s + g_2}{\pmb d}'_1. $  (28) 
$g_1$ and $g_2$ should satisfy that the polynomial $s^2 + g_1 s + g_2$ is Hurwitz. Here,we simply choose $g_1 = 2 \omega_0,g_2 = \omega^2_0$. The ESO views both internal uncertainties and external disturbances as the extended state to be estimated and compensated in the controller. Hence,the ESO can reject the influence caused by both internal uncertainties and external disturbances.
The backstepping technique and an auxiliary observer are introduced to design the trajectory tracking controller. We firstly introduce the following variables:
$ \left\{ \begin{aligned} &{\pmb e}_1 = \tilde {\pmb \xi}  {\pmb \alpha},\\ &{\pmb e}_2 = \tilde {\pmb v}  \dot {\pmb \alpha} + k_1 {\pmb e}_1. \\ \end{aligned} \right. $  (29) 
Then,the control input is designed as
$ \left \{ \begin{aligned} &T_{\rm d} = m_0\left\{\pmb e}_3g  \ddot {\pmb \xi}_{\rm d} + k_{\alpha1} \tanh ({\pmb \alpha}) + k_{\alpha2}\tanh (\dot {\pmb \alpha}) \right\,\\ &R_{\rm d}{\pmb e}_3 = \frac{m_0 \left[{\pmb e}_3g  \ddot {\pmb \xi}_{\rm d} + k_{\alpha1}\tanh ({\pmb \alpha}) + k_{\alpha2}\tanh (\dot {\pmb \alpha}) \right]}{T_{\rm d}},\\ \end{aligned} \right. $  (30) 
where $k_{\alpha1}$,$k_{\alpha2}$ are strictly positive matrices, and an auxiliary observer similar to ^{[12]} is given as
$ \begin{aligned} \ddot {\pmb \alpha} =& k_{\alpha1}\tanh ({\pmb \alpha})  k_{\alpha2}\tanh (\dot {\pmb \alpha}) + \hat {\pmb z}_2 + (k_1+k_2){\pmb e}_2 +\\ & (1k^2_1){\pmb e}_1, \end{aligned} $  (31) 
where $k_1$,$k_2$ are positive matrices to be selected.
It is easy to verify that the thrust input is bounded as
$ \T_{\rm d}\ \leq m_0 (\\ddot {\pmb \xi}_{\rm d}\ + k_{\alpha1} + k_{\alpha2}+g),$  (32) 
and for a candidate Lyapunov function $V_1 = \frac{1}{2}({\pmb e}^{\rm T}_1{\pmb e}_1 + {\pmb e}^{\rm T}_2{\pmb e}_2)$,its derivative is given as
$\begin{array}{l} {{\dot V}_1} \le  {\lambda _{\min }}({k_1}){e_1}{^2}  {\lambda _{\min }}({k_2}){e_2}{^2} + \\ \;\;\;\;\;{e_2}({\widetilde d_1} + \frac{{{f_1}}}{{{m_0}}}), \end{array}$  (33) 
where $\tilde {\pmb d}_1 \triangleq {\pmb d}'_1  \hat {\pmb z}_2$.
The derivative of the control input $R_{\rm d}{\pmb e}_3T_{\rm d}$ is described as follows:
$ \left \{\begin{aligned} \dot {\pmb \delta} &= m_0[{\pmb \xi}^{\cdots}_{\rm d} + k_{\alpha1} \frac{{d} \tanh({\pmb \alpha})}{{d} ({\pmb \alpha})} \dot {\pmb \alpha} + k_{\alpha2}\frac{{d} \tanh (\dot {\pmb \alpha})}{{d} (\dot {\pmb \alpha})} \ddot {\pmb \alpha}],\\ \ddot {\pmb \delta} &= m_0[{\pmb \xi}^{(4)}_{\rm d} + k_{\alpha1} \frac{{d} \tanh({\pmb \alpha})}{{d} ({\pmb \alpha})} \ddot {\pmb \alpha} + k_{\alpha1}\frac{{d}^2 \tanh({\pmb \alpha})} {{d} ({\pmb \alpha})^2} \dot {\pmb \alpha}^2 + \\ &k_{\alpha2}\frac{{d} \tanh (\dot {\pmb \alpha})}{{d} (\dot {\pmb \alpha})}\dddot {\pmb \alpha} + k_{\alpha2}\frac{{d}^2 \tanh (\dot {\pmb \alpha})}{{d} (\dot {\pmb \alpha})^2}\ddot {\pmb \alpha}^2],\\ \end{aligned} \right. $  (34) 
with $\dddot {\pmb \alpha}$ being
$ \begin{aligned} \dddot {\pmb \alpha} =& k_{\alpha1}\frac{{d} \tanh({\pmb \alpha})} {{d} ({\pmb \alpha})} \dot {\pmb \alpha}  k_{\alpha2}\frac{{d} \tanh (\dot {\pmb \alpha})}{{d} (\dot {\pmb \alpha})}\ddot {\pmb \alpha} + g_2(\tilde {\pmb v}  \hat {\pmb z}_1)+ \\ & (k_1 + k_2)({\pmb e}_3g  \frac{R{\pmb e}_3T_{\rm d}}{m_0} + \hat {\pmb z}_2  \ddot {\pmb \xi}_{\rm d}  \ddot {\pmb \alpha} +\\ & k_1{\pmb e}_2  k^2_1{\pmb e}^2_1) + (1k^2_2)({\pmb e}_2  k_1{\pmb e}_1). \end{aligned} $  (35) 
By introducing the notation $B = G^{1}(\tilde {\pmb \sigma})$,we have $\tilde {\pmb \omega} = B\dot {\tilde {\pmb \sigma}}$. Then we get
$ \varepsilon_1 \ddot{\tilde {\pmb \sigma}} = \varepsilon_1( J_0 B)^{  1}{\pmb f}_2 + \varepsilon_1{\pmb d}'_2 + (J_0 B)^{1}{\pmb \tau}'_{\rm d}. $  (36) 
In order to extract the bounded controller,we rewrite the compound disturbances as
$ \begin{align} & {{{{d}'}}_{2}}={{({{J}_{0}}B)}^{1}}[{{d}_{2}}\Delta J\dot{\tilde{\omega }}\omega \times J\omega  \\ & \ \ \ \ \ \ \ J(\dot{B}\dot{\tilde{\sigma }}\tilde{R}[\tilde{\omega }\times]{{\omega }_{\text{d}}}+\tilde{R}{{{\dot{\omega }}}_{\text{d}}})]. \\ \end{align} $  (37) 
Then,the system dynamics can be rewritten as
$ \ddot {\tilde {\pmb \sigma}} = \frac{(J_0B)^{1}{\pmb \tau}'_{\rm d}} {\varepsilon_1} + {\pmb d}'_2,\dot {\pmb d}'_2 = {\pmb h}_2(t),$  (38) 
where ${\pmb h}_2(t)$ denotes the derivative of the compound disturbances ${\pmb d}'_2$.
Then,the secondorder linear ESO for rotational subsystem $\Sigma_2$ is proposed as
$ \left\{ \begin{aligned} & \dot {\hat {\pmb z}}_3 = \frac{{\pmb \tau}'_{\rm d}}{\varepsilon_1} + {\hat {\pmb z}}_4 + g_3 ({\tilde {\pmb \omega}}  {\hat {\pmb z}}_3),\\ & \dot {\hat {\pmb z}}_4 = g_4 ({\tilde {\pmb \omega}}  {\hat {\pmb z}}_3),\\ \end{aligned} \right. $  (39) 
where $g_3$ and $g_4$ are positive constants to be selected.
Transforming (39) to the frequencydomain using the Laplace transform,and substitute (15) into (39),we have:
$ \left\{ \begin{aligned} & s{\hat {\pmb z}}_3 = ({\hat {\pmb z}}_4  {\pmb d}') + s\dot {\tilde {\pmb \sigma}} + g_3 (\dot {\tilde {\pmb \sigma}}  {\hat {\pmb z}}_3),\\ & s{\hat {\pmb z}}_4 = g_4 (\dot {\tilde {\pmb \sigma}}  {\hat {\pmb z}}_3),\\ \end{aligned} \right. $  (40) 
where $s$ is the Laplace operator.
From (40),we finally get
$ \hat {\pmb z}_4 = \frac{g_4}{s^2 + g_3 s + g_4}{\pmb d}'_2. $  (41) 
$g_3$ and $g_4$ should satisfy that the polynomial $s^2 + g_3 s + g_4$ is Hurwitz. Here,we simply choose $g_3 = 2 \omega_1,g_4 = \omega^2_1$.
To design the attitude tracking controller,we introduce the following variables:
$ \left\{ \begin{aligned} &{\pmb e}_3 = \tilde {\pmb \sigma}  {\pmb \beta},\\ &{\pmb e}_4 = \dot {\tilde {\pmb \sigma}}'  \varepsilon_1\dot {\pmb \beta} + k_3 {{\pmb e}_3 }. \\ \end{aligned} \right. $  (42) 
Then,the control input ${\pmb \tau}'_{\rm d}$ is given as:
$ {\pmb \tau}'_{\rm d} = J_0B \left[ k_{\beta1}\tanh({\pmb \beta})  k_{\beta2}\tanh(\dot {\pmb \beta})\right],$  (43) 
where $k_{\beta1}$,$k_{\beta2}$ are strictly positive matrices, and the observer is described as
$ \begin{aligned} \ddot {\pmb \beta} =& \frac{1}{\varepsilon_1}[k_{\beta1}\tanh ({\pmb \beta})  k_{\beta2}\tanh (\dot {\pmb \beta})\!+\!\varepsilon_1\hat {\pmb z}_4 + \frac{(k_3 + k_4)}{\varepsilon_1}{\pmb e}_4 + \\ & \frac{(1k^2_3)}{\varepsilon_1}{\pmb e}_3], \end{aligned} $  (44) 
where $k_3$,$k_4$ are positive matrices to be selected.
The control torque of the attitude tracking problem is finally described as
$ \begin{aligned} {\pmb \tau}_{\rm d} = \frac{J_0B}{\varepsilon_1}[ k_{\beta1} \tanh({\pmb \beta})  k_{\beta2}\tanh(\dot {\pmb \beta})]. \end{aligned} $  (45) 
For a candidate Lyapunov function $V_2 = \frac{1}{2}({\pmb e}^{\rm T}_3{\pmb e}_3 + {\pmb e}^{\rm T}_4{\pmb e}_4)$,its derivative is as follows:
$ \begin{aligned} \dot V_2 \leq & \frac{\lambda_{\min}(k_3)}{\varepsilon_1}\{\pmb e}_3\^2  \frac{\lambda_{\min}(k_4)}{\varepsilon_1}\{\pmb e}_4\^2+ \\ &\varepsilon_1\{\pmb e}_4\\left(\\tilde {\pmb d}_2\ + \frac{\{\pmb f_2}\}{\lambda_{\max}(J_0B)}\right),\\ \end{aligned} $  (46) 
where $\tilde {\pmb d}_2 \triangleq {\pmb d}'_2  \hat {\pmb z}_4$.
Remark 2. From (32),the translational controller is bounded. That is,the output of thrust is saturated. From (45) and the definition of $\tanh$ function,it is clear that the output of torque is saturated.
Ⅴ. STABILITY ANALYSISTheorem 2. Given the error model of a VTOL UAV for trajectory tracking problem in (11),with the compound disturbances shown in (12). Let the thrust input $T_{\rm d}$ and desired attitude information given by (30) and Theorem 1,respectively,with a linear ESO and an auxiliary observer proposed in (26) and (31). Then,let the torque input ${\pmb \tau}_{\rm d}$ in (45),with a linear ESO and an auxiliary observer designed in (39) and (4). There exist the timescale factors such that the proposed control strategy can stabilize the system asymptotically.
Proof Consider the candidate Lyapunov functions $V_1$ and $V_2$ defined in Section IV,we define a new Lyapunov function as: $V = V_1 + V_2 + \frac{1}{2}\tilde T^2 + \frac{1}{2}\tilde {\pmb \tau}^{\rm T} \tilde {\pmb \tau}$. Notice that $\ {I_3  \tilde R} \ = 2 {\sin \frac{\alpha}{2}}  \leq 4\ {\tilde {\pmb \sigma}} \ \leq 4(\{\pmb e}_3\ + \{\pmb \beta}\)$,$\ {{\pmb E}_1} \ \leq \ {\tilde {\pmb \xi}} \ + \ {\tilde {\pmb v}}\$,$\ {R_{\rm d}{\pmb e}_3 \tilde T} \ = \ {\tilde T} \$.
As shown later,$\{\pmb \beta}\$ can converge to zero asymptotically. Then,the derivative of $V$ is shown as
$ \begin{aligned} \dot V \leq &  a_1(\{\pmb e}_1\^2 + \{\pmb e}_2\^2)  \frac{a_2}{\varepsilon_1}(\{\pmb e}_3\^2 + \{\pmb e}_4\^2)  \frac{a_3}{\varepsilon_2}{\tilde T}^2 \\ & \frac{a_4}{\varepsilon_2}\{\tilde {\pmb \tau}}\^2 + \frac{1}{m_0}\{\pmb e}_2\(5{\tilde T} + 4T_{\max}\\tilde {\pmb \sigma}\) + \\ & \frac{\varepsilon_1}{\lambda_{\max}(J_0B)}\{\pmb e}_4\ \\tilde {\pmb \tau}\ + \{\pmb e}_2\\\tilde {\pmb d}_1\ + \varepsilon_1\{\pmb e}_4\\\tilde {\pmb d}_2\, \end{aligned} $  (47) 
where
$ \left\{ \begin{aligned} a_1 &= \min\{\lambda_{\min}(k_1),\lambda_{\min}(k_2)\},\\ a_2 &= \min\{\lambda_{\min}(k_3),\lambda_{\min}(k_4)\},\\ a_3 &= \frac{K_T}{t_T},a_4 = \frac{K_{\pmb \tau}}{t_{\pmb \tau}}. \\ \end{aligned} \right. $  (48) 
Define the general error vector of the VTOL UAV as: ${\pmb E}^{\rm T} = \left[{\begin{array}{*{20}c} {\pmb e}^{\rm T}_1 & {\pmb e}^{\rm T}_2 & {\pmb e}^{\rm T}_3 & {\pmb e}^{\rm T}_4 & \tilde T & \tilde {\pmb \tau}^{\rm T}\\ \end{array} } \right]$,and notice that:
$ \left\{ \begin{aligned} b_1 =& \frac{2\lambda_{\max } (P_1)T_{\max}}{m_0},\\ b_2 =& \frac{5\lambda_{\max} (P_1)}{2m_0},\\ b_3 =& \frac{1}{2\lambda_{\max}(J_0B)}. \\ \end{aligned} \right. $  (49) 
Then,the derivative of $V_2$ can be rewritten as
$ \dot V \leq  {\pmb E}^{\rm T}\Gamma{\pmb E} + \{\pmb e}_2\\\tilde {\pmb d}_1\ + \varepsilon_1\{\pmb e}_4\\\tilde {\pmb d}_2\,$  (50) 
where
$ \Gamma = \left[{\begin{array}{*{20}c} a_1 & 0 & 0 & 0 & 0 & 0 \\ 0 & a_1 & b_1 & 0 &  b_2 & 0 \\ 0 & b_1 & \frac{a_2}{\varepsilon_1} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{a_2}{\varepsilon_1} & 0 & \varepsilon_1b_3 \\ 0 & b_2 & 0 & 0 & \frac{a_3}{\varepsilon_2} & 0 \\ 0 & 0 & 0 & \varepsilon_1b_3 & 0 & \frac{a_4}{\varepsilon_2} \end{array} } \right]. $  (51) 
Since $a_1$ is a positive constant,the first two minors of matrix $\Gamma$ is positive. Then,we should find the scopes of $\varepsilon_1,\varepsilon_2$ such that the third to sixth minors are positive. For simplicity,let $\Gamma_i$ be the matrix$'$s minor of $i$. Then,the third to sixth minors are given as
$ \left \{\begin{aligned} \det(A_3) = &  a_1 b_1^2 + \frac{a^2_1a_2}{\varepsilon_1},\\ \det(A_4) = & \frac{a^2_1a^2_2}{\varepsilon^2_1}  \frac{a_1a_2b^2_1}{\varepsilon_1} ,\\ \det(A_5) = & \frac{a^2_1a^2_2a_3}{\varepsilon^2_1\varepsilon_2}  \frac{a_1a^2_2b^2_2}{\varepsilon^2_1}  \frac{a_1a_2a_3b^2_1}{\varepsilon_1\varepsilon_2},\\ \det(A_6) = & a_1a_2b^2_2b^2_3\varepsilon_1 + \frac{a^2_1a^2_2a_3a_4} {\varepsilon^2_1\varepsilon^2_2}  \frac{a^2_1a_2a_3b^2_3\varepsilon_1}{\varepsilon_2} \\ &\frac{a_1a^2_2a_4b^2_2}{\varepsilon^2_1\varepsilon_2} + \frac{a_1a_3b^2_1b^2_3\varepsilon^2_1}{\varepsilon_2}  \frac{a_1a_2a_3a_4b^2_1}{\varepsilon_1\varepsilon^2_2}. \\ \end{aligned} \right. $  (52) 
From $\det(A_3) > 0$ and $\det (A_4) > 0$,we have:
$ \varepsilon _1 < \frac{a_1 a_2}{b_1^2}. $  (53) 
Since $\det(A_5) > 0$,it follows that:
$ \varepsilon_2 \leq \frac{a_1 a_2 a_3  a_3 b_1^2 \varepsilon_1}{a_2 b_2^2}. $  (54) 
If $\det(A_6) > 0$,the following should be satisfied:
$ A\varepsilon_2^2 + B\varepsilon_2 + C > 0,$  (55) 
where $A = a_1 a_2 b_2^2 b_3^2\varepsilon^3_1$,$B = a_1a_3b_1^2b_3^2\varepsilon^4_1  a_1^2a_2a_3b_3^2\varepsilon^3_1  a_1a_2^2a_4b_2^2$,$C = a_1^2a_2^2a_3a_4  a_1a_2a_3a_4b_1^2\varepsilon_1$.
Considering that $\varepsilon_1 < \frac{a_1 a_2}{b_1^2}$,we have $A > 0$,$C > 0$,$B < 0$. $B^24AC = (a_1a_3b_1^2b_3^2\varepsilon^4_1  a_1^2a_2a_3b_3^2\varepsilon^3_1 + a_1a_2^2a_4b_2^2)^2 > 0$.
Define the variable $\varepsilon_2^ * = \frac{ B  \sqrt{B^2  4AC}}{2A}$,and consider that $B^2  4AC < B^2$. Hence, $\varepsilon _2^*$ is positive. At this time,$\varepsilon _2 < \varepsilon _2^ *$ can make $\det(A_6) > 0$. $\varepsilon_2$ is finally determined as
$ \varepsilon_2 \leq \min \left \{ \frac{a_1 a_2 a_3  a_3 b_1^2 \varepsilon_1}{a_2 b_2^2},\varepsilon_2^* \right \}. $  (56) 
If the timescale factors $\varepsilon_1$ and $\varepsilon_2$ are selected based on the above requirements,matrix $\Gamma$ is positive. The unforced system is exponentially stable,that is, the system is inputtostate stable with the input $\\tilde {\pmb d}'_1\$ and $\\tilde {\pmb d}'_2\$.
From (28) and (41),we know that $\hat {\pmb z}_2$ and $\hat {\pmb z}_4$ can converge to ${\pmb d}'_1$ and ${\pmb d}'_2$ asymptotically,that is,$\tilde {\pmb d}'_1$ and $\tilde {\pmb d}'_2$ can converge to zero asymptotically. Then,from Lemma 4.7 of ^{[24]},we know that the cascade overall system is asymptotically stable. Hence,$\mathop {\lim }\nolimits_{t \rightarrow \infty} {\pmb e}_1 = \mathop {\lim }\nolimits_{t \rightarrow \infty} {\pmb e}_2 = \mathop {\lim }\nolimits_{t \rightarrow \infty} {\pmb e}_3 = \mathop {\lim }\nolimits_{t \rightarrow \infty} {\pmb e}_4 = \mathop {\lim }\nolimits_{t \rightarrow \infty} {\tilde T} = \mathop {\lim }\nolimits_{t \rightarrow \infty} \tilde {\pmb \tau} = 0$.
From the descriptions above and the results in ^{[25]},the auxiliary observers in (31) and (44) are asymptotically stable. Consequently, $\mathop {\lim }\nolimits_{t \rightarrow \infty} {\pmb \alpha} = \mathop {\lim }\nolimits_{t \rightarrow \infty} \dot {\pmb \alpha} = \mathop {\lim }\nolimits_{t \rightarrow \infty} {\pmb \beta} = \mathop {\lim }\nolimits_{t \rightarrow \infty} \dot {\pmb \beta} = 0$. Noticing that ${\pmb e}_1$ to ${\pmb e}_4$ are linear diffeomorphism of $\tilde {\pmb \xi}$,$\tilde {\pmb v}$,$\tilde {\pmb \sigma}$,$\tilde {\pmb \omega}$,${\pmb \alpha}$,$\dot {\pmb \alpha}$,${\pmb \beta}$ and $\dot {\pmb \beta}$,hence,we have the following conclusion:
$ \mathop {\lim }\limits_{t \rightarrow \infty} \tilde {\pmb \xi} = \mathop {\lim }\limits_{t \rightarrow \infty} \tilde {\pmb v} = \mathop {\lim }\limits_{t \rightarrow \infty} \tilde {\pmb \sigma} = \mathop {\lim }\limits_{t \rightarrow \infty} \tilde {\pmb \omega} = 0. $  (57) 
Simulations are shown to illustrate the effectiveness of the proposed control strategy. We consider a VTOL UAV model with the parameters being set as: $m = 4$kg,$J_x = J_y = 0.08{\rm kg} \cdot {\rm m}^2$,and $J_z = 0.14\,{\rm kg} \cdot {\rm m}^2$. The initial condition is given as: ${\pmb \xi}(t_0) = [ {\begin{array}{*{20}c} 2 & 3 & 5 \end{array}}]^{\rm T}$m,${\pmb v}(t_0) = [{\begin{array}{*{20}c} 0 & 0 & 0 \end{array}}]^{\rm T}$m/s,${\pmb \sigma}(t_0) = [ {\begin{array}{*{20}c} 0 & 0 & 0 \end{array}}]^{\rm T}$,${\pmb \omega}(t_0) = [{\begin{array}{*{20}c} 0 & 0 & 0 \end{array}}]^{\rm T}$rad/s. The parameters of the controller are given as follows: $k_2 = 3.5$,$k_3 = 2$,$k_4 = 0.2$,$k_{\alpha1} = k_{\alpha2} = 1.5$,$k_{\beta1} = k_{\beta2} = 2.5$, $\varepsilon_1 = 0.1$,$\varepsilon_2 = 0.05$.
Tracking of a spiral rising trajectory with the existence of perturbation of parameters and unknown disturbance is accomplished in Matlab/Simulink. The desired trajectory is as follows:
$ {\pmb \xi}_{\rm d} = [{\begin{array}{*{20}c} 0.5t & 5\sin(\frac{\pi t}{25}) & 5\cos(\frac{\pi t}{25})  2 \end{array}} ]^{\rm T} {\rm m}. $  (58) 
The external disturbances acting on the translational and rotational dynamics are given as:
$ {{d}_{1}}=\left[\begin{matrix} {} & 0.2\sin (\pi t)+1.5\cos (\frac{\pi t}{10}) \\ {} & 0.1\sin (\pi t)+1\cos (\frac{\pi t}{10}) \\ {} & 0.1\sin (\pi t)+2.5\cos (\frac{\pi t}{10}) \\ \end{matrix} \right]\text{N} $  (59) 
$ {{d}_{2}}=\left[\begin{matrix} {} & 0.1\sin (\pi t)+0.1\cos (\frac{\pi t}{10}) \\ {} & 0.1\sin (\pi t)+0.1\cos (\frac{\pi t}{10}) \\ {} & 0.1\sin (\pi t)+0.1\cos (\frac{\pi t}{10}) \\ \end{matrix} \right]\text{N}\cdot \text{m},$  (60) 
Simulation results are illustrated in Figs.3~7. The trajectory tracking effect of the VTOL UAV is illustrated in Fig. 3. The tracking errors of position,velocity,MRPs and angular velocity are shown in Figs.4 and 5,while Fig. 6 shows the estimation effects of linear ESO for both translational and rotational subsystems. The changing tendency of roll,pitch and yaw angles during the trajectory tracking are indicated in Fig. 7.
Download:


Fig. 3. Trajectory tracking effect. 
Download:


Fig. 4. Tracking error of position and velocity. 
Download:


Fig. 5. Tracking error of MRPs and angular velocity. 
Download:


Fig. 6. Disturbance estimation performance. 
Download:


Fig. 7. Equivalent control effect of Euler angles. 
In the simulation,we find that with the existence of both external disturbances and internal uncertainties,the proposed controller can enable the VTOL UAV to track a timevarying trajectory quickly and accurately. The linear ESO can estimate the disturbances and compensate for them in the control scheme to improve the control performance,whereas the proposed controller can enable the VTOL UAV to track a desired trajectory effectively. We also carry out the adaptive backstepping method in ^{[19]} for comparison. In Table I,the rootmeansquare (RMS) error of the proposed control strategy is compared with that of adaptive backstepping method. It is shown that with timevarying disturbances and internal uncertainties,the control performance of the adaptive backstepping is not as well as our proposed control strategy. It is shown in Fig. 6 that the linear ESO estimates the disturbance accurately,and the estimation error converges quickly. Figs.4 and 5 also show that the proposed control strategy have good tracking performance.
Ⅶ. CONCLUSIONIn this paper,the trajectory tracking control of a VTOL UAV is investigated. The MRPs based system error model is established and the hierarchical control strategy is introduced based on the timescale property. Then,a practical disturbance rejection controller is proposed with linear ESO for both translational and rotational subsystems,respectively. The auxiliary observer is implemented to guarantee the boundedness of the control output. At last,stability conclusion of the overall system is given based on singular perturbation theory. Simulation results verify that the proposed control strategy can successfully enable the VTOL UAV to track a desired trajectory. The designed linear ESO can also estimate the compound disturbances caused by both external and internal uncertainties for higher accuracy of tracking.
[1]  Brockett R W. Asymptotic stability and feedback stabilization. Differential Geometric Control Theory. Boston:Birkhäuser, 1983. 181191 
[2]  Bouabdallah S, Siegwart R. Backstepping and slidingmode techniques applied to an indoor micro quadrotor. In:Proceedings of the 2005 IEEE International Conference on Robotics and Automation. Barcelona, Spain:IEEE, 2005.22472252 
[3]  Mian A A, Wang D B. Modeling and backsteppingbased nonlinear control strategy for a 6 DOF quadrotor helicopter. Chinese Journal of Aeronautics, 2008, 21(3):261268 
[4]  Zuo Z. Trajectory tracking control design with commandfiltered compensation for a quadrotor. IET Control Theory and Applications, 2010, 4(11):23432355 
[5]  Lee D, Kim H J, Sastry S. Feedback linearization vs. adaptive sliding mode control for a quadrotor helicopter. International Journal of Control, Automation and Systems, 2009, 7(3):419428 
[6]  Xu R, Ö zgüner Ü. Sliding mode control of a class of underactuated systems. Automatica, 2008, 44(1):233241 
[7]  Efe M Ö. Battery power loss compensated fractional order sliding mode control of a quadrotor UAV. Asian Journal of Control, 2012, 14(2):413425 
[8]  Raffo G V, Ortega M G, Rubio F R. An integral predictive/nonlinear control structure for a quadrotor helicopter. Automatica, 2010, 46(1):2939 
[9]  Efe M O. Neural network assisted computationally simple Pl^{λ} D^{μ} control of a quadrotor UAV. IEEE Transactions on Industrial Informatics, 2011, 7(2):354361 
[10]  Zemalache K M, Maaref H. Controlling a drone:comparison between a based model method and a fuzzy inference system. Applied Soft Computing, 2009, 9(2):553562 
[11]  Bertrand S, Guénard N, Hamel T, PietLahanier H, Eck L. A hierarchical controller for miniature VTOL UAVs:design and stability analysis using singular perturbation theory. Control Engineering Practice, 2011, 19(10):10991008 
[12]  Abdessameud A, Tayebi A. Global trajectory tracking control of VTOLUAVs without linear velocity measurements. Automatica, 2010, 46(6):10531059 
[13]  Nicol C, Macnab C J B, RamirezSerrano A. Robust adaptive control of a quadrotor helicopter. Mechatronics, 2011, 21(6):927938 
[14]  Ginsberg J. Engineering Dynamics. Cambridge:Cambridge University Press, 2008. 
[15]  Abdessameud A, Tayebi A. Formation control of VTOL Unmanned Aerial Vehicles with communication delays. Automatica, 2011, 47(11):23832394 
[16]  Wang L, Jia H M. The trajectory tracking problem of quadrotor UAV:global stability analysis and control design based on the cascade theory. Asian Journal of Control, 2014, 16(2):574588 
[17]  Esteban S, Gordillo F, Aracil J. Threetime scale singular perturbation control and stability analysis for an autonomous helicopter on a platform. International Journal of Robust and Nonlinear Control, 2013, 23(12):13601392 
[18]  Michael N, Mellinger D, Lindsey Q, Kumar V. The GRASP multiple microUAV testbed. IEEE Robotics and Automation Magazine, 2010, 17(3):5665 
[19]  Cabecinhas D, Cunha R, Silvestre C. A nonlinear quadrotor trajectory tracking controller with disturbance rejection. Control Engineering Practice, 2014, 16:110 
[20]  Wang L, Su J B. Global trajectory tracking of VTOL UAV based on disturbance rejection control. In:Proceedings of the 32nd Chinese Control Conference. Xi'an, China:IEEE, 2013. 42704275 
[21]  Tsiotras P. Further passivity results for the attitude control problem. IEEE Transactions on Automatic Control, 1998, 43(11):15971600 
[22]  Cong B L, Liu X D, Chen Z. Distributed attitude synchronization of formation flying via consensusbased virtual structure. Acta Astronautica, 2011, 68(1112):19731986 
[23]  Zheng Q, Dong L L, Lee D H, Gao Z Q. Active disturbance rejection control for MEMS Gyroscopes. IEEE Transactions on Control Systems Technology, 2009, 17(6):14321438 
[24]  Khalil H K. Nonlinear Systems (Third edition). Upper Saddle River, New Jersey:Prentice Hall, 2002. 
[25]  OlfatiSaber R. Global configuration stabilization for the VTOL aircraft with strong input coupling. IEEE Transactions on Automatic Control, 2002, 47(11):19491952 