Attitude Stabilization with Real-time Experiments of a ... - Springer Link
J Intell Robot Syst (2012) 65:123–136 DOI 10.1007/s10846-011-9584-2
Attitude Stabilization with Real-time Experiments of a Tail-sitter Aircraft in Horizontal Flight Octavio Garcia · Pedro Castillo · K. C. Wong · Rogelio Lozano
Received: 15 February 2011 / Accepted: 18 April 2011 / Published online: 27 August 2011 © Springer Science+Business Media B.V. 2011
Abstract This paper focusses on the attitude stabilization of a mini tail-sitter aircraft, considering aerodynamic effects. The main characteristic of this vehicle is that it operates in either the hover mode for launch and recovery, or the horizontal mode during cruise. The dynamic model is obtained using the Euler–Lagrange formulation, and aerodynamic effects are obtained by studying the propeller effects. A nonlinear saturated Proportional-Integral-Derivative (SPID) control with compensation of aerodynamic moments is proposed in order to achieve the asymptotic sta-
O. Garcia (B) · R. Lozano Laboratoire Franco-Mexicain d’Informatique et Automatique, LAFMIA UMI 3175, CNRS-CINVESTAV, Mexico City, Mexico e-mail: [email protected] R. Lozano e-mail: [email protected], [email protected] P. Castillo · R. Lozano Laboratoire Heudiasyc UMR CNRS 6599, Université de Technologie de Compiègne, Compiègne, France P. Castillo e-mail: [email protected] K. C. Wong School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, Sydney 2006, NSW, Australia e-mail: [email protected]
bilization of the vehicle in horizontal mode. In addition, a homemade inertial measurement unit (HIMU) is built for operating the complete operational range of the vehicle (including vertical and horizontal modes). Finally, simulation results are presented for validating the control law, and practical results are obtained in real-time during the flight. Keywords Tail-sitter · VTOL UAV · Dynamic model · Control technique · Lyapunov function
1 Introduction The use of unmanned aerial vehicles (UAVs) for surveillance and combat operations is likely to increase rapidly in the next 10–20 years as the performance requirements for military aircraft go beyond the capabilities of human pilots and the costs of manned flight become prohibitive. UAVs help eliminate these bottlenecks by allowing the pilot to remain on the ground (or in another aircraft), controlling the UAV remotely. This approach represents an enormous challenge to control theory since many of the functions of a pilot will now have to be taken over by control systems. In addition, UAVs will need to be more maneuverable and more modular, necessitating improved tools for control law design.
124
UAVs of the future will likely undergo frequent configuration and mission changes. This will require the ability to rapidly design, redesign and implement inner-loop control laws; however, most of the current design techniques (especially for nonlinear systems) are not up to this task. Therefore, researchers need to develop automated design tools to facilitate the use of recent and new control design techniques. In order to accommodate the rapid prototyping that will be required for UAVs, these tools should take full advantage of the auto-code generation capabilities of many of the higher level programming languages that controls researchers currently use. The spirit of this effort is to develop tools that automate the entire control design process from model development to hardware-in-the-loop simulation. Many practical autopilot systems are designed based on classical design techniques. The control structure can be found in many flight control or avionics references such as [4, 16, 17]. Modern control techniques are also popular in aircraft autopilot applications [2, 13]. Fuzzy logic and adaptive control for non- minimum phase autopilot designs are discussed in [5] and [1], respectively. In [21], the authors use modern control techniques based on typical pitch attitude control loop for autopilot analysis under icing condition. Even though there are so many modern control techniques that can be adopted for autopilot application, classical designs remain a very interesting and effective approach for autopilot designs. In this paper, a nonlinear control strategy based on Saturated-Proportional-Integral and Derivative technique is presented. In order to obtain this control law, a nonlinear dynamic model with aerodynamic effects is derived. The outline of this paper is the following: In Section 2, an operational description of the mini VTOL (Vertical Take Off and Landing) aircraft is introduced in horizontal flight. In addition, the dynamic model, by employing the Euler–Lagrange formulation, with an aerodynamic analysis of the vehicle is also obtained in this section. The stability analysis, based on the Lyapunov theory, is employed to obtain the control law and the procedure is presented in Section 3. In order to validate the proposed control strategy, several simulations and real-time experiments are carried
J Intell Robot Syst (2012) 65:123–136
out, the more representatives graphs are illustrated in Sections 4 and 5. Finally, the conclusions are discussed in Section 6.
2 The C-Plane II Vehicle This section addresses the operational description and the equations of motion of the aerial vehicle (C-Plane II UAV). Aerodynamic and wind effects of the longitudinal dynamics of this aircraft are studied, for this, the nonlinear mathematical equations to represent this system are derived using the Euler–Lagrange formulation. 2.1 Operational Description The C-Plane II vehicle is a tail-sitter VTOL UAV with a propulsion system having a pair of brushless motors in coaxial configuration (see Fig. 1). This vehicle is capable of performing vertical flight as a helicopter and operating horizontal flight with the same effectiveness as a conventional airplane. Moving through the air, the vehicle utilizes the lift produced by outer body lift surfaces in order to attain the horizontal flight. For exploring the forward flight stability, the mini aircraft possesses aerodynamic control surfaces that allow rotational degrees of freedom in order to control pitch, roll and yaw motion. The aileron-elevon system controls the pitch and roll coupled motion, and the differential speed control of the two rotors also
Fig. 1 C-Plane II UAV
J Intell Robot Syst (2012) 65:123–136
125
2.2 Modeling
spectively. Thus, RIB represents the transformation matrix from the body frame to the inertial frame ⎛ ⎞ cθ cψ sφ sθ cψ − cφ sψ cφ sθ cψ + sφ sψ RIB = ⎝ cθ sψ sφ sθ sψ + cφ cψ cφ sθ sψ − sφ cψ ⎠ −sθ sφ cθ cφ cθ
The dynamic model for forward flight of this mini VTOL aircraft, considering the aerodynamic effects, is obtained by employing the Euler– Lagrange formulation. This formulation is introduced as follows.
where the shorthand notation of sa = sin(a) and ca = cos(a) is used. For this matrix, the order of the rotations is considered as yaw, pitch and roll (ψ, θ, φ) [20]. The attitude kinematics is described as
2.2.1 Preliminaries
ˆ ˙ IB = RIB R
regulates the roll motion whereas the direction of the vehicle is manipulated by the rudder. The vehicle thrust is regulated by the velocity of the propulsion system [12].
Consider an inertial fixed frame and a body frame fixed attached to the center of gravity of the aircraft denoted by I = {xI , yI , zI } and B = {xB , yB , zB }, respectively (see Fig. 2). The stability frame S = {xS , yS , zS } and the wind frame W = {xW , yW , zW } are considered during the rotation of the wind velocity vector, [22]. Assume the generalized coordinates of the mini UAV as q = (x, y, z, ψ, θ, φ)T ∈ R6 , where ξ = (x, y, z)T ∈ R3 represents the translation coordinates relative to the inertial frame, and η = (ψ, θ, φ)T ∈ R3 describes the vector of three Euler angles with rotations around z, y, x axes. These angles ψ, θ, and φ are called yaw, pitch and roll, respectively. Assume the translational velocity and the angular velocity in the body frame as ν = (u, v, w)T ∈ R3 and = ( p, q, r)T ∈ R3 , re-
ˆ is a skew-symmetric matrix such that where ˆ a = × a. The relation between the Euler angles and angular velocity is written as = X (η)η˙ where X (η) is called the Euler matrix and is represented as ⎛ ⎞ −sθ 0 1 X (η) = ⎝ cθ sφ cφ 0 ⎠ cθ cφ −sφ 0 Noting that this matrix X (η) presents a singularity when the pitch angle θ reaches ±90 deg, an inertial measurement unit was designed and built at the laboratory in order to specifically solve this singularity problem. For the aerodynamic analysis, two rotation matrices BSB and BW B are used. The rotation matrix S BB represents the transformation of a vector from the body frame B to the stability frame S , and the rotation matrix BW B describes the transformation of a vector from the body frame B to the wind frame W . These rotation matrices are written as ⎛
⎞ ⎛ ⎞ cα 0 sα cβ sβ 0 ⎝−sβ cβ 0⎠ BSB = ⎝ 0 1 0 ⎠ , BW B = −sα 0 cα 0 0 1
Fig. 2 Schematic of C-Plane II vehicle
where α is the angle of attack and β are the sideslip angle [6, 7, 20, 22].
126
J Intell Robot Syst (2012) 65:123–136
Therefore, the equations of motion of the aircraft are obtained using the Euler–Lagrange formulation ˙ ˙ d ∂ L (q, q) ∂ L (q, q) − =τ (1) dt ∂ q˙ ∂q where τ = (F, )T ∈ R6 denotes the forces and ˙ = moments acting on the body frame, L(q, q) ˙ − U (q) describes the Lagrangian equation K(q, q) ˙ which consists of the total kinetic energy K(q, q) and the potential energy U (q) of the system [8]. Moreover, K and U are defined as
4. 5. 6. 7.
˙ Sat(x) ≤ k4 ˙x
T 1 ˙ x 2 M(y) − C(q, z) x = 0 ∀x, y, z ∈ Rn ˙ ∀x, y ∈ Rn C(x, y) = 12 M(y) C(x, y)z ≤ kc yz ∀x, y, z ∈ Rn
where k1 , k2 , k3 , k4 , and kc are positive constants. In the stability analysis of this paper, the notation λmin {A} and λMax {A} are used to indicate the smallest and largest eigenvalues, respectively, of a symmetric positive definite matrix A(x), for n any x ∈ √R . The norm of vector x is defined as x = x T x and the induced norm of a matrix is also defined as A = λMax {AT A}.
K = Kt + Kr
2.2.2 Forces
U = mgz Kt =
1 ˙T ˙ ξ mξ 2 1 Kr = η˙ T M(η)η˙ 2
The forces acting on the aircraft include those of the propulsion system F p and aerodynamic effects F a . These forces are described as follows
where m ∈ R is the mass of the vehicle, Kt is the translational kinetic energy and Kr is the rotational kinetic energy with M(η) = X (η)T I X (η). I ∈ R3×3 denotes the moments of inertia of the mini VTOL aircraft. The saturation function Sat(x) is defined as follows
F = F p + Fa
⎧ x> ⎨ for Sat(x) = x for − ≤ x ≤ ⎩ − for x < −
(2)
with ⎛
⎞ Sat(x1 ) ⎜ ⎟ .. Sat(x) = ⎝ ⎠ . Sat(xn )
with ⎛
⎞ ⎛ ⎞ Tc −D B ⎝ F p = ⎝ 0 ⎠ , F a = BB Y ⎠ S BW 0 L where Tc is the thrust force of the two motors (Tc = T1 + T2 ). The lift force L, sideforce Y and drag force D are defined as aerodynamic forces, and g is the acceleration due to gravity, (see Fig. 2) [12]. In this analysis, the thrust force is oriented parallel to the axis xB of the body frame. 2.2.3 Moments
where x ∈ Rn and > 0. Now, some properties about the saturation function Sat(x), inertia matrix M(η), and coriolis ˙ need to be and centrifugal forces matrix C(η, η) introduced [9, 11]:
The moments generated on the mini VTOL aircraft are due to actuators (actuator moment act , reaction moment rot and gyroscopic moment gyro ), and the aerodynamic effects a . These moments are defined as follows
1. Sat(x) ≤ k1 x 2. Sat(x) ≤ k2 3. Sat(x)2 ≤ k3 Sat(x)T x
⎞
L = ⎝ M ⎠ = act + rot + gyro + a
N ⎛
J Intell Robot Syst (2012) 65:123–136
127
where
with ⎛
⎞
⎛
⎞
τφ Irot1 ω˙ r1 − Irot2 ω˙ r2 ⎠, act = ⎝ τθ ⎠ , rot = ⎝ 0 τψ 0 ⎛
gyro
⎞
1 T η˙ M(η)η˙ 2
(7)
From Eqs. 6 and 7
⎞ ¯ L 0 ¯ ⎠ = ⎝ r(Ir1 ωr1 − Ir2 ωr2 ) ⎠ , a = ⎝ M ¯ q(−Ir1 ωr1 + Ir2 ωr2 ) N
⎛
˙ (η) η˙ − M (η) η¨ + M
where τφ = a ( fe2 − fe1 ), τθ = e ( fe1 + fe2 ) and τψ = e fr are the control inputs with a and e that represent the distance from the center of mass to the forces fe1 and fe2 . ωri denotes the angular velocity of the rotor, Iri is the inertia moment of the propeller and Iroti is the moment of inertia of ¯ M ¯ and the rotor around its axis for i = 1, 2. L, N¯ are aerodynamic rolling, pitching and yawing moments, respectively (see Fig. 2) [12].
1 ∂ T η˙ M (η) η˙ = 2 ∂η
(8)
From Eq. 8, the Coriolis and Centrifugal vector is defined as ˙ η˙ = C(η, η)
1 ∂ T ˙ η˙ M(η) η˙ M(η) − 2 ∂η
(9)
Thus, the dynamic model for the rotational motion is rewritten as ˙ η˙ = M(η)¨η + C(η, η)
2.2.4 Translational and Rotational Dynamics Since the Lagrangian equation contains no cross˙ terms in the kinetic energy combining ξ˙ with η, the Euler–Lagrange Eq. 1 can be partitioned into dynamics for ξ coordinates and η coordinates [3]. The translational motion can be obtained using the following expression d ∂ Lt ∂ Lt =F (3) − dt ∂ ξ˙ ∂ξ with 1 T Lt = ξ˙ mξ˙ − mgz 2
Lr =
(10)
˙ are represented as where M(η) and C(η, η)
M(η) ⎛ =⎝
Ixx s2θ + I yy c2θ s2φ + Izz c2θ c2φ cθ cφ sφ (I yy − Izz ) −Ixx sθ cθ cφ sφ (I yy − Izz )
I yy c2θ + Izz s2φ
0
−Ixx sθ
0
Ixx
⎞ ⎠
(11) ⎛
(4)
Thus, after some computations, the translation motion of this vehicle is described as ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ mx¨ Tc −D B ⎝ ⎝ m y¨ ⎠ = RIB ⎝ 0 ⎠ + RIB BB Y ⎠ S BW m¨z 0 L ⎛ ⎞ 0 +⎝ 0 ⎠ (5) −mg Similarly, the rotational motion is described as d ∂ Lr ∂ Lr − = (6) dt ∂ η˙ ∂η
⎞ c11 c12 c13 ˙ = ⎝ c21 c22 c23 ⎠ C(η, η) c31 c32 c33
(12)
with ˙ θ cθ + I yy −θs ˙ θ cθ s2φ + φc ˙ 2θ sφ cφ c11 = Ixx θs ˙ 2θ sφ cφ −Izz θ˙ sθ cθ c2φ + φc ˙ θ cθ c12 = Ixx ψs ˙ θ s2φ − φc ˙ θ c2φ + ψs ˙ θ cθ s2φ −I yy θ˙ sθ sφ cφ + φc ˙ θ c2φ − ψs ˙ θ cθ c2φ + θs ˙ θ sφ cφ ˙ θ s2φ − φc +Izz φc ˙ 2θ sφ cφ − Izz ψc ˙ 2θ sφ cφ c13 = −Ixx θ˙ cθ + I yy ψc
128
J Intell Robot Syst (2012) 65:123–136
˙ θ cθ + I yy ψs ˙ θ cθ s2φ + Izz ψs ˙ θ cθ c2φ c21 = −Ixx ψs ˙ φ cφ + Izz φs ˙ φ cφ c22 = −I yy φs ˙ φ cφ + ψc ˙ θ c2φ − ψc ˙ θ s2φ ˙ θ + I yy −θs c23 = Ixx ψc ˙ θ s2φ − ψc ˙ θ c2φ + θs ˙ φ cφ +Izz ψc ˙ 2θ sφ cφ + Izz ψc ˙ 2θ sφ cφ c31 = −I yy ψc ˙ θ + I yy θs ˙ φ cφ + ψc ˙ θ s2φ − ψc ˙ θ c2φ c32 = −Ixx ψc ˙ θ s2φ − ψc ˙ θ c2φ + θs ˙ φ cφ −Izz ψc c33 = 0
(13)
2.3 Aerodynamic Analysis This section presents the aerodynamic analysis of the mini VTOL aircraft in horizontal flight. For simplicity, a two-dimensional (2-D) scheme of C-Plane II vehicle is considered to study the aerodynamic effects in presence of wind. Since the stability of this mini VTOL aircraft is affected by the propulsion system, the propeller effect (Propeller Momentum Theory) is analyzed to obtain the aerodynamic behavior of the vehicle in horizontal flight [14, 15, 19]. In Fig. 3, αs is the slipstream angle of attack, α denotes the wing angle of attack, α p describes
the propeller angle of attack, θ is pitch angle and δ expresses the angle of the ailerons. vo is the freestream velocity, vi represents the induced velocity and is directed opposite to the thrust, vs is the slipstream velocity, vds denotes the far downstream velocity, and finally, vr is the resultant velocity in the propeller slipstream. This Figure shows the wing-propeller combination where the wing is submerged in the air slipstream, and the C-Plane II vehicle is turned in the direction of the air freestream. From Fig. 3, the following aerodynamic equations can be established
vr =
2 2 2vi + vo cos α p + vo sin α p
α = tan−1
vo sin α p 2vi + vo cos α p
αs = sin−1
(14)
(15)
2vi sin α p vr
(16)
Taking, the thrust equation Tc and the downstream velocity vds Tc = ρ Avs vds − vo cos α p
(17)
and vds = vo cos α p + 2vi
(18)
where A represents the area of the rotor disc, ρ denotes the air density. Now, considering the case where air freestream velocity vo is different from zero, then the resultant velocity vr is equal to slipstream velocity vs far behind of the propeller, Eqs. 14, 17 and 18 are rewritten in a quartic form for the induced velocity as
Fig. 3 Longitudinal motion of the C-Plane II UAV
4vi4 + 4vo vi3 cos α p + vo2 vi2 =
Tc 2ρ A
2 (19)
J Intell Robot Syst (2012) 65:123–136
and vh =
129
namic moments. Since the roll motion is mechanically stable using contra-rotating propellers, the reaction and gyroscopic moments ( rot and gyro ) will essentially be zero. The nonlinear control law is described as t ˙˜ ˜ − Ki ˜ γ Sat (η(s)) + η(s) ds act = −K p Sat (η)
Tc 2ρ A
consequently Eq. 19 yields 4 3 vi vi vo 4 +4 cos α p vh vh vh 2 2 vo vi + =1 vh vh
0
(20)
Finally, the aerodynamic forces and moments are written by D = Y = L = ¯ = L ¯ = M N¯ =
1 2 ρv SC D 2 r 1 2 ρv SCY 2 r 1 2 ρv SC L 2 r 1 2 ρv Sb Cl 2 r 1 2 ¯ m ρv ScC 2 r 1 2 ρv Sb Cn 2 r
− K d Sat η˙˜ + M (η) η¨ d ˙ η˙ d − a + C(η, η)
(22)
where K p , K i , K d , are the symmetric definite positive matrices, and are proportional, integral and derivative gains of the controller, and γ > 0. The control for trajectory tracking can be given as the following. Let ηd and η˙ d denote the desired trajectories and define the error functions η˜ 1 = η − ηd η˜ 2 = η˙ − η˙ d
(21)
where S represents the wing area, c¯ is the wing chord, and b is the wing span. C D , CY and C L are aerodynamical non-dimensional coefficients of drag, sideforce and lift. Cl , Cm and Cn are aerodynamical non-dimensional coefficients of rolling, pitching and yawing moments.
3 Control Technique This section presents the control strategy for the attitude stabilization of the C-Plane II UAV in horizontal flight. A nonlinear Saturated Proportional-Integral-Derivative (SPID) control with compensation of aerodynamic moments is proposed in order to stabilize the vehicle during the flight. The stability analysis is based on the Lyapunov theory. The controller is capable of preventing aggressive spiral and phugoid motion during cruise, due to the compensation of aerody-
(23)
In this stability analysis, the variables η = η1 , η˙ = η2 , η˜ = η˜ 1 and η˙˜ = η˜ 2 are used. The integral action of the control law in Eq. 22 adds a state variable η0 in the system [10]. Using Eqs. 10 and 22, the closed-loop system is writ ten as η˜˙ 0 = γ Sat η˜ 1 + η˜ 2 η˙˜ 1 = η˜ 2
η˙˜ 2 = M(η1 )−1 −K p Sat η˜ 1 − K i η˜ 0 −K d Sat η˜ 2 − C(η1 , η2 )η˜ 2
(24)
In order to achieve the asymptotic stability of the system in Eq. 24, a Lyapunov candidate function is proposed as 1 1 V t, η˜ 0 , η˜ 1 , η˜ 2 = η˜ 2T M (η1 ) η˜ 2 + η˜ 1T K p η˜ 1 2 2 T 1 T + η˜ 0 K i η˜ 0 +γ Sat η˜ 1 M (η1 ) η˜ 2 2 (25) Notice that, 2 1 1 T η˜ 2 Mη˜ 2 ≥ λmin {M} η˜ 2 2 2
(26)
130
J Intell Robot Syst (2012) 65:123–136
2 1 T 1 η˜ 1 K p η˜ 1 ≥ λmin K p η˜ 1 2 2 2 1 T 1 η˜ 0 K i η˜ 0 ≥ λmin K i η˜ 0 2 2
(27)
(28)
Computing the time derivative of Eq. 25 and employing properties 5 and 6, it gets V˙ t, η˜ 0 , η˜ 1 , η˜ 2 = −η˜ 2T K d Sat η˜ 2 ˙ η˜ 1 T M (η1 ) η˜ 2 + γ Sat
and from property 1, the last term of Eq. 25 yields
T − γ Sat η˜ 1 K p Sat η˜ 1
T γ Sat η˜ 1 Mη˜ 2 ≤ γ k1 λMax {M} η˜ 1 η˜ 2 T −γ Sat η˜ 1 Mη˜ 2 ≥ −γ k1 λMax {M} η˜ 1 η˜ 2
T − γ Sat η˜ 1 K d Sat η˜ 2 (29)
T + γ Sat η˜ 1 C (η1 , η2 ) η˜ 2
Taking Eqs. 26, 27, 28, and 29, Eq. 25 yields
In order to demonstrate that V˙ < 0, each term of Eq. 32 is analyzed. Using properties 1 and 3, it can be written as
2 1 V t, η˜ 0 , η˜ 1 , η˜ 2 > λmin {M} η˜ 2 2 2 1 + λmin K p η˜ 1 2 2 1 + λmin {K i } η˜ 0 2 −γ k1 λMax {M} η˜ 1 η˜ 2
η˜ 2T K d Sat η˜ 2
(30) The Eq. 30 can be written as ⎡ ˜ ⎤T ⎡ η˜ 0 ⎤ 1 η0 ⎦ Q ⎣ η˜ 1 ⎦ η˜ 1 V t, η˜ 0 , η˜ 1 , η˜ 2 > ⎣ 2 η˜ 2 η˜ 2 where ⎡ ⎢ ⎢ Q=⎢ ⎣
(32)
≥
2 1 λmin {K d } Sat η˜ 2 k3
2 1 −η˜ 2T K d Sat η˜ 2 ≤ − λmin {K d } Sat η˜ 2 k3 ≤−
(31)
2 1 2 k1 λmin {K d } η˜ 2 k3
Considering property 4 for the second term of Eq. 32, it results ˙ η˜ 1 T Mη˜ 2 ≤ γ k4 λMax {M} η˜ 2 2 γ Sat
⎤
λmin {K i } 0 0 0 λmin K p −γ k1 λMax {M} ⎥ ⎥ 0 −γ k1 λMax {M} λmin {M} ⎥ ⎦
Now, applying the theorem of Sylvester [18], it implies λmin {K i } > 0 λmin {K i } λmin K p > 0 γ 2 k21 λMax {M}2 λmin K p > λmin {M} Therefore, V > 0 is a positive definite function.
(33)
(34)
Employing the property 1 in the third and fourth terms of Eq. 32, it leads T γ Sat η˜ 1 K p Sat η˜ 1
2 ≥ γ k21 λmin K p η˜ 1
T 2 −γ Sat η˜ 1 K p Sat η˜ 1 ≤ −γ k21 λmin K p η˜ 1 (35) T γ Sat η˜ 1 K d Sat η˜ 2 ≥ γ k21 λmin {K d } η˜ 1 η˜ 2 T −γ Sat η˜ 1 K d Sat η˜ 2 ≤ −γ k21 λmin {K d } η˜ 1 η˜ 2 (36)
J Intell Robot Syst (2012) 65:123–136
131
Finally, using properties 1, 2 and 7 for the last term of Eq. 32, it yields T γ Sat η˜ 1 Cη˜ 2 ≤ γ kc η˜ 2 η2 Sat η˜ 1 ≤ γ kc η˜ 2 η˜ 2 + ηd2 Sat η˜ 1 2 ≤ γ k2 kc η˜ 2 +γ kc ηd2 η˜ 2 Sat η˜ 1 2 ≤ γ k2 kc η˜ 2 +γ k1 kc ηd2 η˜ 2 η˜ 1
(37)
Then, employing Eqs. 33, 34, 35, 36 and 37, the Eq. 32 is rewritten as 2 1 V˙ t, η˜ 0 , η˜ 1 , η˜ 2 ≤ − k21 λmin {K d } η˜ 2 k3 2 +γ k4 λMax {M} η˜ 2 2 −γ k21 λmin K p η˜ 1 −γ k21 λmin {K d } η˜ 1 η˜ 2 2 +γ k2 kc η˜ 2 +γ k1 kc ηd2 η˜ 2 η˜ 1 (38) The Eq. 38 can be written as $ %T $ % η˜ 1 η˜ 1 V˙ t, η˜ 0 , η˜ 1 , η˜ 2 ≤ − P η˜ 2 η˜ 2
(39)
where ⎡ ⎢ P=⎢ ⎣
γ k21 λmin K p γ k21 2
λmin {K d } −
γ k1 kc ηd2 2
γ k21 γ k1 kc ηd2 λmin {K d } − 2 2 1 2 k λmin {K d } − γ k4 λMax {M} − γ k2 kc k3 1
Applying the theorem of Sylvester [18], it leads
λmin {K p }
$
> γ k21
%2 γ k1 kc λmin {K d } − ηd2 2 2
E = η˜ 1 ∈ R3 , η˜ 2 ∈ R3 | V˙ = 0
γ k21
1 2 k λmin {K d } − γ k4 λMax {M} − γ k2 kc k3 1
⎥ ⎥ ⎦
Then, V˙ ≤ 0 is a negative semidefinite function. Now, using the theorem of LaSalle in Eq. 39 [10], it implies
λmin {K p } > 0
$
⎤
%
(40)
Therefore, the closed-loop system is asymptotically stable about the origin. Figure 4 illustrates the block diagram of the closed-loop system of this vehicle.
4 Simulation Results -
SPID CONTROLLER
+ +
Fig. 4 Closed-loop system
UAV
In this section, the results obtained by simulation are presented in order to observe the performance in horizontal flight of the proposed control law.
132
J Intell Robot Syst (2012) 65:123–136 40
ñ
1
20
ñ0
10 0 −10 −20 5
10
15
0 −50 dñ2/dt
−100
dñ1/dt dñ /dt 0
−150 0
5
10
time [s]
time [s]
(a)
(b)
0.5
12
0
10 L/D ratio
Control law [Nm]
−30 0
Fig. 6 Control law and L/D ratio
50
ñ2
30
velocity errors [deg/s]
Position errors [deg]
Fig. 5 Angular position and velocity of the vehicle
−0.5 −1
8 6 4
−1.5
2
−2 0
0 0
time [s]
15 α [deg]
(a)
(b)
5
10
15
5
10
20
Fig. 8 Microcontroller
1 Aerodynamic coefficients
CD CL
0.8 0.6 0.4 0.2 0
0
5
10
15 α [deg]
15
20
25
30
Fig. 7 Aerodynamic coefficients C D and C L
Table 1 Parameters of the mini VTOL aircraft Parameter
Value
Airfoil shape ¯ Wing chord (c) Wing span (b ) Aspect ratio (AR) Thickness Wing area (S) Vehicle mass (m) Ixx I yy Izz
Flat plate 0.5 m 0.80 m 1.6 0.006 m 0.4 m2 0.8 kg 1.0775 kg/m2 0.4374 kg/m2 0.06399 kg/m2
Fig. 9 The homemade inertial measurement unit
25
30
J Intell Robot Syst (2012) 65:123–136
133
High pass filter
VTOL aircraft was published in [12]). This section also shows practical results during cruise. Experiments of this vehicle were realized in autonomous and semi-autonomous horizontal flight. In the semi-autonomous case, the pilot only provided the throttle control.
Rate gyro
+
+ Accelerometer
−
output Angle estimation
accel g
5.1 Experimental platform
Low pass filter
Fig. 10 Block diagram of the analog complementary filter
The airframe of C-Plane II vehicle is built from a polystyrene foam sheet (Depron) and carbon fiber-tubes. The vehicle is powered by two brushless motors driving contra-rotating propellers 10 × 4.5 in (Coaxial configuration). Table 1 shows some parameters of this vehicle. The embedded system used in this VTOL vehicle consists of a microcontroller RCM3400 having a Rabbit microprocessor which receives the data from sensors, computes the control law and then sends the control signals to actuators in order to obtain a suitable control (see Fig. 8). The Homemade Inertial Measurement Unit (HIMU) is built using three rate gyros (ADXRS150 Analog Devices), three accelerometers (ADXL203 Analog Devices), and two magnetometers (KMZ52 Philips), see Fig. 9. These sensors are placed in orthogonal positions in order to measure the aircraft motion. The main characteristic of this HIMU is that it can operate in both hover and horizontal flight solving the singularity problem. In addition, this HIMU incorporates an analog complementary filter and signal conditioning circuits in order to estimate the angle measurements (see Fig. 10). These analog filter and signal conditioning circuits are designed using active components (precision
Figure 5 shows the angular position and velocity errors during cruise. Observe that, the error signals converge to zero and the closed-loop system is asymptotically stable. The control law that stabilize the mini VTOL aircraft during the horizontal flight is depicted in the Fig. 6a. Figure 6b also illustrates the lift-to-drag ratio for some values of α. The aerodynamic coefficients C D and C L are shown in the Fig. 7. The parameters involved in the control of this vehicle are taken as k1 = 1, k2 = 1, k3 = 1, k4 = 1, kc = 1, γ = 1, and the matrices are defined as K p = diag{20, 20, 20}, K d = diag{12, 12, 12} and K i = diag{15, 15, 15}.
5 Experimental Results This section describes the embedded system and the mini VTOL aircraft (C-Plane II) developed at the Laboratory LAFMIA UMI 3175 CNRS Cinvestav Mexico (the first version of this mini
Fig. 11 Angle estimation in low and high frequencies
50
15 Angle estimation Accelerometer
10 Degrees
Degrees
Angle estimation Accelerometer
0
5 0 −5 −10
−50
0
5
10 time [s]
(a)
15
−15 0
5
10 time [s]
(b)
15
20
134
J Intell Robot Syst (2012) 65:123–136
Fig. 12 Block diagram of the embedded system
Power Interface
RC RECEIVER UAV
SPID CONTROLLER
XBEE MODEM
Microcontroller
Homemade Inertial Measurement Unit PITOT TUBE
ACCELEROMETERS
ANALOG
SIGNAL
COMPLEMENTARY
CONDITIONING
FILTER
CIRCUIT
M
Actuators
+
RATE GYROS
MAGNETOMETERS BRUSHLESS MOTORS
Fig. 13 Roll angle and control
10 Roll control [PWM]
Roll angle [deg]
10 5 0 −5
−10 0
20
40
60
80 100 120 140 160 180 time [s]
5 0 −5
−10 0
20
40
60
(a)
5 Pitch control [PWM]
5 0 −5 −10 0
20
40
60
0
−5 0
80 100 120 140 160 180 time [s]
20
40
60
(a)
Fig. 15 Yaw angle and control
6 Yaw control [PWM]
Yaw angle [deg]
80 100 120 140 160 180 time [s]
(b)
10 5 0 −5
−10
80 100 120 140 160 180 time [s]
(b)
10 Pitch angle [deg]
Fig. 14 Pitch angle and control
RC SERVOMOTORS
0
20
40
60
80 100 120 140 160 180 time [s]
(a)
4 2 0 2 −4 −6 0
20
40
60
80 100 120 140 160 180 time [s]
(b)
J Intell Robot Syst (2012) 65:123–136
operational amplifiers and integrated circuits) and passive components (resistors, capacitors, and inductors). Estimation responses in real-time of the analog filter complementary are illustrated in the Fig. 11. This responses are obtained in low and high frequencies [23]. On the other hand, one magnetometer measures the magnetic field in hover flight, while the other operates in the horizontal flight. For computing the airspeed of the mini VTOL aircraft in horizontal flight, a small pitot-static tube is added to the embedded system. This pitotstatic tube measures the difference between static and total pressures, allowing the freestream airspeed to be determined. Figure 12 shows the embedded system developed for the C-Plane II vehicle. Figures 13, 14 and 15 show the real-time responses of the vehicle orientation. Figure 13 demonstrates the effective behavior in roll motion. Pitch and yaw motion are illustrated in Figs. 14 and 15 respectively.
6 Conclusions and Future Work 6.1 Conclusions A C-Plane II aircraft has been presented in order to perform autonomous horizontal flight. The dynamic model of the vehicle has been obtained using the Euler–Lagrange formulation, and the aerodynamic terms including the wind are described by considering the propeller analysis. A nonlinear saturated Proportional-IntegralDerivative (SPID) controller with compensation of aerodynamic moments has been proposed for stabilizing the vehicle in horizontal flight. The stability analysis based on Lyapunov theory has demonstrated that the closed-loop system is asymptotically stable in the origin. In addition, a homemade inertial measurement unit (HIMU) has been developed in order to operate in both vertical and horizontal flight avoiding the singularity problem. The sensor fusion of this HIMU is obtained by using an analog complementary filter. Finally, simulations and experimental results of the closed-loop system have shown an effective behavior of the proposed control law.
135
6.2 Future work Future work for C-Plane II UAV includes incorporating additional systems like GPS and visionbased sensors to navigate autonomously. The aim of this additional instrumentation is to achieve a complete navigation flight.
References 1. Barkana, I.: Classical and simple adaptive control for nonminimum phase autopilot design. AIAA J. Guid. Control Dyn. 28(4), 631–638 (2005) 2. Blight, J.D., Gangsaas, D., Richardson, T.M.: Control law synthesis for an airplane with relaxed static stability. AIAA J. Guid. Control Dyn. 9(5), 546–554 (1986) 3. Castillo, P., Lozano, R., Dzul, A.: Modelling and Control of Miniflying Machines. Springer, London (2005) 4. Collinson, R.P.G.: Introduction to Avionics. Chapman & Hall, London (1996) 5. Cohen, K., Bossert, D.E.: Fuzzy Logic Non-minimum Phase Autopilot Design. AIAA Paper, pp. 2003–5550, Seoul (2003) 6. Cook, M.V.: Flight Dynamics Principles, 2nd edn. Butterworth-Heinemann, Oxford (2007) 7. Etkin, B., Reid, L.D.: Dynamics of Flight: Stability and Control, 3rd edn. Wiley, New York (1996) 8. Goldstein, H., Poole, C.P., Safko, J.L.: Classical Mechanics, 2nd edn. Addison Wesley Series in Physics, Addison-Wesley, Reading (1983) 9. Kelly, R.: Global Positioning of robot manipulators via PD control plus a class of nonlinear integral actions. IEEE Trans. Automat. Contr. 43(7) 934–938 (1998). doi:10.1109/9.701091 10. Khalil, H.: Nonlinear Systems, 3rd edn. Prentice Hall, Englewood Cliffs (2002) 11. Koditschek, D.E.: Natural motion for robot arms. In: IEEE Conference on Decision and Control (CDC 1984), Las Vegas, NV (1984) 12. Lozano, R.: Unmanned Aerial Vehicles Embedded Control. John Wiley-ISTE Ltd (2010) 13. Magni, J.F., Bennani, S., Terlouw, J.: Robust Flight Control: a Design Challenge. Lecture Notes in Control and Information Sciences, number 224. Springer, New York (1997) 14. McCormick, B.W.: Aerodynamics of V/STOL Flight. Dover, New York (1967) 15. McCormick, B.W.: Aerodynamics, Aeronautics and Flight Mechanics. Wiley, New York (1995) 16. McLean, D.: Automatic Flight Control Systems. Prentice Hall, Englewood Cliffs (1990) 17. Nelson, R.C.: Flight Stability and Automatic Control, 2nd edn. McGraw-Hill, New York (1998) 18. Ogata, K.: Modern Control Engineering, 5th edn. Prentice Hall, Englewood Cliffs (2009) 19. Phillips, W.F.: Mechanics of Flight. John Wiley and Sons, USA (2004)
136 20. Stengel, R.F.: Flight Dynamics. Princeton University Press, Princeton (2004) 21. Sharma, V., Voulgaris, P.G., Frazzoli, E.: Aircraft Autopilot Analysis and Envelope Protection for Operation Under Icing Conditions. AIAA J. Guid. Control Dyn. 27(3), 454–465 (2004)
J Intell Robot Syst (2012) 65:123–136 22. Stevens, B.L., Lewis, F.L.: Aircraft Control and Simulation, 2nd edn. Wiley-Interscience, New York (2003) 23. Zimmermann, M., Sulzer, W.: High bandwidth orientation measurement and control based on complementary filtering. In: Proceeding of SYROCO’91, IFAC Symposium on Robot Control, Vienna, Austria (1991)
Attitude Stabilization with Real-time Experiments of a Tail-sitter Aircraft in Horizontal Flight Octavio Garcia · Pedro Castillo · K. C. Wong · Rogelio Lozano
Received: 15 February 2011 / Accepted: 18 April 2011 / Published online: 27 August 2011 © Springer Science+Business Media B.V. 2011
Abstract This paper focusses on the attitude stabilization of a mini tail-sitter aircraft, considering aerodynamic effects. The main characteristic of this vehicle is that it operates in either the hover mode for launch and recovery, or the horizontal mode during cruise. The dynamic model is obtained using the Euler–Lagrange formulation, and aerodynamic effects are obtained by studying the propeller effects. A nonlinear saturated Proportional-Integral-Derivative (SPID) control with compensation of aerodynamic moments is proposed in order to achieve the asymptotic sta-
O. Garcia (B) · R. Lozano Laboratoire Franco-Mexicain d’Informatique et Automatique, LAFMIA UMI 3175, CNRS-CINVESTAV, Mexico City, Mexico e-mail: [email protected] R. Lozano e-mail: [email protected], [email protected] P. Castillo · R. Lozano Laboratoire Heudiasyc UMR CNRS 6599, Université de Technologie de Compiègne, Compiègne, France P. Castillo e-mail: [email protected] K. C. Wong School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, Sydney 2006, NSW, Australia e-mail: [email protected]
bilization of the vehicle in horizontal mode. In addition, a homemade inertial measurement unit (HIMU) is built for operating the complete operational range of the vehicle (including vertical and horizontal modes). Finally, simulation results are presented for validating the control law, and practical results are obtained in real-time during the flight. Keywords Tail-sitter · VTOL UAV · Dynamic model · Control technique · Lyapunov function
1 Introduction The use of unmanned aerial vehicles (UAVs) for surveillance and combat operations is likely to increase rapidly in the next 10–20 years as the performance requirements for military aircraft go beyond the capabilities of human pilots and the costs of manned flight become prohibitive. UAVs help eliminate these bottlenecks by allowing the pilot to remain on the ground (or in another aircraft), controlling the UAV remotely. This approach represents an enormous challenge to control theory since many of the functions of a pilot will now have to be taken over by control systems. In addition, UAVs will need to be more maneuverable and more modular, necessitating improved tools for control law design.
124
UAVs of the future will likely undergo frequent configuration and mission changes. This will require the ability to rapidly design, redesign and implement inner-loop control laws; however, most of the current design techniques (especially for nonlinear systems) are not up to this task. Therefore, researchers need to develop automated design tools to facilitate the use of recent and new control design techniques. In order to accommodate the rapid prototyping that will be required for UAVs, these tools should take full advantage of the auto-code generation capabilities of many of the higher level programming languages that controls researchers currently use. The spirit of this effort is to develop tools that automate the entire control design process from model development to hardware-in-the-loop simulation. Many practical autopilot systems are designed based on classical design techniques. The control structure can be found in many flight control or avionics references such as [4, 16, 17]. Modern control techniques are also popular in aircraft autopilot applications [2, 13]. Fuzzy logic and adaptive control for non- minimum phase autopilot designs are discussed in [5] and [1], respectively. In [21], the authors use modern control techniques based on typical pitch attitude control loop for autopilot analysis under icing condition. Even though there are so many modern control techniques that can be adopted for autopilot application, classical designs remain a very interesting and effective approach for autopilot designs. In this paper, a nonlinear control strategy based on Saturated-Proportional-Integral and Derivative technique is presented. In order to obtain this control law, a nonlinear dynamic model with aerodynamic effects is derived. The outline of this paper is the following: In Section 2, an operational description of the mini VTOL (Vertical Take Off and Landing) aircraft is introduced in horizontal flight. In addition, the dynamic model, by employing the Euler–Lagrange formulation, with an aerodynamic analysis of the vehicle is also obtained in this section. The stability analysis, based on the Lyapunov theory, is employed to obtain the control law and the procedure is presented in Section 3. In order to validate the proposed control strategy, several simulations and real-time experiments are carried
J Intell Robot Syst (2012) 65:123–136
out, the more representatives graphs are illustrated in Sections 4 and 5. Finally, the conclusions are discussed in Section 6.
2 The C-Plane II Vehicle This section addresses the operational description and the equations of motion of the aerial vehicle (C-Plane II UAV). Aerodynamic and wind effects of the longitudinal dynamics of this aircraft are studied, for this, the nonlinear mathematical equations to represent this system are derived using the Euler–Lagrange formulation. 2.1 Operational Description The C-Plane II vehicle is a tail-sitter VTOL UAV with a propulsion system having a pair of brushless motors in coaxial configuration (see Fig. 1). This vehicle is capable of performing vertical flight as a helicopter and operating horizontal flight with the same effectiveness as a conventional airplane. Moving through the air, the vehicle utilizes the lift produced by outer body lift surfaces in order to attain the horizontal flight. For exploring the forward flight stability, the mini aircraft possesses aerodynamic control surfaces that allow rotational degrees of freedom in order to control pitch, roll and yaw motion. The aileron-elevon system controls the pitch and roll coupled motion, and the differential speed control of the two rotors also
Fig. 1 C-Plane II UAV
J Intell Robot Syst (2012) 65:123–136
125
2.2 Modeling
spectively. Thus, RIB represents the transformation matrix from the body frame to the inertial frame ⎛ ⎞ cθ cψ sφ sθ cψ − cφ sψ cφ sθ cψ + sφ sψ RIB = ⎝ cθ sψ sφ sθ sψ + cφ cψ cφ sθ sψ − sφ cψ ⎠ −sθ sφ cθ cφ cθ
The dynamic model for forward flight of this mini VTOL aircraft, considering the aerodynamic effects, is obtained by employing the Euler– Lagrange formulation. This formulation is introduced as follows.
where the shorthand notation of sa = sin(a) and ca = cos(a) is used. For this matrix, the order of the rotations is considered as yaw, pitch and roll (ψ, θ, φ) [20]. The attitude kinematics is described as
2.2.1 Preliminaries
ˆ ˙ IB = RIB R
regulates the roll motion whereas the direction of the vehicle is manipulated by the rudder. The vehicle thrust is regulated by the velocity of the propulsion system [12].
Consider an inertial fixed frame and a body frame fixed attached to the center of gravity of the aircraft denoted by I = {xI , yI , zI } and B = {xB , yB , zB }, respectively (see Fig. 2). The stability frame S = {xS , yS , zS } and the wind frame W = {xW , yW , zW } are considered during the rotation of the wind velocity vector, [22]. Assume the generalized coordinates of the mini UAV as q = (x, y, z, ψ, θ, φ)T ∈ R6 , where ξ = (x, y, z)T ∈ R3 represents the translation coordinates relative to the inertial frame, and η = (ψ, θ, φ)T ∈ R3 describes the vector of three Euler angles with rotations around z, y, x axes. These angles ψ, θ, and φ are called yaw, pitch and roll, respectively. Assume the translational velocity and the angular velocity in the body frame as ν = (u, v, w)T ∈ R3 and = ( p, q, r)T ∈ R3 , re-
ˆ is a skew-symmetric matrix such that where ˆ a = × a. The relation between the Euler angles and angular velocity is written as = X (η)η˙ where X (η) is called the Euler matrix and is represented as ⎛ ⎞ −sθ 0 1 X (η) = ⎝ cθ sφ cφ 0 ⎠ cθ cφ −sφ 0 Noting that this matrix X (η) presents a singularity when the pitch angle θ reaches ±90 deg, an inertial measurement unit was designed and built at the laboratory in order to specifically solve this singularity problem. For the aerodynamic analysis, two rotation matrices BSB and BW B are used. The rotation matrix S BB represents the transformation of a vector from the body frame B to the stability frame S , and the rotation matrix BW B describes the transformation of a vector from the body frame B to the wind frame W . These rotation matrices are written as ⎛
⎞ ⎛ ⎞ cα 0 sα cβ sβ 0 ⎝−sβ cβ 0⎠ BSB = ⎝ 0 1 0 ⎠ , BW B = −sα 0 cα 0 0 1
Fig. 2 Schematic of C-Plane II vehicle
where α is the angle of attack and β are the sideslip angle [6, 7, 20, 22].
126
J Intell Robot Syst (2012) 65:123–136
Therefore, the equations of motion of the aircraft are obtained using the Euler–Lagrange formulation ˙ ˙ d ∂ L (q, q) ∂ L (q, q) − =τ (1) dt ∂ q˙ ∂q where τ = (F, )T ∈ R6 denotes the forces and ˙ = moments acting on the body frame, L(q, q) ˙ − U (q) describes the Lagrangian equation K(q, q) ˙ which consists of the total kinetic energy K(q, q) and the potential energy U (q) of the system [8]. Moreover, K and U are defined as
4. 5. 6. 7.
˙ Sat(x) ≤ k4 ˙x
T 1 ˙ x 2 M(y) − C(q, z) x = 0 ∀x, y, z ∈ Rn ˙ ∀x, y ∈ Rn C(x, y) = 12 M(y) C(x, y)z ≤ kc yz ∀x, y, z ∈ Rn
where k1 , k2 , k3 , k4 , and kc are positive constants. In the stability analysis of this paper, the notation λmin {A} and λMax {A} are used to indicate the smallest and largest eigenvalues, respectively, of a symmetric positive definite matrix A(x), for n any x ∈ √R . The norm of vector x is defined as x = x T x and the induced norm of a matrix is also defined as A = λMax {AT A}.
K = Kt + Kr
2.2.2 Forces
U = mgz Kt =
1 ˙T ˙ ξ mξ 2 1 Kr = η˙ T M(η)η˙ 2
The forces acting on the aircraft include those of the propulsion system F p and aerodynamic effects F a . These forces are described as follows
where m ∈ R is the mass of the vehicle, Kt is the translational kinetic energy and Kr is the rotational kinetic energy with M(η) = X (η)T I X (η). I ∈ R3×3 denotes the moments of inertia of the mini VTOL aircraft. The saturation function Sat(x) is defined as follows
F = F p + Fa
⎧ x> ⎨ for Sat(x) = x for − ≤ x ≤ ⎩ − for x < −
(2)
with ⎛
⎞ Sat(x1 ) ⎜ ⎟ .. Sat(x) = ⎝ ⎠ . Sat(xn )
with ⎛
⎞ ⎛ ⎞ Tc −D B ⎝ F p = ⎝ 0 ⎠ , F a = BB Y ⎠ S BW 0 L where Tc is the thrust force of the two motors (Tc = T1 + T2 ). The lift force L, sideforce Y and drag force D are defined as aerodynamic forces, and g is the acceleration due to gravity, (see Fig. 2) [12]. In this analysis, the thrust force is oriented parallel to the axis xB of the body frame. 2.2.3 Moments
where x ∈ Rn and > 0. Now, some properties about the saturation function Sat(x), inertia matrix M(η), and coriolis ˙ need to be and centrifugal forces matrix C(η, η) introduced [9, 11]:
The moments generated on the mini VTOL aircraft are due to actuators (actuator moment act , reaction moment rot and gyroscopic moment gyro ), and the aerodynamic effects a . These moments are defined as follows
1. Sat(x) ≤ k1 x 2. Sat(x) ≤ k2 3. Sat(x)2 ≤ k3 Sat(x)T x
⎞
L = ⎝ M ⎠ = act + rot + gyro + a
N ⎛
J Intell Robot Syst (2012) 65:123–136
127
where
with ⎛
⎞
⎛
⎞
τφ Irot1 ω˙ r1 − Irot2 ω˙ r2 ⎠, act = ⎝ τθ ⎠ , rot = ⎝ 0 τψ 0 ⎛
gyro
⎞
1 T η˙ M(η)η˙ 2
(7)
From Eqs. 6 and 7
⎞ ¯ L 0 ¯ ⎠ = ⎝ r(Ir1 ωr1 − Ir2 ωr2 ) ⎠ , a = ⎝ M ¯ q(−Ir1 ωr1 + Ir2 ωr2 ) N
⎛
˙ (η) η˙ − M (η) η¨ + M
where τφ = a ( fe2 − fe1 ), τθ = e ( fe1 + fe2 ) and τψ = e fr are the control inputs with a and e that represent the distance from the center of mass to the forces fe1 and fe2 . ωri denotes the angular velocity of the rotor, Iri is the inertia moment of the propeller and Iroti is the moment of inertia of ¯ M ¯ and the rotor around its axis for i = 1, 2. L, N¯ are aerodynamic rolling, pitching and yawing moments, respectively (see Fig. 2) [12].
1 ∂ T η˙ M (η) η˙ = 2 ∂η
(8)
From Eq. 8, the Coriolis and Centrifugal vector is defined as ˙ η˙ = C(η, η)
1 ∂ T ˙ η˙ M(η) η˙ M(η) − 2 ∂η
(9)
Thus, the dynamic model for the rotational motion is rewritten as ˙ η˙ = M(η)¨η + C(η, η)
2.2.4 Translational and Rotational Dynamics Since the Lagrangian equation contains no cross˙ terms in the kinetic energy combining ξ˙ with η, the Euler–Lagrange Eq. 1 can be partitioned into dynamics for ξ coordinates and η coordinates [3]. The translational motion can be obtained using the following expression d ∂ Lt ∂ Lt =F (3) − dt ∂ ξ˙ ∂ξ with 1 T Lt = ξ˙ mξ˙ − mgz 2
Lr =
(10)
˙ are represented as where M(η) and C(η, η)
M(η) ⎛ =⎝
Ixx s2θ + I yy c2θ s2φ + Izz c2θ c2φ cθ cφ sφ (I yy − Izz ) −Ixx sθ cθ cφ sφ (I yy − Izz )
I yy c2θ + Izz s2φ
0
−Ixx sθ
0
Ixx
⎞ ⎠
(11) ⎛
(4)
Thus, after some computations, the translation motion of this vehicle is described as ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ mx¨ Tc −D B ⎝ ⎝ m y¨ ⎠ = RIB ⎝ 0 ⎠ + RIB BB Y ⎠ S BW m¨z 0 L ⎛ ⎞ 0 +⎝ 0 ⎠ (5) −mg Similarly, the rotational motion is described as d ∂ Lr ∂ Lr − = (6) dt ∂ η˙ ∂η
⎞ c11 c12 c13 ˙ = ⎝ c21 c22 c23 ⎠ C(η, η) c31 c32 c33
(12)
with ˙ θ cθ + I yy −θs ˙ θ cθ s2φ + φc ˙ 2θ sφ cφ c11 = Ixx θs ˙ 2θ sφ cφ −Izz θ˙ sθ cθ c2φ + φc ˙ θ cθ c12 = Ixx ψs ˙ θ s2φ − φc ˙ θ c2φ + ψs ˙ θ cθ s2φ −I yy θ˙ sθ sφ cφ + φc ˙ θ c2φ − ψs ˙ θ cθ c2φ + θs ˙ θ sφ cφ ˙ θ s2φ − φc +Izz φc ˙ 2θ sφ cφ − Izz ψc ˙ 2θ sφ cφ c13 = −Ixx θ˙ cθ + I yy ψc
128
J Intell Robot Syst (2012) 65:123–136
˙ θ cθ + I yy ψs ˙ θ cθ s2φ + Izz ψs ˙ θ cθ c2φ c21 = −Ixx ψs ˙ φ cφ + Izz φs ˙ φ cφ c22 = −I yy φs ˙ φ cφ + ψc ˙ θ c2φ − ψc ˙ θ s2φ ˙ θ + I yy −θs c23 = Ixx ψc ˙ θ s2φ − ψc ˙ θ c2φ + θs ˙ φ cφ +Izz ψc ˙ 2θ sφ cφ + Izz ψc ˙ 2θ sφ cφ c31 = −I yy ψc ˙ θ + I yy θs ˙ φ cφ + ψc ˙ θ s2φ − ψc ˙ θ c2φ c32 = −Ixx ψc ˙ θ s2φ − ψc ˙ θ c2φ + θs ˙ φ cφ −Izz ψc c33 = 0
(13)
2.3 Aerodynamic Analysis This section presents the aerodynamic analysis of the mini VTOL aircraft in horizontal flight. For simplicity, a two-dimensional (2-D) scheme of C-Plane II vehicle is considered to study the aerodynamic effects in presence of wind. Since the stability of this mini VTOL aircraft is affected by the propulsion system, the propeller effect (Propeller Momentum Theory) is analyzed to obtain the aerodynamic behavior of the vehicle in horizontal flight [14, 15, 19]. In Fig. 3, αs is the slipstream angle of attack, α denotes the wing angle of attack, α p describes
the propeller angle of attack, θ is pitch angle and δ expresses the angle of the ailerons. vo is the freestream velocity, vi represents the induced velocity and is directed opposite to the thrust, vs is the slipstream velocity, vds denotes the far downstream velocity, and finally, vr is the resultant velocity in the propeller slipstream. This Figure shows the wing-propeller combination where the wing is submerged in the air slipstream, and the C-Plane II vehicle is turned in the direction of the air freestream. From Fig. 3, the following aerodynamic equations can be established
vr =
2 2 2vi + vo cos α p + vo sin α p
α = tan−1
vo sin α p 2vi + vo cos α p
αs = sin−1
(14)
(15)
2vi sin α p vr
(16)
Taking, the thrust equation Tc and the downstream velocity vds Tc = ρ Avs vds − vo cos α p
(17)
and vds = vo cos α p + 2vi
(18)
where A represents the area of the rotor disc, ρ denotes the air density. Now, considering the case where air freestream velocity vo is different from zero, then the resultant velocity vr is equal to slipstream velocity vs far behind of the propeller, Eqs. 14, 17 and 18 are rewritten in a quartic form for the induced velocity as
Fig. 3 Longitudinal motion of the C-Plane II UAV
4vi4 + 4vo vi3 cos α p + vo2 vi2 =
Tc 2ρ A
2 (19)
J Intell Robot Syst (2012) 65:123–136
and vh =
129
namic moments. Since the roll motion is mechanically stable using contra-rotating propellers, the reaction and gyroscopic moments ( rot and gyro ) will essentially be zero. The nonlinear control law is described as t ˙˜ ˜ − Ki ˜ γ Sat (η(s)) + η(s) ds act = −K p Sat (η)
Tc 2ρ A
consequently Eq. 19 yields 4 3 vi vi vo 4 +4 cos α p vh vh vh 2 2 vo vi + =1 vh vh
0
(20)
Finally, the aerodynamic forces and moments are written by D = Y = L = ¯ = L ¯ = M N¯ =
1 2 ρv SC D 2 r 1 2 ρv SCY 2 r 1 2 ρv SC L 2 r 1 2 ρv Sb Cl 2 r 1 2 ¯ m ρv ScC 2 r 1 2 ρv Sb Cn 2 r
− K d Sat η˙˜ + M (η) η¨ d ˙ η˙ d − a + C(η, η)
(22)
where K p , K i , K d , are the symmetric definite positive matrices, and are proportional, integral and derivative gains of the controller, and γ > 0. The control for trajectory tracking can be given as the following. Let ηd and η˙ d denote the desired trajectories and define the error functions η˜ 1 = η − ηd η˜ 2 = η˙ − η˙ d
(21)
where S represents the wing area, c¯ is the wing chord, and b is the wing span. C D , CY and C L are aerodynamical non-dimensional coefficients of drag, sideforce and lift. Cl , Cm and Cn are aerodynamical non-dimensional coefficients of rolling, pitching and yawing moments.
3 Control Technique This section presents the control strategy for the attitude stabilization of the C-Plane II UAV in horizontal flight. A nonlinear Saturated Proportional-Integral-Derivative (SPID) control with compensation of aerodynamic moments is proposed in order to stabilize the vehicle during the flight. The stability analysis is based on the Lyapunov theory. The controller is capable of preventing aggressive spiral and phugoid motion during cruise, due to the compensation of aerody-
(23)
In this stability analysis, the variables η = η1 , η˙ = η2 , η˜ = η˜ 1 and η˙˜ = η˜ 2 are used. The integral action of the control law in Eq. 22 adds a state variable η0 in the system [10]. Using Eqs. 10 and 22, the closed-loop system is writ ten as η˜˙ 0 = γ Sat η˜ 1 + η˜ 2 η˙˜ 1 = η˜ 2
η˙˜ 2 = M(η1 )−1 −K p Sat η˜ 1 − K i η˜ 0 −K d Sat η˜ 2 − C(η1 , η2 )η˜ 2
(24)
In order to achieve the asymptotic stability of the system in Eq. 24, a Lyapunov candidate function is proposed as 1 1 V t, η˜ 0 , η˜ 1 , η˜ 2 = η˜ 2T M (η1 ) η˜ 2 + η˜ 1T K p η˜ 1 2 2 T 1 T + η˜ 0 K i η˜ 0 +γ Sat η˜ 1 M (η1 ) η˜ 2 2 (25) Notice that, 2 1 1 T η˜ 2 Mη˜ 2 ≥ λmin {M} η˜ 2 2 2
(26)
130
J Intell Robot Syst (2012) 65:123–136
2 1 T 1 η˜ 1 K p η˜ 1 ≥ λmin K p η˜ 1 2 2 2 1 T 1 η˜ 0 K i η˜ 0 ≥ λmin K i η˜ 0 2 2
(27)
(28)
Computing the time derivative of Eq. 25 and employing properties 5 and 6, it gets V˙ t, η˜ 0 , η˜ 1 , η˜ 2 = −η˜ 2T K d Sat η˜ 2 ˙ η˜ 1 T M (η1 ) η˜ 2 + γ Sat
and from property 1, the last term of Eq. 25 yields
T − γ Sat η˜ 1 K p Sat η˜ 1
T γ Sat η˜ 1 Mη˜ 2 ≤ γ k1 λMax {M} η˜ 1 η˜ 2 T −γ Sat η˜ 1 Mη˜ 2 ≥ −γ k1 λMax {M} η˜ 1 η˜ 2
T − γ Sat η˜ 1 K d Sat η˜ 2 (29)
T + γ Sat η˜ 1 C (η1 , η2 ) η˜ 2
Taking Eqs. 26, 27, 28, and 29, Eq. 25 yields
In order to demonstrate that V˙ < 0, each term of Eq. 32 is analyzed. Using properties 1 and 3, it can be written as
2 1 V t, η˜ 0 , η˜ 1 , η˜ 2 > λmin {M} η˜ 2 2 2 1 + λmin K p η˜ 1 2 2 1 + λmin {K i } η˜ 0 2 −γ k1 λMax {M} η˜ 1 η˜ 2
η˜ 2T K d Sat η˜ 2
(30) The Eq. 30 can be written as ⎡ ˜ ⎤T ⎡ η˜ 0 ⎤ 1 η0 ⎦ Q ⎣ η˜ 1 ⎦ η˜ 1 V t, η˜ 0 , η˜ 1 , η˜ 2 > ⎣ 2 η˜ 2 η˜ 2 where ⎡ ⎢ ⎢ Q=⎢ ⎣
(32)
≥
2 1 λmin {K d } Sat η˜ 2 k3
2 1 −η˜ 2T K d Sat η˜ 2 ≤ − λmin {K d } Sat η˜ 2 k3 ≤−
(31)
2 1 2 k1 λmin {K d } η˜ 2 k3
Considering property 4 for the second term of Eq. 32, it results ˙ η˜ 1 T Mη˜ 2 ≤ γ k4 λMax {M} η˜ 2 2 γ Sat
⎤
λmin {K i } 0 0 0 λmin K p −γ k1 λMax {M} ⎥ ⎥ 0 −γ k1 λMax {M} λmin {M} ⎥ ⎦
Now, applying the theorem of Sylvester [18], it implies λmin {K i } > 0 λmin {K i } λmin K p > 0 γ 2 k21 λMax {M}2 λmin K p > λmin {M} Therefore, V > 0 is a positive definite function.
(33)
(34)
Employing the property 1 in the third and fourth terms of Eq. 32, it leads T γ Sat η˜ 1 K p Sat η˜ 1
2 ≥ γ k21 λmin K p η˜ 1
T 2 −γ Sat η˜ 1 K p Sat η˜ 1 ≤ −γ k21 λmin K p η˜ 1 (35) T γ Sat η˜ 1 K d Sat η˜ 2 ≥ γ k21 λmin {K d } η˜ 1 η˜ 2 T −γ Sat η˜ 1 K d Sat η˜ 2 ≤ −γ k21 λmin {K d } η˜ 1 η˜ 2 (36)
J Intell Robot Syst (2012) 65:123–136
131
Finally, using properties 1, 2 and 7 for the last term of Eq. 32, it yields T γ Sat η˜ 1 Cη˜ 2 ≤ γ kc η˜ 2 η2 Sat η˜ 1 ≤ γ kc η˜ 2 η˜ 2 + ηd2 Sat η˜ 1 2 ≤ γ k2 kc η˜ 2 +γ kc ηd2 η˜ 2 Sat η˜ 1 2 ≤ γ k2 kc η˜ 2 +γ k1 kc ηd2 η˜ 2 η˜ 1
(37)
Then, employing Eqs. 33, 34, 35, 36 and 37, the Eq. 32 is rewritten as 2 1 V˙ t, η˜ 0 , η˜ 1 , η˜ 2 ≤ − k21 λmin {K d } η˜ 2 k3 2 +γ k4 λMax {M} η˜ 2 2 −γ k21 λmin K p η˜ 1 −γ k21 λmin {K d } η˜ 1 η˜ 2 2 +γ k2 kc η˜ 2 +γ k1 kc ηd2 η˜ 2 η˜ 1 (38) The Eq. 38 can be written as $ %T $ % η˜ 1 η˜ 1 V˙ t, η˜ 0 , η˜ 1 , η˜ 2 ≤ − P η˜ 2 η˜ 2
(39)
where ⎡ ⎢ P=⎢ ⎣
γ k21 λmin K p γ k21 2
λmin {K d } −
γ k1 kc ηd2 2
γ k21 γ k1 kc ηd2 λmin {K d } − 2 2 1 2 k λmin {K d } − γ k4 λMax {M} − γ k2 kc k3 1
Applying the theorem of Sylvester [18], it leads
λmin {K p }
$
> γ k21
%2 γ k1 kc λmin {K d } − ηd2 2 2
E = η˜ 1 ∈ R3 , η˜ 2 ∈ R3 | V˙ = 0
γ k21
1 2 k λmin {K d } − γ k4 λMax {M} − γ k2 kc k3 1
⎥ ⎥ ⎦
Then, V˙ ≤ 0 is a negative semidefinite function. Now, using the theorem of LaSalle in Eq. 39 [10], it implies
λmin {K p } > 0
$
⎤
%
(40)
Therefore, the closed-loop system is asymptotically stable about the origin. Figure 4 illustrates the block diagram of the closed-loop system of this vehicle.
4 Simulation Results -
SPID CONTROLLER
+ +
Fig. 4 Closed-loop system
UAV
In this section, the results obtained by simulation are presented in order to observe the performance in horizontal flight of the proposed control law.
132
J Intell Robot Syst (2012) 65:123–136 40
ñ
1
20
ñ0
10 0 −10 −20 5
10
15
0 −50 dñ2/dt
−100
dñ1/dt dñ /dt 0
−150 0
5
10
time [s]
time [s]
(a)
(b)
0.5
12
0
10 L/D ratio
Control law [Nm]
−30 0
Fig. 6 Control law and L/D ratio
50
ñ2
30
velocity errors [deg/s]
Position errors [deg]
Fig. 5 Angular position and velocity of the vehicle
−0.5 −1
8 6 4
−1.5
2
−2 0
0 0
time [s]
15 α [deg]
(a)
(b)
5
10
15
5
10
20
Fig. 8 Microcontroller
1 Aerodynamic coefficients
CD CL
0.8 0.6 0.4 0.2 0
0
5
10
15 α [deg]
15
20
25
30
Fig. 7 Aerodynamic coefficients C D and C L
Table 1 Parameters of the mini VTOL aircraft Parameter
Value
Airfoil shape ¯ Wing chord (c) Wing span (b ) Aspect ratio (AR) Thickness Wing area (S) Vehicle mass (m) Ixx I yy Izz
Flat plate 0.5 m 0.80 m 1.6 0.006 m 0.4 m2 0.8 kg 1.0775 kg/m2 0.4374 kg/m2 0.06399 kg/m2
Fig. 9 The homemade inertial measurement unit
25
30
J Intell Robot Syst (2012) 65:123–136
133
High pass filter
VTOL aircraft was published in [12]). This section also shows practical results during cruise. Experiments of this vehicle were realized in autonomous and semi-autonomous horizontal flight. In the semi-autonomous case, the pilot only provided the throttle control.
Rate gyro
+
+ Accelerometer
−
output Angle estimation
accel g
5.1 Experimental platform
Low pass filter
Fig. 10 Block diagram of the analog complementary filter
The airframe of C-Plane II vehicle is built from a polystyrene foam sheet (Depron) and carbon fiber-tubes. The vehicle is powered by two brushless motors driving contra-rotating propellers 10 × 4.5 in (Coaxial configuration). Table 1 shows some parameters of this vehicle. The embedded system used in this VTOL vehicle consists of a microcontroller RCM3400 having a Rabbit microprocessor which receives the data from sensors, computes the control law and then sends the control signals to actuators in order to obtain a suitable control (see Fig. 8). The Homemade Inertial Measurement Unit (HIMU) is built using three rate gyros (ADXRS150 Analog Devices), three accelerometers (ADXL203 Analog Devices), and two magnetometers (KMZ52 Philips), see Fig. 9. These sensors are placed in orthogonal positions in order to measure the aircraft motion. The main characteristic of this HIMU is that it can operate in both hover and horizontal flight solving the singularity problem. In addition, this HIMU incorporates an analog complementary filter and signal conditioning circuits in order to estimate the angle measurements (see Fig. 10). These analog filter and signal conditioning circuits are designed using active components (precision
Figure 5 shows the angular position and velocity errors during cruise. Observe that, the error signals converge to zero and the closed-loop system is asymptotically stable. The control law that stabilize the mini VTOL aircraft during the horizontal flight is depicted in the Fig. 6a. Figure 6b also illustrates the lift-to-drag ratio for some values of α. The aerodynamic coefficients C D and C L are shown in the Fig. 7. The parameters involved in the control of this vehicle are taken as k1 = 1, k2 = 1, k3 = 1, k4 = 1, kc = 1, γ = 1, and the matrices are defined as K p = diag{20, 20, 20}, K d = diag{12, 12, 12} and K i = diag{15, 15, 15}.
5 Experimental Results This section describes the embedded system and the mini VTOL aircraft (C-Plane II) developed at the Laboratory LAFMIA UMI 3175 CNRS Cinvestav Mexico (the first version of this mini
Fig. 11 Angle estimation in low and high frequencies
50
15 Angle estimation Accelerometer
10 Degrees
Degrees
Angle estimation Accelerometer
0
5 0 −5 −10
−50
0
5
10 time [s]
(a)
15
−15 0
5
10 time [s]
(b)
15
20
134
J Intell Robot Syst (2012) 65:123–136
Fig. 12 Block diagram of the embedded system
Power Interface
RC RECEIVER UAV
SPID CONTROLLER
XBEE MODEM
Microcontroller
Homemade Inertial Measurement Unit PITOT TUBE
ACCELEROMETERS
ANALOG
SIGNAL
COMPLEMENTARY
CONDITIONING
FILTER
CIRCUIT
M
Actuators
+
RATE GYROS
MAGNETOMETERS BRUSHLESS MOTORS
Fig. 13 Roll angle and control
10 Roll control [PWM]
Roll angle [deg]
10 5 0 −5
−10 0
20
40
60
80 100 120 140 160 180 time [s]
5 0 −5
−10 0
20
40
60
(a)
5 Pitch control [PWM]
5 0 −5 −10 0
20
40
60
0
−5 0
80 100 120 140 160 180 time [s]
20
40
60
(a)
Fig. 15 Yaw angle and control
6 Yaw control [PWM]
Yaw angle [deg]
80 100 120 140 160 180 time [s]
(b)
10 5 0 −5
−10
80 100 120 140 160 180 time [s]
(b)
10 Pitch angle [deg]
Fig. 14 Pitch angle and control
RC SERVOMOTORS
0
20
40
60
80 100 120 140 160 180 time [s]
(a)
4 2 0 2 −4 −6 0
20
40
60
80 100 120 140 160 180 time [s]
(b)
J Intell Robot Syst (2012) 65:123–136
operational amplifiers and integrated circuits) and passive components (resistors, capacitors, and inductors). Estimation responses in real-time of the analog filter complementary are illustrated in the Fig. 11. This responses are obtained in low and high frequencies [23]. On the other hand, one magnetometer measures the magnetic field in hover flight, while the other operates in the horizontal flight. For computing the airspeed of the mini VTOL aircraft in horizontal flight, a small pitot-static tube is added to the embedded system. This pitotstatic tube measures the difference between static and total pressures, allowing the freestream airspeed to be determined. Figure 12 shows the embedded system developed for the C-Plane II vehicle. Figures 13, 14 and 15 show the real-time responses of the vehicle orientation. Figure 13 demonstrates the effective behavior in roll motion. Pitch and yaw motion are illustrated in Figs. 14 and 15 respectively.
6 Conclusions and Future Work 6.1 Conclusions A C-Plane II aircraft has been presented in order to perform autonomous horizontal flight. The dynamic model of the vehicle has been obtained using the Euler–Lagrange formulation, and the aerodynamic terms including the wind are described by considering the propeller analysis. A nonlinear saturated Proportional-IntegralDerivative (SPID) controller with compensation of aerodynamic moments has been proposed for stabilizing the vehicle in horizontal flight. The stability analysis based on Lyapunov theory has demonstrated that the closed-loop system is asymptotically stable in the origin. In addition, a homemade inertial measurement unit (HIMU) has been developed in order to operate in both vertical and horizontal flight avoiding the singularity problem. The sensor fusion of this HIMU is obtained by using an analog complementary filter. Finally, simulations and experimental results of the closed-loop system have shown an effective behavior of the proposed control law.
135
6.2 Future work Future work for C-Plane II UAV includes incorporating additional systems like GPS and visionbased sensors to navigate autonomously. The aim of this additional instrumentation is to achieve a complete navigation flight.
References 1. Barkana, I.: Classical and simple adaptive control for nonminimum phase autopilot design. AIAA J. Guid. Control Dyn. 28(4), 631–638 (2005) 2. Blight, J.D., Gangsaas, D., Richardson, T.M.: Control law synthesis for an airplane with relaxed static stability. AIAA J. Guid. Control Dyn. 9(5), 546–554 (1986) 3. Castillo, P., Lozano, R., Dzul, A.: Modelling and Control of Miniflying Machines. Springer, London (2005) 4. Collinson, R.P.G.: Introduction to Avionics. Chapman & Hall, London (1996) 5. Cohen, K., Bossert, D.E.: Fuzzy Logic Non-minimum Phase Autopilot Design. AIAA Paper, pp. 2003–5550, Seoul (2003) 6. Cook, M.V.: Flight Dynamics Principles, 2nd edn. Butterworth-Heinemann, Oxford (2007) 7. Etkin, B., Reid, L.D.: Dynamics of Flight: Stability and Control, 3rd edn. Wiley, New York (1996) 8. Goldstein, H., Poole, C.P., Safko, J.L.: Classical Mechanics, 2nd edn. Addison Wesley Series in Physics, Addison-Wesley, Reading (1983) 9. Kelly, R.: Global Positioning of robot manipulators via PD control plus a class of nonlinear integral actions. IEEE Trans. Automat. Contr. 43(7) 934–938 (1998). doi:10.1109/9.701091 10. Khalil, H.: Nonlinear Systems, 3rd edn. Prentice Hall, Englewood Cliffs (2002) 11. Koditschek, D.E.: Natural motion for robot arms. In: IEEE Conference on Decision and Control (CDC 1984), Las Vegas, NV (1984) 12. Lozano, R.: Unmanned Aerial Vehicles Embedded Control. John Wiley-ISTE Ltd (2010) 13. Magni, J.F., Bennani, S., Terlouw, J.: Robust Flight Control: a Design Challenge. Lecture Notes in Control and Information Sciences, number 224. Springer, New York (1997) 14. McCormick, B.W.: Aerodynamics of V/STOL Flight. Dover, New York (1967) 15. McCormick, B.W.: Aerodynamics, Aeronautics and Flight Mechanics. Wiley, New York (1995) 16. McLean, D.: Automatic Flight Control Systems. Prentice Hall, Englewood Cliffs (1990) 17. Nelson, R.C.: Flight Stability and Automatic Control, 2nd edn. McGraw-Hill, New York (1998) 18. Ogata, K.: Modern Control Engineering, 5th edn. Prentice Hall, Englewood Cliffs (2009) 19. Phillips, W.F.: Mechanics of Flight. John Wiley and Sons, USA (2004)
136 20. Stengel, R.F.: Flight Dynamics. Princeton University Press, Princeton (2004) 21. Sharma, V., Voulgaris, P.G., Frazzoli, E.: Aircraft Autopilot Analysis and Envelope Protection for Operation Under Icing Conditions. AIAA J. Guid. Control Dyn. 27(3), 454–465 (2004)
J Intell Robot Syst (2012) 65:123–136 22. Stevens, B.L., Lewis, F.L.: Aircraft Control and Simulation, 2nd edn. Wiley-Interscience, New York (2003) 23. Zimmermann, M., Sulzer, W.: High bandwidth orientation measurement and control based on complementary filtering. In: Proceeding of SYROCO’91, IFAC Symposium on Robot Control, Vienna, Austria (1991)
