Available online at www.sciencedirect.com
Journal of the Franklin Institute 351 (2014) 5494–5510 www.elsevier.com/locate/jfranklin
Robust motion control of quadrotors Hao Liua,b, Jianxiang Xic,n, Yisheng Zhongb a
School of Astronautics, Beihang University, Beijing 100191, PR China Department of Automation, TNList, Tsinghua University, Beijing 100084, PR China c High-Tech Institute of Xi'an, Xi'an 710025, PR China
b
Received 15 June 2014; received in revised form 22 August 2014; accepted 7 October 2014 Available online 30 October 2014
Abstract In this paper, the robust motion control problem is investigated for quadrotors. The proposed controller includes two parts: an attitude controller and a position controller. Both the attitude and position controllers include a nominal controller and a robust compensator. The robust compensators are introduced to restrain the influence of uncertainties such as nonlinear dynamics, coupling, parametric uncertainties, and external disturbances in the rotational and translational dynamics. It is proven that the position tracking errors are ultimately bounded and the boundaries can be specified by choosing controller parameters. Experimental results on the quadrotor demonstrate the effectiveness of the robust control method. & 2014 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.
1. Introduction Unmanned rotorcrafts have attracted increasing attention in recent years due to their various functions [1–3]. These aerial robotic vehicles can carry sufficient payload to implement numerous indoor and outdoor missions and the improvements on the battery capacity and other technology are rapidly increasing their commercial applications [4]. The rotorcrafts possess six degrees of freedom (6DOF) and have the ability to move in a three-dimensional space. Therefore, compared to wheeled vehicles, the aerial vehicles have advantages on mapping, indoor exploration, surveillance, inspection, etc. The unmanned rotorcraft discussed in this paper n
Corresponding author. Tel./fax: þ86 29 84744111. E-mail addresses:
[email protected] (H. Liu),
[email protected] (J. Xi),
[email protected] (Y. Zhong). http://dx.doi.org/10.1016/j.jfranklin.2014.10.003 0016-0032/& 2014 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
Eby
5495
T2
T1 Ebz
Ebx
T3 T4 Ez Ey
Ex
Fig. 1. The schematic of the quadrotor.
is a quadrotor aerial vehicle as shown in Fig. 1. This kind of rotorcraft has some advantages over the regular rotorcraft due to its four independent rotors. By varying the angular speeds of the rotors, the quadrotor can change its aerodynamic forces and moments without the need of a swashplate to control the blade pitch angles. Besides, the need for a tail rotor is eliminated because the rotors are paired to rotate in opposite directions and thereby the resultant moment in the conventional rotorcraft is canceled [5]. In fact, quadrotors have become ideal experimental platforms for the aerial vehicle control research worldwide (see, e.g., [4–6]). Proportional-integral-derivative (PID) controllers were applied to achieve the automatic flight of the quadrotors in [4–6]. Then, nested saturation algorithms based on the Lyapunov analysis and flatness-based feedback control approach were discussed in [7,8], respectively, to achieve the trajectory tracking for four-rotor mini rotorcrafts. However, the quadrotor is a multiple-input multiple-output (MIMO) underactuated system containing various uncertainties such as parametric uncertainties, unmodeled uncertainties, nonlinear dynamics, coupling, and external disturbances. Although the control methods described in [4–8] are promising, the robust tracking performance cannot be guaranteed, because the uncertainty problem is not discussed in their controller design. Moreover, many efforts have been devoted to the robust attitude controller design to reduce the influence of the uncertainties in the rotational dynamics. Fuzzy controller [9], sliding-mode based feedback controller [10], PD2 feedback controller [11], quaternion-based feedback controller [12], and switching model predictive controller [13] were applied for the rotorcrafts to stabilize the attitude. However, for the methods described in [11,12], accurate nonlinear models of the quadrotors were needed. In [13], the considered uncertainties were limited to external time-invariant disturbances. In [10], the proposed sliding mode observer required a transitional process to estimate the uncertainties and the transitional properties cannot be specified. Furthermore, the robust attitude control methods depicted in [9–13] focused on the uncertainties in the rotational dynamics, whereas how to restrain the uncertainties in the translational dynamics was not further discussed. The two-stage control scheme is a common method for the quadrotors [4]. By this approach, the attitude and position controllers can be designed separately to achieve the attitude and position control. This method relies on the time-scale separation assumption that the closed-loop translational dynamics converges much slower than the closed-loop rotational dynamics. Twoloop nonlinear control approaches by input–output linearization and command-filtered backstepping method were applied in [14,15] to achieve the motion control for the quadrotors. Their approaches mainly focused on the uncertainties in the rotational dynamics, whereas the disturbances in the translational dynamics were ignored. In addition, in [16], a partial state feedback controller by the singular perturbation theory was designed for a miniature quadrotor
5496
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
without discussing the effects of uncertainties. In [17–19], data-driven ways were discussed to obtain the linear models of modern complicated industrial systems and model-free control approaches were proposed. In [20,21], new kinds of switching controller and robust controller were designed for a class of slowly switched systems and continuous-time systems, respectively. In [22–24], new robust control methods were proposed for the robot manipulator control. In this paper, a robust control scheme is proposed for quadrotors. This control method is developed based on the robust compensation technique (RCT) as defined in [25], which has already been successfully applied to the attitude regulation of the 3DOF helicopter (see, e.g., [26,27]). However, the motion control problem for the quadrotor is more challenging because this vehicle possesses 6DOF. Furthermore, the rotorcraft is underactuated in the longitudinal and latitudinal directions; that is, the pitch and roll angles are required to be changed properly in order to achieve the longitudinal and the latitudinal motion control. In this paper, a robust controller is designed with an attitude controller and a position controller to address the rotational and translational control problems, respectively. The position controller stabilizes the height of the quadrotor and generates the desired pitch and roll angle references based on the position tracking errors in the longitudinal and latitudinal directions. Besides, the attitude controller achieves the desired reference tracking of the pitch, roll, and yaw angles. Both the attitude and position controllers consist of a nominal controller and a robust compensator to restrain the influence of the uncertainties. Compared to previous studies on the motion control problem of the quadrotor, the uncertainties including parametric perturbations, coupling, nonlinear dynamics, and external disturbances involved in the rotational and translational dynamics can be restrained. The rotorcraft can achieve the motion control with high precision and the position tracking errors of the quadrotor are proven to be ultimately bounded with a priori set boundary. Another advantage is that this method does not depend on the time-scale separation assumption and the attitude tracking errors are considered as uncertainties rather than being ignored in the discussions of the translational dynamics. In addition, this control method results in a linear time-invariant and decentralized controller and it is easily implementable in practical applications. Actually, since the models of conventional helicopters and fixed-wing aircrafts share some special features of the quadrotor model, the proposed nonlinear robust control method can be applied to their robust controller design. Furthermore, the proposed robust control method can also be extended to the closed-loop control system design for uncertain high-order systems. The remaining parts of this paper are organized as follows: the nonlinear mathematical model of the quadrotor is described in Section 2; the design procedure of the robust motion controller is presented in Section 3; Section 4 analyzes the robust tracking properties of the closed-loop system; Experimental results are shown in Section 5; Section 6 draws the concluding remarks. 2. Nonlinear model of quadrotor As depicted in Fig. 1, the quadrotor has four control inputs: the rotational speeds of the four rotors, and has 6DOF: the pitch, roll, and yaw angles, the longitudinal and latitudinal positions, and the height. The longitudinal (front and back) rotors rotate clockwise whereas the latitudinal (left and right) rotors spin in the opposite direction. The rotorcraft will pitch and hence move longitudinally if the angular speed of the front rotor is increased and that of the back rotor is decreased. A roll and a latitudinal motion of the quadrotor results from changing the rotational velocities of the left and right rotors in the opposite directions. Increasing (reducing) the angular speeds of the longitudinal rotors while reducing (increasing) those of the latitudinal rotors causes
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
5497
the yaw movement. The sum of thrusts produced from each rotor results in the up-down (height) motion of the rotorcraft. Let α ¼ Ex E y E z denote an earth-fixed inertial frame and β ¼ Ebx Eby Ebz a frame attached to the quadrotor with origin in the center of mass of the vehicle as shown in Fig. 1. The T vector ξ ¼ ξx ξy ξz indicates the position of the origin of the body-fixed frame β with respect to the inertial frame α. ξx denotes the longitudinal position of the vehicle, ξy the latitudinal position, and ξz the height. Moreover, let η ¼ ½θ ϕ ψ T indicate the three Euler angles, the pitch angle θ, the roll angle ϕ, and the yaw angle ψ, which determine the rotation matrix R from the inertial frame α to the body-fixed frame β: 2 3 C θ C ψ C ψ Sϕ Sθ C ϕ S ψ Sϕ Sψ þ C ϕ C ψ Sθ 6 7 R ¼ 4 C θ Sψ C ϕ C ψ þ Sϕ Sθ Sψ C ϕ Sθ Sψ C ψ Sϕ 5 ; Sθ
C θ Sϕ
Cϕ Cθ
where C i ¼ cos ðiÞ and Si ¼ sin ðiÞ. The quadrotor motion equations can be expressed as follows (see also in [28]): ξ€ x ¼ ax1 ξ_ x þ f ð cos ϕ cos ψ sin θ þ sin ϕ sin ψÞ=m; ξ€ y ¼ ay1 ξ_ y f ð cos ψ sin ϕ cos ϕ sin θ sin ψÞ=m; ξ€ z ¼ az1 ξ_ z þ f cos θ cos ϕ=m g; θ€ ¼ aθ1 θ_ þ cθ ðη; η_ Þθ_ þ aθ2 τθ ; ϕ€ ¼ aϕ1 ϕ_ þ cϕ ðη; η_ Þϕ_ þ aϕ2 τϕ ; ψ€ ¼ aψ1 ψ_ þ cψ ðη; η_ Þψ_ þ aψ2 τψ ;
ð1Þ
where m denotes the rotorcraft mass, g is the gravity constant, aj1 ðj ¼ x; y; z; θ; ϕ; ψÞ and ai2 ði ¼ θ; ϕ; ψÞ are the rotorcraft parameters, ci ðη; η_ Þði ¼ θ; ϕ; ψÞ are the T Coriolis terms, which contain the gyroscopic and centrifugal terms, and f and τ ¼ τθ τϕ τψ are the external bodyfixed frame force and torque acting on the aerial vehicle, respectively. f and τ can be obtained by f ¼ T 1 þ T 2 þ T 3 þ T 4 ; τθ ¼ lmc ðT 1 T 3 Þ; τϕ ¼ lmc ðT 2 T 4 Þ; τψ ¼ kf ðT 1 T 2 þ T 3 T 4 Þ; where lmc is the distance from each motor to the center of mass, kf 40 denotes the force-tomoment scaling factor, and T i ði ¼ 1; 2; 3; 4Þ are the thrusts generated by the four rotors. The thrusts T i ði ¼ 1; 2; 3; 4Þ take the following forms: T i ¼ k ω ω2i ;
i ¼ 1; 2; 3; 4;
where kω is a positive constant and ωi ði ¼ 1; 2; 3; 4Þ are the rotational velocities of the four rotors. Let ui ði ¼ z; θ; ϕ; ψÞ denote the control inputs and have the following forms: uz ¼ ω21 þ ω22 þ ω23 þ ω24 ;
uθ ¼ ω21 ω23 ;
uϕ ¼ ω22 ω24 ;
uψ ¼ ω21 ω22 þ ω23 ω24 :
The control inputs ui ði ¼ z; θ; ϕ; ψÞ, which are proportional to the force f and the torques τθ ; τϕ , and τψ , respectively, are defined to simplify the expressions of the control inputs in the following mathematical derivation. In practical applications, the real control inputs would be ωi ði ¼ 1; 2; 3; 4Þ or the control inputs of the motor control electronics rather than ui ði ¼ z; θ; ϕ; ψÞ. Furthermore, let bx ¼ g;
by ¼ g;
bz ¼ k ω =m;
bθ ¼ lmc kω aθ2 ;
bϕ ¼ lmc kω aϕ2 ;
bψ ¼ kf k ω aψ2 :
5498
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
Then from Eq. (1), the following expressions can be obtained: ξ€ x ¼ bNx ðθ þ qx Þ; ξ€ y ¼ bN ðϕ þ qy Þ; y
ξ€ z ¼ bNz ðuz þ qz Þ g; θ€ ¼ bNθ ðuθ þ qθ Þ; ϕ€ ¼ bN ðuϕ þ qϕ Þ; ϕ
ψ€ ¼ bNψ ðuψ þ qψ Þ;
ð2Þ
where the superscript N denotes the nominal parameter and qi ði ¼ x; y; z; θ; ϕ; ψÞ take the following forms: qx ¼ ð ax1 ξ_ x þ cos ϕ cos ψ sin θk ω uz =m bNx θ; þ sin ϕ sin ψkω uz =mÞ=bNx þ d x ; qy ¼ ð ay1 ξ_ y þ cos ψ sin ϕk ω uz =m bN ϕ; cos ϕ sin θ sin ψk ω uz =mÞ=bN þ dy ; y
qz ¼ ð az1 ξ_ z þ cos θ cos ϕbz uz bNz uz Þ=bNz þ dz ; qi ¼ ð ai1 _i þ ci ðη; η_ Þ_i þ bi ui bNi ui Þ=bNi þ d i ; i ¼ θ; ϕ; ψ;
y
ð3Þ
and di ði ¼ x; y; z; θ; ϕ; ψÞ are external disturbances. In Eq. (2), qi ði ¼ x; y; z; θ; ϕ; ψÞ are called equivalent disturbances. By ignoring the equivalent disturbances, model (2) indicates the nominal simplified linear model of the quadrotor. It should be noted that nonlinear dynamics and coupling are also considered as uncertainties and included in qi ði ¼ x; y; z; θ; ϕ; ψÞ. Assumption 1. The pitch and roll angles satisfy that θ A ½ π=2 þ δθ ; π=2 δθ and ϕA ½ π=2 þ δϕ ; π=2 δϕ , where δθ and δϕ are positive constants. Remark 1. The use of the propagation equations to update the Euler rotations of the helicopter body with respect to the inertial frame α is limited, because of the Euler angles singularities. In order to prevent singularities in the Euler angle, the quadrotor is required to avoid overturning during the flight [29]. Assumption 2. The sum of thrusts satisfies that f 4δf with δf being a positive constant. Remark 2. The above assumption is made to avoid f¼ 0. If the initial conditions in the vertical channel are bounded with boundaries as shown in [30] and the value of j€r z j is limited so that the descent maneuvers are slower than a free fall condition, f 40 can be guaranteed. Assumption 3. The uncertain parameters kω ; k f ; lmc ; m, aj1 ðj ¼ x; y; z; θ; ϕ; ψÞ, and ai2 ði ¼ θ; ϕ; ψÞ are bounded. The nominal parameters bNi ði ¼ x; y; z; θ; ϕ; ψÞ are positive and satisfy that jbi bNi jobNi ði ¼ z; θ; ϕ; ψÞ. Define constants ρi ði ¼ x; y; z; θ; ϕ; ψÞ as ρx ¼ maxjbNx cos ϕ cos ψ sin θf =m=θj=bNx , ρy ¼ maxj cos ψ sin ϕf =m=ϕ bNy j=bNy , ρz ¼ maxj cos θ cos ϕbz bNz j=bNz , and ρi ¼ max jbi bNi j=bNi ði ¼ θ; ϕ; ψÞ. Assumption 4. ρi ði ¼ x; yÞ satisfy that ρi o1. Remark 3. The above assumption can be guaranteed, if sufficiently large bNi ði ¼ x; yÞ are applied and Assumptions 1 and 2 hold. From Assumptions 1 through 4, it follows that 0 r ρi o1ði ¼ x; y; z; θ; ϕ; ψÞ.
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
5499
Assumption 5. The external disturbances di ði ¼ x; y; z; θ, ϕ; ψÞ are bounded. In the motion control problem discussed here, the three components of the positions and the yaw angle ξx ξy ξz ψ are chosen as outputs. This paper will investigate the robust motion controller design problem to achieve trajectory tracking of the prescribed references r i ði ¼ x; y; z; ψÞ for ξx ; ξy ; ξz , and ψ, respectively. Assumption 6. The prescribed desired references and their derivatives r ðkÞ i ði ¼ x; y; k ¼ 0; 1; 2; 3; 4Þ and rðkÞ i ði ¼ z; ψ; k ¼ 0; 1; 2Þ are piecewise uniformly bounded. Assumption 7. There exist positive constants λciη0 and λciη1 ði ¼ θ; ϕ; ψÞ such that Jci ðη; η_ ÞJ 1 r λciη0 Jη J 1 þ λciη1 J η_ J 1 . 3. Robust motion controller design The attitude and position controllers are designed in order to achieve the motion control of the quadrotor. From the rotorcraft model (2), it can be seen that the height ξz , the pitch angle θ, the roll angle ϕ, and the yaw angle ψ can be controlled directly by the inputs uz ; uθ ; uϕ , and uψ , respectively. Moreover, the positions of the pitch and roll angles result in the motions in the longitudinal and latitudinal directions. Therefore, the desired references of the pitch and roll angles are generated by the position controller based on the tracking errors of the longitudinal and latitudinal positions. The position controller also stabilizes the height. Besides, the attitude controller aims at achieving the desired reference tracking for the three attitude angles. 3.1. Position controller design The position controller is designed based on the following subsystems: ξ€ x ¼ bNx ðθ þ qx Þ; ξ€ y ¼ bN ðϕ þ qy Þ; y
ξ€ z ¼ bNz ðuz þ qz Þ g: Firstly, the virtual control inputs θ^ and ϕ^ are designed for the longitudinal and latitudinal channels with two parts: the nominal control inputs uNi ði ¼ x; yÞ and the robust compensating inputs uRC ði ¼ x; yÞ: i θ^ ¼ ux ¼ uNx þ uRC x ;
ϕ^ ¼ uy ¼ uNy þ uRC y ;
ð4Þ
where uNx ¼ ðk x1 ξx k x1 r x þ kx2 ξ_ x k x2 r_ x r€ x Þ=bNx ; uNy ¼ ðky1 ξy ky1 r y þ ky2 ξ_ y ky2 r_ y r€ y Þ=bNy ; ^ i ðsÞ; i ¼ x; y; uRC i ðsÞ ¼ F i ðsÞq ^ q^ y ¼ qy þ ϕ ϕ; ^ q^ x ¼ qx þ θ θ;
ð5Þ
and s is the Laplace operator, F i ðsÞ ¼ F if ðsÞF ig ðsÞ ði ¼ x; yÞ, and F if ðsÞ ¼ f i =ðs þ f i Þ and F ig ðsÞ ¼ gi =ðs þ gi Þ are called robust filters. kij ði ¼ x; y; j ¼ 1; 2Þ; f i and gi ði ¼ x; yÞ are positive constants to be determined. The position controller generates the desired references r θ ¼ θ^ and r ϕ ¼ ϕ^ for the pitch and roll angles to track.
5500
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
The nominal controller inputs uNx and uNy are designed for the nominal linear system (2) without considering the equivalent disturbances q^ x and q^ y , so that the closed-loop subsystem has the desired poles or modes. Moreover, the robust compensators are introduced to restrain the influence of q^ x and q^ y . In fact, the robust filters F if ðsÞ and F ig ðsÞ ði ¼ x; yÞ are constructed as the second order low pass filters and have certain frequency bandwidths, within which the gains of the robust filters are 1. If f i and gi ði ¼ x; yÞ are sufficiently large, one can expect that F if ðsÞ and F ig ðsÞ ði ¼ x; yÞ would have sufficiently wide frequency bandwidths. In this case, the robust filters can pass the low frequency primary components of the involved signals. Therefore, uRC x ^ x and q^ y , respectively and the influence of the uncertainties q^ x and uRC y would approximate q and q^ y can be reduced. The robust compensators used in the position controller are the second order filters. They are constructed without linear velocities. In fact, if the linear velocities are available and reliable, the first order filters can be used instead. However, in practical applications, linear velocities can contain significant noise, which may affect the tracking performance of the closed-loop system. Therefore, the second order robust filters are instead used here. Remark 4. It should be noted that the proposed robust controller does not rely on the assumption that the closed-loop rotational dynamics converges much faster than the closed-loop translational dynamics. The tracking errors of the pitch and roll angles are considered as uncertainties in q^ i ði ¼ x; yÞ rather than being ignored in the discussions of the translational dynamics. Similarly, the design of the control law for input uz is composed of two parts: uz ¼ uNz þ uRC z ;
uNz ¼ ðkz1 ξz k z1 r z þ kz2 ξ_ z kz2 r_ z r€ z gÞ=bNz ; uRC z ðsÞ ¼ F z ðsÞqz ðsÞ;
ð6Þ
where F z ðsÞ ¼ f z gz =ðs þ f z Þ=ðs þ gz Þ and k z1 ; kz2 ; f z , and gz are positive constants to be determined. 3.2. Attitude controller design Now consider the following subsystems: θ€ ¼ bNθ ðuθ þ qθ Þ; ϕ€ ¼ bN ðuϕ þ qϕ Þ; ϕ
ψ€ ¼ bNψ ðuψ þ qψ Þ: The attitude controller is designed to track the references r i ði ¼ θ; ϕ; ψÞ for the pitch angle θ, the roll angle ϕ, and the pitch angle ψ. Similarly, the attitude control inputs ui ði ¼ θ; ϕ; ψÞ have two parts: the nominal control input uNi and the robust compensating input uRC i : ui ¼ uNi þ uRC i ;
i ¼ θ; ϕ; ψ:
ð7Þ
The attitude nominal control laws are defined as uNi ¼ ðki1 i ki1 r þ ki2 _i ki2 r_ i r€ i Þ=bNi ;
ð8Þ
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
5501
where kij ði ¼ θ; ϕ; ψ; j ¼ 1; 2Þ are positive constants to be determined. In practical applications, r_ i ði ¼ θ; ϕÞ can be estimated by the change rates of r i ði ¼ θ; ϕÞ, respectively. Besides, r€ i ði ¼ θ; ϕÞ can be obtained in a similar way. The robust compensating inputs can be obtained as follows: uRC i ðsÞ ¼ F i ðsÞqi ðsÞ;
i ¼ θ; ϕ; ψ;
ð9Þ
where F i ðsÞ ¼ f i gi =ðs þ f i Þ=ðs þ gi Þ ði ¼ θ; ϕ; ψÞ with positive constants f i and gi to be determined. Because the uncertainties qi ði ¼ z; θ; ϕ; ψÞ and q^ j ðj ¼ x; yÞ in Eqs. (5), (6), and (9) cannot be measured directly in practical applications. Then, from (2), it follows that € N uθ : qθ ¼ θ=b θ Therefore, uRC θ can be implemented with states zθ1 and zθ2 as z_ θ1 ¼ gθ zθ1 g2θ θ þ bNθ uθ ; z_ θ2 ¼ f θ zθ2 þ ðf θ þ gθ Þθ þ zθ1 ; N uRC θ ¼ f θ gθ ðzθ2 θÞ=bθ : From the above expressions, one can see that uRC θ does not depend on the equivalent disturbance qθ . The robust compensating inputs uRC ði ¼ x; y; z; ϕ; ψÞ can be implemented in similar ways. i Remark 5. The proposed control scheme results in a linear time-invariant controller, which is easily implementable in practical applications. The block diagram of the whole robust control system is depicted in Fig. 2. 4. Robustness property analysis In this section, the proposed closed-loop system will be proven to be finite-gain L2 stable. It will analyze the robust tracking performance of the closed-loop system constructed as in the third section in two steps: the tracking errors of the height channel will be proven to converge into the given neighborhoods of the origin in a finite time; then, it will be shown that the tracking errors of other channels will also ultimately converge to the a priori set neighborhood of the origin, based on the analysis of the height channel.
(i
ri , ri , ri rj , rj , rj (j
Nominal Controller
x, y )
u
N j
u
uj
rk , rk , rk
RC j
Robust Compensator
(k
Robust Compensator
, z) Nominal Controller
Nominal Controller
, )
uiRC
ui
z
N i
u
ukN
, z,
,
uk
Quadrotor
u kRC Robust Compensator
Fig. 2. The block diagram of the robust control system.
x
, x, , ,
y
,
y
, ,
5502
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
4.1. The robustness properties of the height channel Define the error vectors of the height channel as ez ¼ ½ez1 ez2 T , where ez1 ¼ ξz r z and ez2 ¼ e_ z1 . From Eqs. (2) and (6), the following expression can be obtained: e_ z ¼ Azc ez þ Bz ðuRC z þ qz Þ; where
"
Azc ¼
# 1 ; k z2
0 kz1
ð10Þ "
Bi ¼
0
#
bNz
:
The parameters k z1 and kz2 are selected so that Azc is a Hurwitz matrix. From Eqs. (6) and (10), the following expression can be obtained: Jez J 1 r λez ð0Þ þ γ z J qz J 1 ;
ð11Þ
where γ z ¼ J ðsI 22 Azc Þ 1 Bz ð1 F z ÞJ 1 , λez ð0Þ is a positive constant satisfying that λez ð0Þ Z J eAzc t ez ð0Þ J 1 , and I 22 is a 2 2 unit matrix. It is shown in [27] that γ z can be made as small as desired, if f z and gz have sufficiently large values and f z is much larger than gz . Furthermore, there exists a positive constant λγz such that γ z oλγz =gz . Theorem 1. If Assumptions 1, 3, 5, and 6 hold, the closed-loop system of the height channel has robust tracking performance, i.e., for the given positive constant εz and the given bounded initial state ez ð0Þ, there exist positive constants t nz , f z , gz and satisfying that f z and gz are sufficiently large and f z is much larger than gz , such that the tracking error ez is bounded and satisfies jez ðtÞj r εz ; 8 t Z t nz . Proof. From (3), there exist positive constants λqez and λqcz such that Jqz J 1 r ρz Juz J 1 þ λqez Jez J 1 þ λqcz :
ð12Þ
From (6), there exist positive constants λuez and λucz such that Juz J 1 r λuez J ez J 1 þ Jqz J 1 þ λucz :
ð13Þ
Substituting (13) into (12), there exist positive constants λez and λcz such that Jqz J 1 r λez Jez J 1 þ λcz :
ð14Þ
From Eqs. (11) and (14), the following expressions can be obtained: Jqz J 1 r
λez λez ð0Þ þ λcz ; 1 λez γ z
Jez J 1 r
λez ð0Þ þ λcz γ z : 1 λez γ z
ð15Þ
If γ z is sufficiently small, it can be seen that qz and ez are bounded and thus uz is bounded; that is, there exist positive constants ηqz ; ηez , and ηuz satisfying that Jqz J 1 r ηqz ;
Jez J 1 r ηez ;
J uz J 1 r ηuz :
Furthermore, from Eqs. (6), (10), and (16), it follows that maxjezj ðtÞj r maxcTzj eAzc t ez ð0Þ þ γ z ηqz ; j
j
where czj is a 2 1 vector with one on the jth row and zero elsewhere.
ð16Þ
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
5503
From the above inequality, for the given positive constant εz and the given initial conditions, there exist a positive constant t nz and sufficiently large positive parameters f z and gz satisfying that f z is much larger than gz , such that gz 42ηqz λrz =ε and ez ðtÞr εz ; 8 t Z t nz . □ Remark 6. From the above robustness property analysis, it can be seen that the three attitude angles and the longitudinal and latitudinal positions do not affect the robust tracking properties of the height channel. Therefore, the parameters of the feedback controller in the height channel can be determined before analyzing the robustness properties of other channels.
4.2. The robustness properties of the attitude angles and the longitudinal and latitudinal channels h iT Define the error vector as e ¼ eTx eTy eTθ eTϕ eTψ , where ei ¼ ½ei1 ei2 T ði ¼ x; y; θ; ϕ; ψÞ, ex1 ¼ ξx r x ; ey1 ¼ r y ξy , ej1 ¼ j r j ðj ¼ θ; ϕ; ψÞ, and ei2 ¼ e_ i1 ði ¼ x; y; θ; ϕ; ψÞ. From Eqs. (2), (4), (5), and (7)–(9), it follows that e_ ¼ Ac e þ BU;
ð17Þ
where h iT RC RC ^ x ðuRC ^ y Þ uRC U ¼ uRC ; x þq y þq θ þ qθ uϕ þ qϕ uψ þ qψ Ac ¼ diagðAxc ; Ayc ; Aθc ; Aϕc ; Aψc Þ; " Aic ¼
0
1
ki1
k i2
#
B ¼ diagðBx ; By ; Bθ ; Bϕ ; Bψ Þ;
"
;
# 0 Bi ¼ N ; bi
i ¼ x; y; θ; ϕ; ψ:
The parameters k i1 and ki2 ði ¼ x; y; θ; ϕ; ψÞ are selected so that Aic ði ¼ x; y; θ; ϕ; ψÞ are Hurwitz matrices, respectively. Moreover, from Eqs. (5), (9), and (17), it follows that J eJ 1 r λeð0Þ þ γ J Δ J 1 ; ð18Þ T where Δ ¼ q^ x q^ y qθ =f x =gx qϕ =f y =gy qψ , γ ¼ J ðsI 1010 Ac Þ 1 BF J 1 , F ¼ diagðð1 F x Þ; ð1 F y Þ, f x gx ð1 F θ Þ; f y gx ð1 F ϕ Þ; ð1 F ψ ÞÞ, λeð0Þ is a positive constant satisfying that λeð0Þ Z J eAc t eð0Þ J 1 , and I 1010 is a 10 10 unit matrix. Lemma 1. If f x and gx are sufficiently large and f x is much larger than gx , there exist positive constants ηrieθ and ηricθ ði ¼ 0; 1; 2Þ such that the pitch reference and its first and second derivatives r θ ; r_ θ , and r€ θ satisfy that J r θ J 1 r ηr0eθ Je J 1 þ ηr0cθ ; J r_ θ J 1 r ηr1eθ gx Je J 1 þ ηr1cθ gx ; J r€ θ J 1 r ηr2eθ gx f x Je J 1 þ ηr2cθ gx f x :
ð19Þ
Proof. From Eqs. (4) and (5), the following expression can be obtained: r θ ¼ ðk x1 ex1 þ kx2 ex2 r€ x Þ=bNx þ uRC x :
ð20Þ
5504
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
From Eqs. (3), (5), (16), and (20), there exist positive constants ηr0eθ and ηr0cθ such that Jr θ J 1 r ηr0eθ J eJ 1 þ ηr0cθ :
ð21Þ
From Eqs. (2), (4), and (5), the following expression can be obtained: ^ x Þ: e_ x2 ¼ kx1 ex1 kx2 ex2 þ bNx ðuRC x þq
ð22Þ
Differentiating both sides of Eq. (20) leads to N _ RC r_ θ ¼ ðkx1 ex2 þ kx2 e_ x2 r ð3Þ x Þ=bx þ u x :
ð23Þ
It follows that there exist positive constants λr1eθ and λr1cθ such that J r_ θ J 1 r λr1eθ J eJ 1 þ J1 F x J 1 J q^ x J 1 þ gx J1 F xg J 1 J q^ x J 1 þ λr1cθ :
ð24Þ
If gx is sufficiently large, from Eqs. (3), (5), (16), (21), and (24), there exist positive constants ηr1eθ and ηr1cθ such that J r_ θ J 1 r ηr1eθ gx Je J 1 þ ηr1cθ gx :
ð25Þ
Similarly, if f x and gx are sufficiently large, there exist positive constants ηr2eθ and ηr2cθ such that J r€ θ J 1 r ηr2eθ gx f x J eJ 1 þ ηr2cθ gx f x :
ð26Þ
From Eqs. (21), (25), and (26), it can be seen that Lemma 1 holds. Theorem 2. If Assumptions 1 through 7 are met, for a priori set positive constant ε and the given bounded initial state eð0Þ, there exist a positive constant t n and positive parameters f i ; gi ði ¼ x; y; θ; ϕ; ψÞ with sufficiently large values and satisfying that f i ði ¼ x; y; θ; ϕ; ψÞ are much larger than gi ði ¼ x; y; θ; ϕ; ψÞ, respectively, and f j ; gj ðj ¼ x; yÞ are much smaller than f k ; gk ðk ¼ θ; ϕÞ, such that the tracking error e is bounded and satisfy that eðtÞ r ε; 8t Z t n . Proof. From Eqs. (3), (4), (5), (16), and (19), there exist positive constants λqex and λqcx such that J q^ x J 1 r ρx J ux J 1 þ λqex Je J 1 þ λqcx :
ð27Þ
From Eqs. (4), (5), and (27), there exist positive constants λex and λcx such that J q^ x J 1 r λex J eJ 1 þ λcx :
ð28Þ
Furthermore, if Assumption 7 is met, then from Eqs. (3), (7), (8), and (9), there exist positive constants λr1θ ; λr2θ ; λreθ ; λrcθ ; λreθ2 ; λr01θ ; λr11θ ; λr0θe , and λr1θe such that Jqθ J 1 r λr1θ J r_ θ J 1 þ λr2θ J r€ θ J 1 þ λreθ Je J 1 þ λrcθ þ λreθ2 ‖e‖21 þλr01θ J r θ J 1 J r_ θ J 1 þ λr11θ J r_ θ J 21 þ λr0θe Jr θ J 1 Je J 1 þλr1θe J r_ θ J 1 J eJ 1 :
ð29Þ
If Lemma 1 holds and f x and gx are sufficiently large satisfying that f x c gx , there exist positive constants λeθ ; λe2θ , and λcθ such that Jqθ J 1 r λeθ gx f x Je J 1 þ λe2θ g2x ‖e‖21 þ λcθ gx f x :
ð30Þ
Similarly, there exist positive constants λey ; λcy ; λeϕ ; λe2ϕ , λcϕ ; λeψ ; λe2ψ , and λcψ such that J q^ y J 1 r λey J eJ 1 þ λcy ; Jqϕ J 1 r λeϕ gy f y Je J 1 þ λe2ϕ g2y ‖e‖21 þ λcϕ gy f y ; Jqψ J 1 r λeψ Je J 1 þ λe2ψ ‖e‖21 þ λcψ :
ð31Þ
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
5505
The Tsinghua Autonomous Quadrotor System Onboard avionic electronic system
Ground station
GPS
Ultrasonic
IMU
Camera
Flight control computer
Zigbee wireless
Zigbee wireless Fig. 3. Configuration of the TAQS.
From Eqs. (28), (30), and (31), it follows that J Δ J 1 r λe J eJ 1 þ λe2 ‖e‖21 þ λc ;
ð32Þ
where λe ¼ λex þ λey þ λeθ þ λeϕ þ λeψ , λe2 ¼ λe2θ þ λe2ϕ þ λe2ψ , and λc ¼ λcx þ λcy þ λcθ þ λcϕ þ λcψ . If f i ; gi ði ¼ x; y; θ; ϕ; ψÞ have sufficiently large values and satisfying that f i ði ¼ x; y; θ; ϕ; ψÞ are much larger than gi ði ¼ x; y; θ; ϕ; ψÞ, respectively, and f j ; gj ðj ¼ x; yÞ are much smaller than f k ; gk ðk ¼ θ; ϕÞ, γ can also be made as small as desired. If γ is sufficiently small and satisfies that pffiffiffi ðλe þ λe2 J eJ 1 Þð γ þ γÞ r 1; ð33Þ then from Eqs. (18) and (32), it follows that pffiffiffi pffiffiffi J Δ J 1 r ðλeð0Þ þ γ λc þ γλc Þ= γ :
ð34Þ
It follows that J eJ 1 r λeð0Þ þ
pffiffiffi γ ηΔ ;
where ηΔ is a positive constant satisfying that ηΔ Z λeð0Þ þ From Eq. (33) follows the attractive region of e as 1 pffiffiffi 1 e : Je J 1 r λe2 =ð γ þ γÞ λe2 λe :
ð35Þ pffiffiffi γ λc þ γλc . ð36Þ
Thus, e can remain inside the attractive region, if e starts from it and the following inequality holds: pffiffiffi 1 pffiffiffi 1 λeð0Þ þ γ ηΔ r λe2 =ð γ þ γÞ λe2 λe : ð37Þ For the given initial state eð0Þ, if γ is made sufficiently small, Eq. (37) can be met. In this case,
5506
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
Eq. (33) holds, if the initial state satisfies 1 pffiffiffi 1 Jeð0ÞJ 1 oλe2 =ð γ þ γÞ λe2 λe :
ð38Þ
From Eqs. (5), (9), (17) and (35), it follows that pffiffiffi maxjej ðtÞjr maxjcTj eAc t eð0Þj þ γ ηΔ ; j
j
Table 1 Controller parameters in experimental test. Value
Parameter
Value
Parameter
Value
kx1 kθ1 kx2 kθ2 fx fθ gx gθ
0.05 6 0.5 2 2 20 0.4 5
k y1 k ϕ1 k y2 k ϕ2 fy fϕ gy gϕ
0.05 6 0.5 2 2 20 0.4 4
k z1 k ψ1 k z2 k ψ2 fz fψ gz gψ
2.5 4 45 2 3 30 0.6 8
Longitudinal position (m)
Parameter
4 Reference signal Longitudinal position
3 2 1 0 0
50
100
150
200
Latitudinal position (m)
Time (s)
4 Reference signal Latitudinal position
3 2 1 0 0
50
100
150
200
Time (s)
Height (m)
1.5 Reference signal Height 1
0.5 0
50
100
150
200
Time (s)
Fig. 4. Responses of the longitudinal, latitudinal, and height positions. (a) Longitudinal response, (b) Latitudinal response and (c) Height response.
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
5507
Yaw angle (deg)
4 Reference signal Yaw angle
2 0 −2
0
50
100
150
200
Time (s)
Pitch angle (deg)
20 Reference signal Pitch angle
10 0 −10
0
50
100
150
200
Time (s)
Roll angle (deg)
20 Reference signal Roll angle
10 0 −10
0
50
100
150
200
Time (s)
Fig. 5. Responses of the yaw, pitch, and roll angles. (a) Yaw response, (b) Pitch response and (c) Roll response.
where cj is a 10 1 vector with one on the jth row and zeros elsewhere. Then, it leads to the similar robust tracking properties with the height channel; that is, for a given positive constant ε and the given initial condition eð0Þ, there exist a positive constant t n and positive parameters f i ; gi ði ¼ x; y; θ; ϕ; ψÞ with sufficiently large values and satisfying that f i ði ¼ x; y; θ; ϕ; ψÞ are much larger than gi ði ¼ x; y; θ; ϕ; ψÞ, respectively, and f j ; gj ðj ¼ x; yÞ are much smaller than f k ; gk ðk ¼ θ; ϕÞ, such that jeðtÞj r ε; 8 t Z t n . □
5. Experimental results The Tsinghua Autonomous Quadrotor System (TAQS) has been developed using the mechanical frame of X-aircraft X650 to check the effectiveness of the proposed robust control method. The experimental platform includes an onboard avionic electronic system, which consists of an inertial measurement unit module, a GPS module, a sonar sensor, a camera, and a flight control computer. The C code is programmed based on the designed controller, then is compiled and downloaded to the flight control computer. The attitude loop can run at a rate of 100 Hz, which is also the updating rate of the IMU data collection and data fusion. The position loop runs at 33 Hz. The hardware configuration of the unmanned quadrotor system is depicted in Fig. 3. The nominal parameters of the TAQS are bNx ¼ 9:81; bNy ¼ 9:81; bNz ¼ 1; bNθ ¼ 9:7; bNϕ ¼ 9:2, and bNψ ¼ 15:99. The controller parameters are selected as shown in Table 1.
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
Latitudinal error (m)
Longitudinal error (m)
5508
1.5 Error with nominal controller Error with robust controller
1 0.5 0 −0.5
0
50
100 Time (s)
150
200
1.5 Error with nominal controller Error with robust controller
1 0.5 0 −0.5
0
50
100
150
200
Time (s)
Height error (m)
0.2 Error with nominal controller Error with robust controller
0.1 0 −0.1
0
50
100
150
200
Time (s)
Fig. 6. Comparisons of the position tracking errors. (a) Longitudinal response, (b) Latitudinal response and (c) Height response.
The indoor task that the quadrotor has to carry out is to follow references with amplitudes of 2 m, 2 m, and 0.5 m for the longitudinal, latitudinal, and vertical positions simultaneously. The yaw angle is required to stay at 01. The corresponding responses of the positions ξx , ξy , and ξz and the three attitude angles ψ, θ, and ϕ are shown in Figs. 4 and 5, respectively. The steady-state errors of the longitudinal, latitudinal, and vertical positions are 0.2 m, 0.3 m, and 0.04 m, respectively in the trajectory tracking mission. Besides, from the model (1), it can be seen that the there exist inter-axis couplings between channels. By the proposed robust control method, good tracking performance can be achieved under the influence of various uncertainties. Because the quadrotor is underactuated in the longitudinal and latitudinal directions, the pitch and roll angles are required to track the desired references obtained by the longitudinal and latitudinal position tracking errors. The responses of the pitch and roll angles are depicted in the second and third subfigures of Fig. 5. Furthermore, the proposed robust controller is compared with the nominal controller. The tracking errors of the three components of the position vector are compared in Fig. 6. From the figure, one can see that the proposed control method achieves better dynamical and steady-state performance, especially in the longitudinal and latitudinal directions. The tracking errors in the height channel are less improved, because the nonlinear model in this channel is simple compared to the models in the longitudinal and latitudinal channels as shown in Eq. (1). Therefore, the proposed control method can reduce the effects of the uncertainties including parametric perturbations, coupling, nonlinear dynamics, and external disturbances on the closed-loop control system, especially in the longitudinal and latitudinal channels, effectively.
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
5509
6. Conclusions A robust controller was proposed to achieve the motion control of a quadrotor. This control approach is a two-stage scheme including an attitude controller and a position controller. Both the attitude and position controllers consist of a nominal controller to achieve the desired tracking performance and a robust compensator to restrain the influence of uncertainties. It is proven that the position tracking errors are ultimately bounded and the boundaries can be specified. Experimental results on the quadrotor demonstrate the effectiveness of the control method. In further, the proposed robust two-loop control methods will be applied on the quadrotor experimental platform to achieve more aggressive missions under external wind disturbances. Acknowledgments This work was supported by the National Natural Science Foundation of China under Grants 61174067, 61203071, 61374054, and 61473324, as well as Shaanxi Province Natural Science Foundation Research Projection under Grants 2013JQ8038.
References [1] M.S. Mahmoud, A.B. Koesdwiady, Improved digital tracking controller design for pilot-scale unmanned helicopter, J. Frankl. Inst. 349 (2012) 42–58. [2] H. Liu, X. Wang, Y. Zhong, Robust position control of a lab helicopter under wind disturbances, IET Control Theory Appl. 8 (2014) 1555–1565. [3] H. Liu, X. Jiang, Y. Zhong, Robust hierarchical control of a laboratory helicopter, J. Frankl. Inst. 351 (2014) 259–276. [4] R. Mahony, V. Kumar, P. Corke, Multirotor aerial vehicles: modeling, estimation, and control of quadrotor, IEEE Robot. Autom. Mag. 19 (2012) 20–32. [5] G.M. Hoffmann, H. Huang, S.L. Waslander, C.J. Tomlin, Precision flight control for a multi-vehicle quadrotor helicopter testbed, Control Eng. Pract. 19 (2011) 1023–1036. [6] E. Altug, J.P. Ostrowski, R. Mahony, Control of a quadrotor helicopter using visual feedback, in: Proceedings of the IEEE International Conference on Robotics and Automation, Washington, DC, USA, 2002, pp. 72–77. [7] P. Castillo, A. Dzul, R. Lozano, Real-time stabilization and tracking of a four-rotor mini rotorcraft, IEEE Trans. Control Syst. Technol. 12 (2004) 510–516. [8] C. Aguilar-Ibanez, H. Sira-Ramirez, M.S. Suarez-Castanon, E. Martinez-Navarro, M.A. Moreno-Armendariz, The trajectory tracking problem for an unmanned four-rotor system: flatness-based approach, Int. J. Control 85 (2012) 69–77. [9] K. Tanaka, H. Ohtake, H.O. Wang, A practical design approach to stabilization of a 3-DOF RC helicopter, IEEE Trans. Control Syst. Technol. 12 (2004) 315–325. [10] R. Zhang, Q. Quan, K.Y. Cai, Attitude control of a quadrotor aircraft subject to a class of time-varying disturbances, IET Control Theory Appl. 5 (2011) 1140–1146. [11] A. Tayebi, S. McGilvray, Attitude stabilization of a VTOL quadrotor aircraft, IEEE Trans. Control Syst. Technol. 14 (2006) 562–571. [12] J.F. Guerrero-Castellanos, N. Marchand, A. Hably, S. Lesecq, J. Delamare, Bounded attitude control of rigid bodies: real-time experimentation to a quadrotor mini-helicopter, Control Eng. Pract. 19 (2011) 790–797. [13] K. Alexis, G. Nikolakopoulos, A. Tzes, Switching model predictive attitude control for a quadrotor helicopter subject to atmosphere disturbances, Control Eng. Pract. 10 (2011) 1195–1207. [14] A. Das, K. Subbarao, F. Lewis, Dynamic inversion with zero-dynamics stabilization for quadrotor control, IET Control Theory Appl. 3 (2009) 303–314. [15] Z. Zuo, Trajectory tracking control design with command-filtered compensation for a quadrotor, IET Control Theory Appl. 4 (2010) 2343–2355.
5510
H. Liu et al. / Journal of the Franklin Institute 351 (2014) 5494–5510
[16] S. Bertrand, N. Guenard, T. Hamel, H. Piet-Lahanier, L. Eck, A hierarchical controller for miniature VTOL UAVs: design and stability analysis using singular perturbation theory, Control Eng. Pract. 19 (2011) 1099–1108. [17] S. Yin, H. Luo, S.X. Ding, Real-time implementation of fault-tolerant control systems with performance optimization, IEEE Trans. Ind. Electron. 61 (2014) 2402–2411. [18] S. Yin, S.X. Ding, A. Haghani, H. Hao, P. Zhang, A comparison study of basic data-driven fault diagnosis and process monitoring methods on the benchmark Tennessee Eastman process, J. Process Control 22 (2012) 1567–1581. [19] S. Yin, S.X. Ding, A.H.A. Sari, H. Hao, Data-driven monitoring for stochastic systems and its application on batch process, Int. J. Syst. Sci. 44 (2013) 1366–2376. [20] X. Zhao, S. Yin, H. Li, B. Niu, Switching stabilization for a class of slowly switched systems, IEEE Trans. Autom. Control (2014), http://dx.doi.org/10.1109/TAC.2014.2322961, in press. [21] X. Zhao, L. Zhang, P. Shi, H.R. Karimi, Robust control of continuous-time systems with state-dependent uncertainties and its application to electronic circuits, IEEE Trans. Ind. Electron. 61 (2014) 4161–4170. [22] M.S. Mahmoud, Robust control of robot arms including motor dynamics, Int. J. Control 58 (1993) 853–873. [23] M.S. Mahmoud, H.Y. Hajeer, A globally convergent adaptive controller for robot manipulators, IEEE Trans. Autom. Control 39 (1994) 148–151. [24] M.S. Mahmoud, Y. Xia, Applied Control Systems Design: State-Space Methods, Springer-Verlag, UK, 2012. [25] Y. Zhong, Robust output tracking control of SISO plants with multiple operating points and with parametric and unstructured uncertainties, Int. J. Control 75 (2002) 219–241. [26] B. Zheng, Y. Zhong, Precision flight control for a multi-vehicle quadrotor helicopter testbed, IEEE Trans. Ind. Electron. 58 (2011) 660–670. [27] H. Liu, G. Lu, Y. Zhong, Robust LQR attitude control of a 3-DOF lab helicopter for aggressive maneuvers, IEEE Trans. Ind. Electron. 60 (2013) 4627–4636. [28] P. Castillo, A. Dzul, R. Lozano, Real-time stabilization and tracking of a four-rotor mini rotorcraft, IEEE Trans. Control Syst. Technol. 12 (2004) 510–516. [29] I.A. Raptis, K.P. Valavanis, W.A. Moreno, A novel nonlinear backstepping controller designed for helicopters using the rotation matrix, IEEE Trans. Control Syst. Technol. 19 (2011) 465–473. [30] H. Sira-Ramirez, R. Castro-Linarez, E. Liceaga-Castro, A Liouvillian systems approach for the trajectory planningbased control of helicopter models, Int. J. Robust Nonlinear Control 10 (2000) 301–320.