AIAA 2005-6266
AIAA Guidance, Navigation, and Control Conference and Exhibit 15 - 18 August 2005, San Francisco, California
Design of a Nonlinear Autopilot for Velocity and Attitude Control using Block Backstepping John W.C. Robinson∗ and Ulrik Nilsson† Swedish Defence Research Agency, Stockholm, SE 172 90, Sweden We present a solution to the three-dimensional attitude and velocity control problem for an autopilot based on backstepping control of the full six-degrees-of-freedom equations of motion. The resulting controller has a very simple structure and requires only a limited amount of a priori information about the aerodynamic properties of the aircraft. It has global stabilizing properties in attitude and local stabilizing properties in velocity. A key element of the controller is the use of a spherical linear interpolation (slerp) to compute a geodesic on the manifold of unit norm quaternions in order to achieve the minimal rotation required to control the attitude. We give a proof of stability of the resulting controller and illustrate its behavior using a version of the ADMIRE model which is a realistic model of a small single engine fighter.
I.
Introduction
n recent years a number of design methodologies for nonlinear design of control systems for aircraft and Iinversion missiles have emerged, including nonlinear dynamic inversion and backstepping. Nonlinear dynamic has gained interest mainly because it is simple to apply and can be combined with established 1–3
4
knowledge from linear design methods. Backstepping has showed great potential, and has also been applied to a number of different design tasks, but has the drawback that a Lyapunov function has to be found in order to apply the method. Here we are going to apply multi-input backstepping, also known as block backstepping, to the problem of controlling the attitude and velocity of an aircraft described by the full nonlinear six-degree-of-freedom dynamics. Part of this problem has been treated before, e.g. by Steinicke and Michalka5 in the context of trajectory control of a missile. Recently, Glad and H¨ arkeg˚ ard6 have given a solution to a variant of the velocity and angular velocity control problem in the context of flight control using multi-input backstepping. They assume, however, that the thrust force is always aligned along the velocity vector which is not the case in general. a Our solution employs a velocity control similar to that of Glad and H¨ arkeg˚ ard but extends it to the case of nonaligned thrust and combines it with a quaternion based control for attitude control. The quaternion based attitude control utilizes spherical linear interpolation (slerp) on the sphere S 3 to compute a geodesic representing the minimal rotation of the body needed to control the attitude to the desired value.
II.
Equations of Motion
The equations of motion for the aircraft are given by the Newton and Euler equations for rigid body motion, formulated around the center of gravity, and combined with the standard quaternion differential ∗ Deputy
Research Director, Department of Autonomous Systems. Department of Autonomous Systems. a After acceptance of this paper the authors were made aware by Dr. H¨ arkeg˚ ard of a manuscript in preparation in which the thrust force alignment requirement used by Glad and H¨ arkeg˚ ard6 is removed. The solution presented is limited to control of the body states however, and does not cover the combined attitude and velocity control problem considered here. † Scientist,
1 of 20 American Institute of Aeronautics and Astronautics Copyright © 2005 by John W.C. Robinson. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.
Nomenclature V v R Q ω f ex τ t = τ ex uτ b m g j m [Cx , Cy , Cz ]T ρ Sref cref bref k·k η˙ q = 21 ρkVk2 P[v0 ] (v), P[v0 ]⊥ (v) Subscripts 0 Superscripts T des
Velocity vector in Earth fixed (Inertial) frame Fe Velocity vector in body fixed (aircraft) frame fb Rotation matrix for transformation fb → Fe Orientation quaternion corresponding to R Angular velocity vector in body frame fb Aerodynamic force vector in body frame fb Unit vector along x-axis (ex = [1, 0, 0]T ) Engine thrust force magnitude Thrust force vector in body frame fb Commanded thrust force magnitude Engine thrust force time constant Moment vector around center of gravity in body frame fb Gravitational acceleration vector in body frame fb Moment of inertia matrix around center of gravity in body frame fb Total mass of aircraft Vector of aerodynamic force coefficients (for the x,y, z axes in the body frame fb ) Air density Reference wing area Reference wing chord Reference wing span Vector 2-norm and associated (induced) matrix 2-norm Time derivative of a variable η Dynamic pressure Projection of the vector v onto the subspace spanned by the vector v0 and projection of v onto the plane orthogonal to v0 , respectively Reference values Transposition of a vector Virtual controls
equation for orientation,7 viz. V
=
˙ Q
=
f +t = m =
Rv, 1 ˜ Q ◦ ω, 2 m(v˙ + ω × v − g), jω˙ + ω × jω,
(1) (2) (3) (4)
˜ = (0, ω) denotes the pure quaternion formed from the body angular velocity vector ω and ◦ denotes where ω quaternion multiplication. The aerodynamic force f = [f1 , f2 , f3 ]T and moment m = [m1 , m2 , m3 ]T vectors
2 of 20 American Institute of Aeronautics and Astronautics
are defined in terms of standard aerodynamic coefficients as f1
=
qSref Cx ,
(5)
f2 f3
= =
qSref Cy , qSref Cz ,
(6) (7)
m1 m2
= =
qSref bref Cl , qSref cref Cm ,
(8) (9)
m3
=
qSref bref Cn , .
(10)
The engine dynamics are modeled using a simple first order linear system as τ˙ = b(τ − uτ ), A.
b > 0.
(11)
Deviations from reference values
In order for (1)–(4) and (11) to be useful for our further developments we must make a change of variables and rewrite these equations in terms of deviations from a reference point. Let [v0 , Q0 , ω 0 , τ0 ]T
(12)
be a vector of reference values for the state variables in (2)–(4) and (11), and let [v, Q, ω, τ ]T
(13)
be the vector of deviations from the reference values. We shall study two types of reference vectors (12). The first type is a reference vector corresponding to an equilibrium to the state vector in (2)–(4) and (11) (i.e. the entire vector in (12) is constant). In this case we must necessarily have ω 0 = 0, which corresponds to straight path flight. The second type of reference vector is the one obtained when we have an equilibrium only for the states in (3), (4) and (11) (i.e. [v0 , ω0 , τ0 ]T is constant), and ω0 is constant but nonzero (so that Q0 is time varying b ). This corresponds to a constant-g turn. In this case we shall assume that the reference trajectory for Q0 is the solution to (2) corresponding to ω 0 . 1.
Change of variables.
From now on we assume that we have made the change of variables v v + v0 Q Q + Q0 → ω ω + ω0 τ τ + τ0
in the equations (2)–(4) and (11), as outlined above. We moreover assume that the aerodynamic forces f are mainly dependent on v + v0 , the aerodynamic moments m are mainly dependent on v + v0 , ω + ω 0 , and the thrust t acts only along the aircraft x-axis (i.e. the x-axis in the body frame fb ). We can make this dependence explicit by writing f = f (v + v0 ),
m = m(v + v0 , ω + ω 0 ),
t = (τ + τ0 )ex .
Likewise, the gravity vector g is only dependent on Q + Q0 and therefore we can write g = g(Q + Q0 ). ˜ ˜ (Q, Q0 ) by It will be convenient to introduce the vector functions ˜f (v, v0 ), m(v, v0 , ω, ω0 ), and g ˜f (v, v0 ) ˜ m(v, v0 , ω, ω0 ) ˜ (Q, Q0 ) g
= f (v + v0 ) − f (v0 ),
= m(v + v0 , ω + ω 0 ) − m(v0 , ω 0 ), = g(Q + Q0 ) − g(Q0 ),
(14) (15) (16)
b The main stability result below covers, as it stands, only case of a constant Q but the theory has practical applicability 0 also to “sufficiently slowly” time varying cases, as can be seen from the simulation example.
3 of 20 American Institute of Aeronautics and Astronautics
and the skew symmetric matrix function C taking values in R3×3 representing the cross product, so that e.g. 0 −v3 v2 ω1 C(v)ω = v3 (17) 0 −v1 ω2 = v × ω. −v2 v1 0 ω3
Moreover, it will be convenient to introduce also a matrix–vector representation for the product of two quaternions. Let Q1 = (a1 , b1 ), Q2 = (a2 , b2 ) be two quaternions with real parts a1 , a2 ∈ R and imaginary parts b1 , b2 ∈ R3 , respectively. Then the quaternion product Q1 ◦ Q2 can be written in terms of an ordinary matrix–vector product as Q1 ◦ Q2 = T(Q1 )Q2 , where the matrix T(Q1 ) is given by T(Q1 ) =
"
a1 b1
#
−bT1 C(b1 ) + a1 I3×3
,
(18)
and C is the matrix in (17) giving the vector product. In the special case that Q2 is a pure quaternion, i.e. a0 = 0, we have T(Q1 )Q2 = B(Q1 )b2 , where the matrix function B is defined by B(Q1 ) =
"
−bT1 C(b1 ) + a1 I3×3
#
.
(19)
For later reference we also note that the product B(Q1 )b2 can be written B(Q1 )b2 = A(b2 )Q1 where the skew symmetric matrix function A is given by " # 0 −bT2 A(b2 ) = . b2 −C(b2 )
(20)
As an application of the representation (20) we recall that since the reference value Q0 is assumed to satisfy (2) when driven by ω0 we must have ˙ 0 = 1 A(ω 0 )Q0 . Q (21) 2 2.
Equations of motion for attitude & velocity control.
Now, if we recall the assumption that a reference point as in (12) is always at least an equilibrium point to the equations (3), (4) and (11) we can, using (14)–(21), rewrite the entire set of state equations (2)–(4) and (11) equivalently as c v˙ ˙ Q ω˙ τ˙
1˜ 1 ˜ (Q, Q0 ) + v × ω 0 + C(v + v0 )ω + τ ex , f (v, v0 ) + g m m 1 1 = A(ω 0 )Q + B(Q + Q0 )ω, 2 2 ˜ = j−1 m(v, v0 , ω, ω 0 ) − (ω + ω 0 ) × jω − ω × jω 0 , = bτ − b˜ uτ ,
=
(22) (23) (24) (25)
where we have used the fact that many terms cancel or vanish at an equilibrium and also introduced u ˜τ as u ˜τ = uτ − τ0 . c Note
that there is no linearization involved in (22)–(25), it is still the full nonlinear equations (2)–(4) and (11) albeit rewritten in terms of deviations from a reference and with omission of terms that sum to zero.
4 of 20 American Institute of Aeronautics and Astronautics
For control of the system (22)–(25) we are going to assume that all states in are directly measurable so that we can employ state feedback control. Moreover, we are going to assume that we can control the ˜ three components of the moment vector m(v, v0 , ω, ω0 ) directly (which is equivalent to controlling the original moment vector components m1 , m3 , m3 ) and that it is the responsibility of the controller to perform the “inversion” from moments to actual control surface commands. Furthermore, we are going to neglect actuator dynamics.
III. A.
A Nonlinear Autopilot
The Simultaneous Attitude and Velocity Control Problem
The three-dimensional attitude-velocity control problem can be cast as the problem of controlling all three components of the aircraft velocity vector V in the frame Fe . If this is done with regard only to the relations between the velocity vector components V1 , V2 , V3 we obtain a natural three-dimensional generalization of the standard flight path angle control problem. However, we are here going to also consider simultaneous control of the magnitude kVk of the velocity vector V in Fe . If we recall the relation (1) between the velocity in the frames Fe and fb respectively as V = Rv, it is clear that the attitude-velocity control problem can be split into two separate problems; (i) the problem of controlling the body components (i.e. in fb ) of the airspeed vector v and (ii) simultaneously controlling the orientation in terms of R of the aircraft (in Fe ). Since we only have to our disposal as control inputs the three moment vector components m1 , m2 , m3 and the thrust command τ it is clear that this is an underactuated control problem. B.
Integrator Backstepping
The simplest form of multi-input backstepping deals with controlled dynamical systems of the generic affine form ˙ x(t) = f (x(t)) + g(x(t))u(t), (26) where f : Rn → Rn and g : Rn → Rn × Rm are Lipschitz continuous u(t) is a vector in Rm of continuous control functions. To make things simple we are going to assume that both f, g are smooth. The object is then to find ℓ : Rn → Rm such that if u in (26) is given by u = ℓ(x)
(27)
˙ x(t) = f (x(t)) + g(x(t))ℓ(x(t))
(28)
the resulting system is stable, in some suitable sense. In integrator backstepping this is accomplished by first augmenting the system with integrators on the inputs so that the system (26) transforms into ˙ x(t) ˙ ξ(t)
= f (x(t)) + g(x(t))ξ(t), = u(t).
(29) (30)
The following step can intuitively be described as finding a control function u to the augmented system (29),(30) such that when it is applied the state ξ behaves as “if it were” a suitable feedback function ℓ as in (27). This desired behavior is called a virtual control and closely associated to it is the concept of control Lyapunov function. A smooth positive definite and radially unbounded function V : Rn → [0, ∞) is called a control Lyapunov function (CLF) for the system (26) if it holds that infm
u∈R
∂V (x) f (x) + g(x)u < 0, ∂x
∀x 6= 0.
The basic assumption about the system (30) that we need in order to actually apply integrator backstepping is one about stabilizability.
5 of 20 American Institute of Aeronautics and Astronautics
Assumption III.1. Consider the system in (26) and assume that there exist a smooth feedback law ℓ as in (27) and a smooth, positive definite and radially unbounded function V : Rn → [0, ∞) such that ∂V (x) f (x) + g(x)ℓ(x) ≤ −W (x) ≤ 0, ∂x
∀x ∈ Rn ,
(31)
for some continuous W : Rn → R which is positive definite.
Under this assumption one can formulate the basic result for integrator backstepping as follows.
Theorem III.2. Consider the system (29), (30) and suppose that (29) satisfies Assumption III.1 with ξ replaced by the control u in (27). If the function W is positive definite, then 1 Va (x, ξ) = V(x) + kξ − ℓ(x)k2 2
(32)
is a CLF for the full system (29), (30) (i.e. ξ plays the role of a control in (32)) and there exists a feedback law u = ℓa (x, ξ) that makes the full system (29), (30) globally asymptotically stable around x = 0, ξ = 0. One such control law is u = −c(ξ − ℓ(x)) +
T ∂ℓ(x) ∂V f (x) + g(x)ξ − (x)g(x) , ∂x ∂x
c > 0.
(33)
A proof of this theorem can be found in e.g. Ref 8. In what follows we are going to apply the method of backstepping as outlined in Thm. III.2 to a slightly augmented version of the system (29),(30) namely ˙ x(t) ˙ ξ(t)
= f (x(t)) + g(x(t))ξ(t),
(34)
= h(ξ(t)) + ku(t),
(35)
where ξ takes values in Rm , the function h : Rm → Rm is smooth and k is an invertible (constant) matrix in Rm×m . Since k in (35) is nonsingular, and we can choose the control u freely, the system in (34), (35) can via a simple change of variables be brought to the form (29),(30) and the backstepping control problem for the two systems is therefore one and the same. C.
Standard Form of the Equations of Motion
To be able to apply the theory for integrator backstepping as it stands the system in question has to be on the standard from (29),(30) or (34),(35). However, a glance at (22)–(25) reveals that this system is already on the required standard form. Indeed, if we make the following associations (here ∼ means “corresponds to”) " # " # " # ˜ v ω m(v, v0 , ω, ω 0 ) x∼ , ξ∼ , u∼ (36) Q τ u˜τ and f∼ g∼
"
C(v + v0 ) + Q0 )
1 2 B(Q
"
# ˜ (Q, Q0 ) + v × ω 0 +g , 1 2 A(ω 0 )Q # " # 1 −1 e j (jω × (ω + ω ) + jω × ω) 0 0 m x , h∼ , 04×1 bτ " # j−1 03×1 k∼ 01×3 −b 1˜ m f (v, v0 )
we see that (22)–(25) is on standard form for integrator backstepping.
6 of 20 American Institute of Aeronautics and Astronautics
(37)
(38)
(39)
D.
The Proposed Controller
The first task when developing a backstepping controller for the system on standard form is to find a suitable Lyapunov function V as in Thm. III.2 and an accompanying feedback law ℓ such that the feedback connected first part of the system, corresponding to (34), is stable with suitable dynamics. When determining what is “suitable” dynamics for the system (22), (23) we must take into account at least two obvious requirements; (i) the need to aerodynamically stabilize the aircraft and (ii) the desire to solve the attitude and velocity control problem outlined in Sec. A. It is intuitively clear that these two requirements can not be dealt with independently since rotating the aircraft body so that the body velocity error vector v becomes zero does not necessarily mean that the aircraft has desired orientation Q0 . We are going to solve this problem by combining two controllers, one for the velocity error and one for the orientation error. 1.
Controlling the velocity.
The problem of aerodynamically stabilizing the aircraft, without regard to its orientation, is not too hard (in principle) once the system has been brought onto the form (22)–(25). For instance, one can control ˜ m(v, v0 , ω, ω0 ) such that the velocity vector v + v0 is rotated into a position aligned with v0 , while simultaneously controlling the thrust setting u ˜τ so that the magnitude becomes right. A simple way of achieving a rotation of the velocity vector v+v0 in the right direction is to use a (virtual) control of the form cv ω des v × v0 , (40) v (v, v0 ) = − kv0 k2 where cv is some positive constant. To give some motivation at this point for the choice (40) of a (virtual) control ω des v (v, v0 ) one can note that − (v + v0 ) ×
1 (v × v0 ) = kv0 k2 =
1 1 (v × v0 ) − v0 × (v × v0 ) kv0 k2 kv0 k2 1 −v × (v × v0 ) − P[v0 ]⊥ (v), kv0 k2
−v ×
(41)
where the first term on the right is perpendicular to v and P[v0 ]⊥ (v) is the projection of v onto the subspace of vectors in R3 that are orthogonal to v0 . The first term is in general much smaller in magnitude than the second and therefore, when inserted instead of ω in (22), the (virtual) control ω des v (v, v0 ) in (40) can act to reduce the error v. However, since the main reduction of the velocity error v is in the component of it that is orthogonal to v0 there is a need to complement the control with some action also in the direction of v0 . The direction of v0 is normally almost the same as the direction in which the thrust acts (here, in the body x-direction) and therefore it is natural to try to achieve control action in the v0 direction by (virtual) thrust control. A simple way to achieve this is to employ a negative velocity feedback along the direction of v0 , for example using the (virtual) thrust control τ des (v, v0 , Q, Q0 ) as τ des (v, v0 , Q, Q0 ) = −
cτ mv1 vT v0 2 cτ mv1
P[v ] (v) 2 , =− 0 2 2 cℓ + v1 kv0 k cℓ + v1
(42)
where cτ , cℓ are positive constants and P[v0 ] (v) is the projection of v onto the one dimensional subspace spanned by v0 . This type of virtual thrust control would, if v is large and mostly aligned with v0 , approximately give a stable linear first order contribution to the dynamics in the x-axis component of the force equation (22). (For small v this control would under the same conditions do essentially nothing.) Therefore, at least when when cτ is close to cv it is clear from (40) and (41) that the combined effect of the virtual des controls ω des (v, v0 , Q, Q0 ) for large v is to give the overall system roughly first order stable v (v, v0 ) and τ (virtual) dynamics for the error v. 2.
Controlling the attitude.
We now turn to the problem of controlling the attitude of the aircraft. When conceiving a solution to this problem we will, in analogy with the approach above, neglect the other part of the control problem, viz. the problem of aerodynamically stabilizing the aircraft. The attitude control can be achieved by rotating the body along the shortest path from the current orientation to the desired orientation on the set of
7 of 20 American Institute of Aeronautics and Astronautics
unit norm quaternions, which we here identify with (one “half” of) S 3 , the unit sphere in ordinary four dimensional space. Such a shortest path is the same as a geodesic (in the ordinary Euclidean metric). A simple parametrization for this type of geodesic called slerp (for Spherical Linear Interpolation) was introduced in 1985 by Shoemake.9 The slerp Q that describes the path from unit norm quaternion Q1 to unit norm quaternion Q0 is given by sin(tθ)Q1 sin((1 − t)θ)Q0 + , t ∈ [0, 1], (43) Q(t) = sin(θ) sin(θ) where θ is defined by cos(θ) = QT0 Q1
and the inner product on the right hand side is calculated as for ordinary vectors in four-space. The timederivative of the slerp Q is easily calculated as ˙ Q(t) =
θ (cos(tθ)Q1 − cos((1 − t)θ)Q0 ), sin(θ)
t ∈ [0, 1],
(44)
˙ and this shows that the motion along the slerp path takes place at constant speed, i.e. kQ(t)k ≡ const. In our application, where we want to design a feedback law based on the slerp formula (43), the start quaternion Q1 will be constantly changing and so we really only use the expression (44) for the slerp velocity vector, and evaluate it at the (changing) starting point. Indeed, at least if we assume that Q0 is constant (i.e. ω0 = 0) it is clear that what we want to achieve with the (virtual) attitude control is ˙ = cQ Q(0), ˙ Q (45) where cQ is a positive constant, Q is the state quaternion in (23) and Q + Q0 , Q0 in (23) are used instead of Q0 , Q1 in (44). This gives a condition for the sought virtual angular velocity ω des Q (Q, Q0 ) for attitude control as 1 θ B(Q + Q0 )ωdes Q0 − cos(θ)(Q + Q0 ) , (46) Q (Q, Q0 ) = cQ 2 sin(θ) where now θ is given by
cos(θ) = (Q + Q0 )T Q0 .
(47)
ω des Q (Q, Q0 )
in (46) since the matrix B is not From (46) it might appear impossible to solve (uniquely) for square, but if we remember that the left hand side of (46) is really just another way of writing the product of two quaternions, one unit norm and one pure, we can determine ω des Q (Q, Q0 ) explicitly. Working through the algebra we get 2θ ωdes ℑ(Qc ◦ Q0 ) (48) Q (Q, Q0 ) = cQ sin(θ) (with θ as in (47)) where ℑ( · ) denotes quaternion imaginary part and ( · )c denotes quaternion conjugation. In case Q0 is time-varying (i.e. ω 0 is constant but nonzero) it is clear, after a moments contemplation, that the same principle for selecting ω des Q (Q, Q0 ) ought to apply, and that the resulting dynamics in this “moving” scenario on S 3 then becomes the same as in (23), if we replace ω there by ω des Q (Q, Q0 ). We then have d ˙ +Q ˙ 0 = 1 B(Q + Q0 )(ω des (Q, Q0 ) + ω 0 ), (Q + Q0 ) = Q (49) Q dt 2 ˙ as in (45) and ω des (Q, Q0 ) as in (48). with Q Q
3.
Total controller.
The complete (virtual) controller corresponding to ℓ in Thm. III.2 is now given by the vector ! ! des ω des ωdes (v, v0 , Q, Q0 ) v (v, v0 ) + ω Q (Q, Q0 ) = , τ des (v, v0 , Q, Q0 ) τ des (v, v0 , Q, Q0 ),
(50)
with the components on the right hand side given by (40),(48) and (42). It should be pointed out, however, that many other solutions are possible. Once the model has been put on the standard form as in (22)–(25) there are many ways of constructing stabilizing controllers using e.g. backstepping. 8 of 20 American Institute of Aeronautics and Astronautics
E.
Stability
As a candidate for the “inner” Lyapunov function V as in Thm.III.2 we shall take V(v, Q) =
γv γQ kvk2 + kQk2 , 2 2
(51)
where γv , γQ are two positive constants (to be determined later) and the norms are ordinary 2-norms in R3 and R4 , respectively. With this choice, the time derivative of V along the solutions of (22),(23), and with ω, τ replaced by ω des , τ des as in (50), is given by d V(v, Q) = dt =
=
=
˙ γv vT v˙ + γQ QT Q 1˜ 1 ˜ (Q, Q0 ) + v × ω 0 + (v + v0 ) × ω des + τ des ex f (v, v0 ) + g m m 1 T 1 des +γQ Q A(ω 0 )Q + B(Q + Q0 )ω 2 2 1 T 1 ˜ ˜ (Q, Q0 ) + v0 × ω des + τ des ex γv v f (v, v0 ) + g m m γQ T des + Q B(Q + Q0 )ω 2 γv T ˜ ˜ (Q, Q0 ) v f (v, v0 ) + γv vT g m v1 des T des +γv vT (v0 × ωdes τ v ) + γv v (v0 × ω Q ) + γv m γQ T γ Q T + Q B(Q + Q0 )ω des Q B(Q + Q0 )ω des v + Q , 2 2 γv v T
(52)
where we have used the fact that A(ω 0 ) is skew symmetric. An important part of our main stability result is the following technical assumption about the aerodynamic forces. Assumption III.3. Assume that the function v 7→ vT ˜f (v, v0 ) is locally negative definite, i.e. for some open set U ⊆ R3 containing the origin the following expansion holds ˜f (v, v0 ) = F(v0 )v + O(kvk2 ),
v ∈ U,
(53)
where the symmetric part 21 (F(v0 ) + FT (v0 )) of F(v0 ) is a negative definite matrix, and that the constants cℓ and cv can be selected such that 1 T˜ v12 v f (v, v0 ) − cv kvk2 + cv 1 − kP[v0 ] (v)k2 ≤ −(σcv + c0 )kvk2 , m cℓ + v12
v ∈ U,
(54)
for some constants σ > π/2, c0 > 0.
The condition of negative definiteness of the symmetric part of F(v0 ) is for a large number of aircraft configurations satisfied under normal flying conditions (cf. Sect IV). To see that also the condition (54) is reasonable to assume we let λ3 ≤ λ2 ≤ λ1 < 0 be the eigenvalues of the symmetric part 21 (F(v0 ) + FT (v0 )) of F(v0 ). Then, locally around v = 0 we have 1 T˜ v12 v f (v, v0 ) − cv kvk2 + cv 1 − kP[v0 ] (v)k2 = m cℓ + v12 1 T v12 v F(v0 )v − cv kvk2 + cv 1 − kP[v0 ] (v)k2 + O(kvk3 ) ≤ m cℓ + v12 λ1 v12 kvk2 − cv kvk2 + cv 1 − kP[v0 ] (v)k2 + O(kvk3 ) = m cℓ + v12 λ1 kvk2 − cv (kvk2 − kP[v0 ] (v)k2 ) + O(kvk3 ) ≤ m λ1 kvk2 + O(kvk3 ), m
9 of 20 American Institute of Aeronautics and Astronautics
(55)
where we have used the fact that kvk2 = kP[v0 ] (v)k2 + kP[v0 ]⊥ (v)k2 .
(56)
The first term on the right hand side in (55) is negative definite and if the product σcv satisfies λ1 < −σcv m
(57)
then there exists c0 > 0 such that (54) is fulfilled in some U ⊆ R3 containing the origin. For later reference we also note that (57) gives an upper bound (since λ1 < 0) on the feasible values for the product σcv , which together with the requirement σ > π/2 gives an upper bound on the feasible values for cv . We have now everything in place for our stability result. Lemma III.4. Suppose that Assumption III.3 holds, for some set U and constants σ > π/2, c0 > 0, assume that cτ = cv + c1 for some c1 > 0, and assume that there exists a positive constant θ0 < π/2 such that q γQ cv (58) γv (40 + cQ πkv0 k) + ≤ 2 γv cv γQ cQ σθ0 cot(θ0 ). 2kv0 k
Then the Lyapunov function candidate in (51) is indeed an “inner” Lyapunov function for the integrator backstepping problem for (22)–(25) over the domain of (v, Q) such that v ∈ U and θ in (47) satisfies |θ| < θ0 , i.e. the Lyapunov function candidate satisfies condition (31) in Assumption III.1 over this domain. Proof. For the time derivative of the Lyapunov function candidate (51) along the solutions to the system (22), (23), with the virtual controls (50) inserted, we get using (52) and the bounds (78)–(83) in Appendix the following relations d V(v, Q) dt
γv T ˜ ˜ (Q, Q0 ) + γv vT (v0 × ωdes v f (v, v0 ) + γv vT g v ) m v1 des +γv vT (v0 × ω des τ Q ) + γv m γQ γQ T + QT B(Q + Q0 )ω des Q B(Q + Q0 )ω des v + Q 2 2 γv T ˜ ˆ v, Q)Q − γv cv kP[v ]⊥ (v)k2 ¯ T D(¯ ≤ v f (v, v0 ) + 40 γv v 0 m 2 ˆ v, Q)Q − γv cτ v1 kP[v ] (v)k2 +γv cQ πkv0 k¯ vT D(¯ 0 cℓ + v12 γQ cv T ˆ θ ¯ D(¯ + v v, Q)Q − γQ cQ θ cot(θ)kQk2 − γQ cQ (QT Q0 )2 2kv0 k sin(θ) γv T ˜ v12 v12 2 = v f (v, v0 ) − γv cv k¯ v k 2 + γv c v 1 − kP (v)k − γ c kP[v0 ] (v)k2 v 1 [v ] 0 m cℓ + v12 cℓ + v12 θ −γQ cQ θ cot(θ)kQk2 − γQ cQ (QT Q0 )2 sin(θ) γQ cv T ˆ ¯ D(¯ + γv (40 + cQ πkv0 k) + v v, Q)Q, (59) 2kv0 k
=
ˆ is the orthogonal matrix function in (77) and we have defined the embedded (in R4 ) velocity error where D ¯ v as " # 0 ¯= v , (60) v and in the last equality in (59) we have used the facts (56) and k¯ vk = kvk as well as the assumption cτ = cv + c1 . If we apply condition (54) we can proceed one step further with (59) to obtain d V(v, Q) dt
≤ −γv cv σk¯ vk2 − γv c0 k¯ v k 2 − γv c 1
v12 kP[v0 ] (v)k2 cℓ + v12
θ (QT Q0 )2 sin(θ) γQ cv T ˆ ¯ D(¯ + γv (40 + cQ πkv0 k) + v v, Q)Q. 2kv0 k
−γQ cQ θ cot(θ)kQk2 − γQ cQ
10 of 20 American Institute of Aeronautics and Astronautics
(61)
From condition (58) we get for |θ| < θ0 that γQ cv T ˆ ¯ D(¯ v v, Q)Q ≤ 2kv0 k γQ cv T ˆ ¯ D(¯ −γv cv σk¯ vk2 − γQ cQ θ0 cot(θ0 )kQk2 + γv (40 + cQ πkv0 k) + v v, Q)Q ≤ 2kv0 k q ˆ v, Q)Q = vT D(¯ −γv cv σk¯ vk2 − γQ cQ θ0 cot(θ0 )kQk2 + 2 γv cv γQ cQ σθ0 cot(θ0 )¯ q √ ˆ v, Q)Qk2 , −k γv cv σ¯ v − γQ cQ θ0 cot(θ0 )D(¯
−γv cv σk¯ vk2 − γQ cQ θ cot(θ)kQk2 + γv (40 + cQ πkv0 k) +
(62)
and if we combine (61) with (62) and invoke condition (54) we therefore have d V(v, Q) dt
≤
≤
θ v12 −γv c0 k¯ v k 2 − γQ c Q (QT Q0 )2 − γv c1 kP[v0 ] (v)k2 sin(θ) cℓ + v12 q √ ˆ v, Q)Qk2 ¯ − γQ cQ θ0 cot(θ0 )D(¯ −k γv cv σ v −γv c0 kvk2 − γQ cQ
θ (QT Q0 )2 , sin(θ)
(63)
for v ∈ U, |θ| ≤ |θ0 |. Thus, if we take −W (v, Q) = −γv c0 kvk2 − γQ cQ
θ (QT Q0 )2 , sin(θ)
(64)
and make the identifications (36), (37) we see that we satisfy condition (31) in Assumption III.1 for v ∈ U, |θ| < |θ0 |. Remark III.5. Even though the Assumption III.1 is global in nature it is clear that it can be applied in a local version, as is done here. It may not be immediately clear that the second term in (64) is negative definite in Q but if we recall that for θ 6= 0 we have 1 > cos(θ) = (Q + Q0 )T Q0 = QT Q0 + 1 and thus QT Q0 < 0, and at the same time θ/ sin(θ) > 0 for |θ| < π/2 this should be clear. Having established the above lemma we can now present our main result.
Theorem III.6. Assume that the conditions in Lemma III.4 are fulfilled. Then the integrator backstepping problem for the attitude-velocity control problem for (22)–(25) is solvable using the standard method in Theorem III.2 (using the identifications in (36)–(39)) and one stabilizing control law is given by (33). Proof. The result follows immediately from Lemma III.4 above. 1.
Selection of the constants γv , γQ , cv , cQ
It is clear that only the ratio γv /γQ of the weighting constants in the Lyapunov function V(v, Q) in (51) is relevant for stability and this is also apparent if we rewrite the condition (58) as r γv cv γv (40 + cQ πkv0 k) + ≤2 cv cQ σθ0 cot(θ0 ). (65) γQ 2kv0 k γQ The magnitude of the constants γv , γQ however determines the relative influence of the two terms on the right of (32) which make up the Lyapunov function for the total system. This can also be inferred from (33), where the effect of the magnitude of γv , γQ on the control signal amplitude (i.e. the aerodynamical moments generated by control surface deflections in our setting) enters in the last term on the right. A design procedure can therefore be to first determine the upper bound on the product σcv given by (57) and the desired region of attitude stability given by θ0 in Lemma. III.4. Then one can determine allowable regions for γv /γQ , cv and cQ from (65). Finally, the magnitude of γv , γQ can be determined on the basis of the tradeoffs in (32) and (33). If we introduce the notation k0 = θ0 cot(θ0 ),
(66)
k1 = σcv ,
(67)
11 of 20 American Institute of Aeronautics and Astronautics
the requirement θ0 < π/2 in Lemma III.4 and the bound for the product σcv in (57) can be expressed as 0 < k0 < 1, 0 < k1 < − λm1 ,
(68) (69)
and by using the quantities k0 , k1 we can give a sufficient condition for (65) to be fulfilled in terms of the following set of inequalities d 2 π k0 k1 , v − −2k0 k1 kv40c , 0 k+cv kv0 kπ
0 < cv < cQ >
(70) (71)
and −40cv + cQ kv0 k(4k0 k1 − cv π) √ − 2 2kv0 k(40 + cQ kv0 kπ)2
s
−
40cQcv k0 k1 − 2c2Q k02 k12 kv0 k + c2Q cv k0 k1 kv0 kπ kv0 k(40 + cQ kv0 kπ)4
≤
s 40cQ cv k0 k1 − 2c2Q k02 k12 kv0 k + c2Q cv k0 k1 kv0 kπ γv −40cv + cQ kv0 k(4k0 k1 − cv π) √ − ≤ + . (72) 2 γQ 2kv0 k(40 + cQ kv0 kπ)2 kv0 k(40 + cQ kv0 kπ)4
The relations (68)–(72) together determine an admissible region for cv , cQ , γv and γQ in order to ensure stability of the proposed controller.
IV.
Example
In this section we present some simulation results based on a realistic aircraft model. A.
The ADMIRE Model
The ADMIRE model is a full nonlinear six-degree-of-freedom model of a small single engine fighter aircraft with delta canard configuration (similar to the JAS39 Gripen). It is implemented in Simulink and is freely available from the Swedish Defence Research Agency’s website at http://www.foi.se/admire. It has been used as classroom model at universities and as benchmark model in research collaboration projects between industry and universities. A simplified version of the model, without actuator dynamics and a simpler engine model, has also been implemented in Modelica, and this version is the one used in the present work. A block diagram overview of the Modelica version of the model with controller is shown in Fig. 1. The simulations have been carried out using the Dymola development environment on a personal computer running GNU/Linux. 1.
The controller
The engine time constant b in (11) is set to 0.5 in the model and the upper bound k1 in (69) on the product σcv , as given by (57), is set to 0.0066. This bound varies between 0.0097 at altitude 3000m to 0.0041 at 10000m. We have selected the constant k0 in (68) to 1.05 (this corresponds to an allowable quaternion error angle θ as (47) in of about 60◦ ). The eigenvalues of the symmetric part of F(v0 ) as in (53) is illustrated in Fig. 2. The rest of the constants in the controller are given in Table 1. It should be pointed out, however, that bounds for the various quantities used in the proof of Lemma III.4 are conservative and simulations have shown that it is possible to use considerably different values (higher feedback gains) of the constants with improved performance (in particular with respect to the velocity errors) and maintained stability. B. 1.
Simulations Mixed maneuver flight path.
A nominal flight trajectory lasting for 140 seconds was programmed which included maneuvering as well as straight path flight. The nominal Mach number is held constant at 0.6 throughout the flight. An outline of d These sufficient conditions were derived using Mathematica’s Reduce function. The one situation missing to make a necessary and sufficient condition for (65) is the case where (68)–(70) are fulfilled, and there is equality in (71) and on the left in (72). We also note that the requirement σ > π/2 (rather than just σ > 0) comes from (70).
12 of 20 American Institute of Aeronautics and Astronautics
Backstepping control setup
X
Y
h height... trimout2trim...
pathTable
deMultiplex4
Z
trimBox
r_rel={0,0,0} 2
M aeroDynamicsBasic
referencePoint
M
1
aeroDynInvMom...
toAC
controlBS
r_rel={0,0,0} g = f(h)
2 1
bodyCtrled
earth
atmosphere
toEngine
engineWithDyn
Figure 1. Modelica implementation of a simplified version of the ADMIRE model. The model is viewed using the Dymola development environment. Modelica is an object oriented language for acausal multi-domain modeling. This is reflected in the figure where one can see icons for the various objects representing dynamics, aerodynamics, environment and propulsion. The diamond like icon in the center is a summing point for the variables declared as “flow” variables, such as the total force and moment acting on the body.
the flight path shown in Fig. 3. The trajectory starts with a initial straight path segment at an altitude of 3000m, after which the trajectory continues with a horizontal constant-g turn to the right. About 40s into the flight the trajectory enters a short transition period after which it changes to a constant-g turn right with simultaneous ascent, resulting in a helix motion. This motion is interrupted about 85s into the flight and the trajectory is then changed into straight path flight at an altitude of about 5600m. 2.
Results.
In Figs. 4–7 below various aspects of the performance of the controller are illustrated. As can be seen from Fig. 4 the attitude errors are kept small throughout the flight. The velocity errors are kept reasonably small as well, as shown in Fig. 5. In Fig. 6 the angle of attack and sideslip are displayed. It can be seen that the angle of attack tracks the desired value closely, except perhaps at the points with the highest load factor. The largest deviation occurs when the trajectory changes from upward helix to level flight. As for the sideslip angle, the largest deviations occurs at the same points in the trajectory, and as expected the sideslip angle is not always small. This is due to the fact that the slerp motion, on which the controller is based, does not in general correspond to an ordinary coordinated turn (roll–pitch). In Fig. 7 the (generic) control surface deflections are shown and it can be seen that the sideslip angle closely follows the rudder deflection. The rudder deflection thus primarely acts as to control the attitude and results in the noncoordinated flight behavior. Part of this behavior is mitigated by the aileron which acts as desired and produces rolling motion when desired.
V.
Conclusion
We have demonstrated the power of multi-input backstepping for nonlinear aircraft controller design and derived a nonlinear autopilot for simultaneous attitude and velocity control. The controller incorporates a 13 of 20 American Institute of Aeronautics and Astronautics
Constant c cv cτ cQ γv γQ
Constants Value 1.0 0.0023 0.1 4 2.5 × 10−9 1.0
in the controller Remark Outer loop error gain Eq. (33) Velocity control gain Eq. (40) Thrust control gain Eq. (11) Attitude control gain Eq. (48) Velocity error weight Eq. (51) Attitude error weight Eq. (51)
Table 1. Constants used in the controller.
1000 0 −1000 −2000
Eigenvalues
−3000 −4000 −5000 −6000 −7000 −8000 −9000 −10000 0.2
0.3
0.4
0.5
0.6 Mach number
0.7
0.8
0.9
1
Figure 2. The eigenvalues of the the symmetric part of the force Jacobian F(v0 ) in (53) at 3000m. The largest eigenvalue is approximately −88 and the corresponding eigenvector is almost parallel with ex . As can be seen, the other two eigenvalues are considerably more negative. Thus, there is more natural stability against velocity deviations in the plane [ex ]⊥ and the thrust control helps enhance the stability in the less stable ex direction.
simple and well-known quaternion formula for interpolation between rotations and is capable of controlling the aircraft over a large region of deviations in attitude and moderate deviations in velocity from the nominal values. It relies on only weak assumptions on the aerodynamic forces.
A. A.
Bounds for the terms in the Lyapunov function
General relations
The quaternion norm N ( · ) is defined as the square of the ordinary 2-norm for vectors in R4 and a fundamental property of the quaternion norm is N (Q1 ◦ Q2 ) = N (Q1 )N (Q2 ).
(73)
It therefore follows that if Q1 is a unit norm quaternion we have for the matrix function T in (18) that kT(Q1 )k = sup kT(Q1 )Q2 k = sup kQ1 ◦ Q2 k = sup kQ2 k=1
kQ2 k=1
kQ2 k=1
1/2 N (Q1 )N (Q2 ) = 1,
14 of 20 American Institute of Aeronautics and Astronautics
(74)
Flight trajectory
5000
z [m]
4000
3000
2000
1000
0
2
4
1.5
x 10
1 0.5 x [m]
0
1000
0
2000
3000
4000
5000
6000
y [m]
Figure 3. Flight trajectory.
where k · k denotes both the ordinary 2-norm on R4 and its induced counterpart on R4×4 , respectively. If u, U ∈ R3 are two vectors and R ∈ R3×3 is a rotation matrix such that U = Ru we have (0, u) = Qc ◦ (0, U) ◦ Q,
(75)
where Q is the unit norm quaternion corresponding to R and ( · )c denotes quaternion conjugation (sign change on the imaginary part). In connection with expressions of the form (75) we shall have reason to consider also expressions of the form ℑ(Q1 ◦ (0, U) ◦ Q2 ), (76)
where Q1 ,Q2 are two arbitrary quaternions. The expression (76) will occur in connection with matrix-vector products of the form " # " # 0 0 01×3 = Q1 ◦ (0, U) ◦ Q2 ℑ(Q1 ◦ (0, U) ◦ Q2 ) 03×1 I3×3 " # 0 01×3 = T(Q1 ◦ (0, U))Q2 03×1 I3×3 " # Q1 U 0 01×3 = kQ1 kkUk T( ◦ (0, ))Q2 kQ1 k kUk 03×1 I3×3 =
kQ1 kkUk S(Q1 , U)Q2 ,
where the matrix S(Q1 , U) is given by S(Q1 , U) =
"
0 01×3 03×1 I3×3
#
T(
Q1 U ◦ (0, )) kQ1 k kUk
From (73) and (74) it follows that we have a bound for the matrix S(Q1 , U) as kS(Q1 , U)k ≤ 1.
Finally we shall make an elementary observation about maximization of bilinear forms. Let x, y ∈ R4 be two arbitrary vectors. By the Cauchy-Schwarz inequality we then have sup D∈R4×4 ,kDk≤1
yT Dx ≤ kxkkyk.
15 of 20 American Institute of Aeronautics and Astronautics
0.5 60
50
Euler angles [degrees]
Quaternion components Reference quaternion components
Yaw angle Pitch angle Roll angle Reference Euler angles
0.4
40
0.3
30
0.2
20
0.1
10
0 0
−0.1 −10
−0.2 −20
−0.3
−30
−40
0
20
40
60
80
100
120
140
−0.4
0
20
40
60
80
100
120
140
Time [s]
Time [s]
(a) Euler angles
(b) Quaternion values
Figure 4. Euler angles and three components (of four) of the quaternion. Good tracking behavior is seen for all components of the Euler angles and quaternion.
ˆ 1 given by It is easy to see that the supremum is attained for the rank-1 matrix D ˆ1 = D
1 yxT . kxkkyk
Thus, if we collect three vectors u1 , u2 , u3 ∈ R4 that are orthonormal and orthogonal to x and likewise collect three more vectors v1 , v2 , v3 ∈ R4 that are orthonormal and orthogonal to y, and form the matrix ˆ D(x, y) given by xT kxk
ˆ D(x, y) =
h
we have sup D∈R4×4 ,kDk≤1
y kyk
v1
v2
i uT v3 1 uT 2 uT3
(77)
ˆ yT Dx = yT D(x, y)x = kxkkyk
ˆ where the maximizer D(x, y) is an orthogonal matrix. B.
Bounds for specific terms
Using the above results we can now establish a series of bounds for the various terms on the right hand side of (52).
16 of 20 American Institute of Aeronautics and Astronautics
15
0.616 Mach number Reference Mach number
0.614
10 0.612
0.61
Mach number [−]
Velocity deviation [m/s]
5
0
0.608
0.606
0.604
−5 0.602
0.6
−10
Velocity deviation in body x Velocity deviation in body y
0.598
Velocity deviation in body z −15
0
20
40
60
80
100
120
140
0.596
0
20
40
60
Time [s]
80
100
120
140
Time [s]
(a) Velocity error components
(b) Mach number
Figure 5. The three velocity error components and the Mach number. The controller tracks the desired velocity components and Mach number well.
1.
Gravity term.
The second term on the right in (52) can, using (75), be bounded as ˜ (Q, Q0 ) = γv vT g = = = ≤ ≤ =
¯ T (Q + Q0 )c ◦ (0, G) ◦ (Q + Q0 ) − Qc0 ◦ (0, G) ◦ Q0 γv v
¯ T Qc ◦ (0, G) ◦ Q + Qc0 ◦ (0, G) ◦ Q − (Qc0 ◦ (0, G) ◦ Q)c γv v ¯ T Qc ◦ (0, G) ◦ Q + 2(0, ℑ(Qc0 ◦ (0, G) ◦ Q)) γv v Qc G G γv kQkkGk¯ vT T ◦ (0, ) Q + 2γv kGk¯ vT S Qc0 , Q kQk kGk kGk Qc G G 20 γv |¯ vT T ◦ (0, ) Q| + 20 γv |¯ vT S Qc0 , Q| kQk kGk kGk 40 γv
¯ T DQ v
sup
D∈R4×4 ,kDk=1
ˆ v, Q)Q, ¯ T D(¯ 40 γv v
(78)
ˆ v, Q) is the matrix function in (77) and we have used the facts that kQ0 k = 1, kQk ≤ 2 and where D(¯ kGk ≤ 10. 2.
Desired velocity control term.
For the term involving v and the virtual control ω des we have similarly v γv vT (v0 × ω des v )
= −γv cv vT
v0 v0 × (v × ) kv0 k kv0 k
= −γv cv vT P[v0 ]⊥ (v)
= −γv cv kP[v0 ]⊥ (v)k2 .
17 of 20 American Institute of Aeronautics and Astronautics
(79)
4
4.5 AoA Reference AoA
3
4
2
Side slip angle [degrees]
Angle of attack [degrees]
3.5
3
2.5
1
0
−1
2
−2
1.5
1
−3
0
20
40
60
80
100
120
140
−4
0
20
40
60
80
100
120
140
Time [s]
Time [s]
(a) Angle of attack
(b) Side slip angle
Figure 6. Angle of attack and sideslip angle. As expected the angle of attack can differ from the desired values at points of high load factor but apart from these regions the tracking is very good. The sideslip angle deviates in the beginning of each change in turn direction since the controller is not designed to enforce coordinated maneuvering.
3.
Undesired velocity control term.
The term in (52) involving both v and ω des Q can be bounded as γv vT (v0 × ωdes Q ) = = = =
γv vT v0 ×
2cQ θ ℑ(Qc ◦ Q0 ) sin(θ)
2θ vT C(v0 )ℑ(Qc0 ◦ Q) sin(θ) " # 2θ 0 01×3 T ¯ T(Qc0 )Q −γv cQ v sin(θ) 03×1 C(v0 ) # " 0 01×3 2θ T −γv cQ T(Qc0 )Q kv0 k¯ v sin(θ) 03×1 C( kvv00 k ) −γv cQ
≤
γv cQ πkv0 k
sup
=
ˆ v, Q)Q, γv cQ πkv0 k¯ vT D(¯
¯ T DQ v
D∈R4×4 ,kDk=1
(80)
where we have used the easily verified fact that kC(v0 )k = kv0 k, the bound (74) for T(Qc0 ) and the bound 1 ≤ θ/ sin(θ) ≤ π/2 for |θ| ≤ π/2. 4.
Virtual thrust term.
The virtual thrust control term in (52) is simply γv
v1 des v12 τ = −γv cτ kP[v0 ] (v)k2 . m cℓ + v12
18 of 20 American Institute of Aeronautics and Astronautics
(81)
−3
1
x 10
5
Lyapunov function Aileron Elevator Rudder
Control surface deflection [degrees]
4
0
3
−1
2
−2
1
−3
0
−4 −1
−5 −2
−6 −3
−7
−4
−5
−8 0
20
40
60
80
100
120
140
0
20
40
60
80
100
120
140
Time [s]
Time [s]
(a) Control surface deflections
(b) Time derivative of Lyapunov function
Figure 7. Control surface deflections and time derivative of the Lyapunov function (51) along the solution trajectory for the system (22)–(25) with controller added. The control surface deflections are very moderate.
5.
Undesired attitude control term.
For the term involving Q and ωdes we start by recalling that by the definition (19) of B we have v B(Q + Q0 )(v × v0 ) = = = =
A(v × v0 )(Q + Q0 ) −A(C(v0 )v)(Q + Q0 )
−(Q + Q0 ) ◦ (0, C(v0 )v) " 0 −kv0 kT(Q + Q0 ) 03×1
01×3 C( kvv00 k )
#
¯. v
With this we obtain the following bound γQ T γQ cv T Q B(Q + Q0 )ω des = − Q B(Q + Q0 )(v × v0 ) v 2 2kv0 k2 " # 0 01×3 γQ cv T ¯ Q T(Q + Q0 ) = v 2kv0 k 03×1 C( kvv00 k ) " # 0 01×3 γQ cv T ¯ = − v TT (Q + Q0 )Q 2kv0 k 03×1 C( kvv00 k ) γQ cv ¯ T DQ sup v ≤ 2kv0 k D∈R4×4 ,kDk=1 γQ cv T ˆ ¯ D(¯ = v v, Q)Q. 2kv0 k 6.
(82)
Desired attitude control term.
From the defining relation (46) for ωdes Q we have γQ T Q B(Q + Q0 )ω des Q 2
θ QT Q0 − cos(θ)(Q + Q0 ) sin(θ) θ = −γQ cQ θ cot(θ)kQk2 + γQ cQ (1 − cos(θ))QT Q0 sin(θ) θ = −γQ cQ θ cot(θ)kQk2 − γQ cQ (QT Q0 )2 , sin(θ)
= γQ c Q
19 of 20 American Institute of Aeronautics and Astronautics
(83)
where we have used the fact that (47) implies 1 − cos(θ) = −QT Q0 for the last equality.
References 1 Lane, S. H. and Stengel, R. F., “Flight Control Design Using Non-linear Inverse Dynamics,” Automatica, Vol. 24, No. 4, 1988, pp. 471–483. 2 Ito, D., Ward, D., and Valasek, J., “Robust Dynamic Inversion Controller Design and Analysis for the X-38,” Proc. AIAA Cuidance Control and Navigation Conference and Exhibit, AIAA, Montreal, Canada, 2001. 3 Smith, P. R., “A Simplified Approach to Nonlinear Dynamic Inversion Based Flight Control,” AIAA Atmospheric Flight Mechanics Conference and Exhibit, Boston, MA, 1998. 4 Sharma, M. and Ward, D., “Flight-Path Angle Control via Neuro-Adaptive Backstepping,” Proc. American Control Conf., Anchorage, AK, 2002. 5 Steinicke, A. and Michalka, G., “Improving Transient Performance of Dynamic Inversion Missile Autopilot By Use Of Backstepping,” Proc. AIAA Cuidance Control and Navigation Conference and Exhibit, AIAA, Monterey, CA, 2002. 6 Glad, T. and H¨ arkeg˚ ard, O., “Backstepping Control of a Rigid Body,” Proc. CDC ’02 , Las Vegas, NV, 2002. 7 Stevens, B. and Lewis, F., Aircraft Control and Simulation, Wiley, Hoboken NJ, 2003. 8 Khalil, H., Nonlinear Systems, 2:nd Ed, Prentice-Hall, Upper Saddle River, NJ, 2002. 9 Shoemake, K., “Animating Rotation with Quaternion Curves,” Proceedings of SIGGRAPH ’85 , ACM, San Francisco, CA, 1985, pp. 245–254.
20 of 20 American Institute of Aeronautics and Astronautics
