# A Unified Approach to Nonlinear Dynamic Inversion Control with Parameter Determination by Eigenvalue Assignment

Authors: Yu-Chi Wang; Donglong Sheu; Chin-E Lin
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This paper presents a unified approach to nonlinear dynamic inversion control algorithm with the parameters for desired dynamics determined by using an eigenvalue assignment method, which may be applied in a very straightforward and convenient way. By using this method, it is not necessary to transform the nonlinear equations into linear equations by feedback linearization before beginning control designs. The applications of this method are not limited to affine nonlinear control systems or limited to minimum phase problems if the eigenvalues of error dynamics are carefully assigned so that the desired dynamics is stable. The control design by using this method is shown to be robust to modeling uncertainties. To validate the theory, the design of a UAV control system is presented as an example. Numerical simulations show the performance of the design to be quite remarkable.


---

## Body

## 1. Introduction

In the development of high-performance aircraft, control difficulties may be encountered over some parts of flight envelope. These difficulties arise from highly nonlinear aerodynamic properties [1] in some flight conditions. In order to solve these control difficulties, nonlinear controllers are required for high-performance aircraft.Among many control methods, nonlinear dynamic inversion (NDI) is very popular and has been widely studied for flight control designs (e.g., [2, 3]). NDI-based control systems are usually divided into fast- and slow-loop control subsystems according to the multiple time-scale method [4]. In each subsystem, Lie derivatives [5] are used to transform the nonlinear equations into linear equations. Then, linear control design methods can be employed and the control inputs are obtained by converting the linear system control variables into the original coordinates. However, the control systems obtained by feedback linearization [6] may have nonminimum phase problems for affine or nonaffine nonlinear system [7] and robust issues in case of model mismatch. A typical nonminimum phase problem may be found in flight dynamics where the altitude-elevator transfer function usually has a right-half zero. The internal state control [8] is often used to overcome these nonminimum phase problems. In addition, the fuzzy logic control [9] was also applied for solving these kinds of problems. Furthermore, to overcome the robust problems, μ-analysis [10] and H
∞ method [11] were applied. Specially, incremental NDI (INDI) [12] was used to increase the robustness to aerodynamic uncertainties by calculating the control surface deflection changes instead of giving inputs directly.To circumvent some aforementioned robust problems, an adaptive nonlinear model inversion control [13] was introduced, in which the design concept is similar to the conventional NDI yet without linearizing the nonlinear system. The model inversion method replaces the motion rates with a P-form or PI-form desired dynamics to negate the original dynamics. The choice of parameters in the desired dynamics is based on the bandwidth of response and time scales. The effects of different types of desired dynamics on the resulted control system were discussed [14].Although the aforementioned NDI approaches are successful in many flight control system designs [14–16] over a large part of the flight envelope, the systems of dynamics in general have to be separated into several subsystems according to the rates of response. There are many cases in which the fast rate and the slow rate might not be distinguished so clearly however. Also, although the pole assignment method had been introduced to determine the parameters in the desired dynamics, very often the control system must first be transformed into a standard feedback control form from which the standard eigenvalue assignment method can be applied.In consideration of the aforementioned problems existing in the current literature on control designs, a unified approach to nonlinear dynamic inversion control is proposed in this paper. The equations of motion will not be necessary to be separated into fast rate and slow rate groups, nor will they be limited to an affine system. Feedback linearizations will not be required to transform the nonlinear equations into linear equations. Nonminimum phase problems are solved by eigenvalue assignments for error dynamics. An iterative method for determining the parameters of the desired dynamics from the assigned eigenvalues of error dynamics is proposed. Analysis of robustness to model uncertainties or disturbances is conducted. This method will be ready for design without simplifying the system of equations based on physical insights once the governing dynamic equations are established and the state variables to be tracked are selected. The theory is to be developed in detail in the following sections. A UAV is introduced and its control system is designed with the developed method. Numerical simulations are conducted to validate this method.


## 2. A Unified Approach to Nonlinear Dynamic Inversion Control

### 2.1. Nonlinear Dynamic Inversion Control

In general, dynamic equations of motion with control inputs can be expressed by(1)
x
˙
=
f
x
,
u
,where x
∈
R
n is the state vector, u
∈
R
m (m
<
n) the control vector, and f the nonlinear function representing the model of dynamics with controls. By extending the concept of dynamic inversion (DI) [4], the control vector u can be assumed to be computed from(2)
x
˙
d
=
f
x
d
,
u
,where x
d
∈
R
n is a desired state vector with its rate of change being designated.In this paper, the desired dynamics is designated as a set of stable first-order differential equations:(3)
x
˙
d
=
Ω
x
d
-
x
,where Ω
∈
R
n
×
n represents a constant matrix with n independent parameters which can be chosen. Substituting (3) into (2) yields(4)
Ω
x
d
-
x
=
f
x
d
,
uwhich constitutes a set of n algebraic equations. Since m
<
n, obviously, u cannot satisfy (4) if all elements of x
d are to be designated. It means that only part of x
d can be designated. So let x
d be divided into two groups, say x
c
∈
R
m and x
r
∈
R
n
-
m, where x
c contains the state variables which are to be controlled or designated and x
r the residual ones. Both u and x
r constitute n unknown variables which are to be determined from (4). To solve a set of nonlinear algebraic equations, the Newton-Raphson iteration method can be employed as follows:(5)
u
x
r
k
+
1
=
u
x
r
k
-
∂
f
e
r
r
∂
u
∂
f
e
r
r
∂
x
r
-
1
f
e
r
r
,where f
e
r
r
≜
f
x
d
,
u
-
Ω
x
d
-
x and k
=
0,1
,
2
,
… represents the iteration number.


### 2.2. Parameter Determination

Now, a question arises whether the state vectorx will asymptotically follow the desired vector x
d if u is determined from (4) and substituted in (1). In order to answer it, let (1) and (4) be examined more carefully as follows.Subtracting (4) from (1) yields(6)
x
˙
-
Ω
x
d
-
x
=
f
x
,
u
-
f
x
d
,
u
.If x is very close to x
d, then the above equation can be linearized as follows:(7)
x
˙
-
Ω
x
d
-
x
≈
f
x
x
d
,
u
x
-
x
d
.With the concept that the desired variables x
d are near constant, say x
˙
d
≈
0, (7) can be rewritten as(8)
x
˙
-
x
˙
d
≈
-
Ω
+
f
x
x
d
,
u
x
-
x
d
.Defining an error vector e
=
x
-
x
d and replacing the approximate sign with the equal sign lead (8) to(9)
e
˙
=
-
Ω
+
f
x
x
d
,
u
e
.Equation (9) is a set of error dynamics in which the error vector e will vanish eventually if all the real parts of the eigenvalues of -
Ω
+
f
x
x
d
,
u are negative. This is possible if Ω is chosen appropriately. It means that the state vector x can approach the desired vector x
d once e approaches 0. Recall that x
d contains x
c, a state vector to be tracked.Notice that the eigenvalueλ in (9) can be determined by(10)
λ
I
+
Ω
-
f
x
x
d
,
u
=
0
,where I is an identity matrix. For simplicity, if Ω is diagonal and all elements are the same, say, Ω
=
σ
I, then (10) can be rewritten as(11)
λ
+
σ
I
-
f
x
x
d
,
u
=
0
.In order to make e bounded, σ can be so chosen that(12)
real
⁡
λ
i
+
σ
<
0for all λ
i
(
i
=
1,2
,
…
,
n
). Although this method is simple, the resulting σ may be unnecessarily large.For more general cases,Ω contains a set of parameters, σ
k
,
(
k
=
1,2
,
3
,
…
,
n
), which may be chosen to determine the eigenvalues and eigenvectors of -
Ω
+
f
x
x
d
,
u. Recall that (3) must be a stable model and, therefore, the simplest way of constructing Ω is to let all its elements vanish except those at the diagonal, which are assumed to be Ω
k
k
=
σ
k
>
0. In fact, some off-diagonal elements can also be allowed to exist. For example, let Ω
k
k
=
Ω
l
l
=
σ
k
>
0, and Ω
k
l
2
=
Ω
l
k
2
=
σ
l
2
<
σ
k
2. It is trivial to prove that the latter case is also a stable model.Now, assume that theith (
i
=
1,2
,
3
,
…
,
n
) eigenvalue and eigenvector are λ
-
i and e
-
i, respectively. Accordingly,(13)
-
Ω
-
+
f
x
x
d
,
u
e
-
i
=
λ
-
i
e
-
i
.In general, the eigenvector e
-
i can be normalized so that(14)
e
-
i
T
e
-
i
=
1
.However, if the eigenvalue λ
i is desired rather than λ
-
i, then it is necessary to adjust the parameters in Ω. Assume that an increment ∂
Ω
/
∂
σ
k
Δ
σ
k to Ω
- is required. Then(15)
-
Ω
-
-
∂
Ω
∂
σ
k
Δ
σ
k
+
f
x
x
d
,
u
e
i
=
λ
i
e
i
,where λ
i and e
i are assumed to be approximated to λ
-
i
+
∂
λ
i
/
∂
σ
k
Δ
σ
k and e
-
i
+
∂
e
i
/
∂
σ
k
Δ
σ
k, respectively, and e
i is also normalized so that(16)
e
-
i
+
∂
e
i
∂
σ
k
Δ
σ
k
T
e
-
i
+
∂
e
i
∂
σ
k
Δ
σ
k
=
1
.Accordingly, the derivative of (13) with respect to σ
k results in(17)
-
∂
Ω
∂
σ
k
e
-
i
+
-
Ω
-
+
f
x
x
d
,
u
∂
e
i
∂
σ
k
=
∂
λ
i
∂
σ
k
e
-
i
+
λ
-
i
∂
e
i
∂
σ
kand the derivative of (14) becomes(18)
e
-
i
T
∂
e
i
∂
σ
k
=
0
.Equations (17) and (18) can be rearranged to(19)
e
-
i
λ
-
i
I
+
Ω
-
-
f
x
x
d
,
u
0
e
-
i
T
∂
λ
i
∂
σ
k
∂
e
i
∂
σ
k
=
-
∂
Ω
∂
σ
k
e
-
i
0from which ∂
λ
i
/
∂
σ
k can be determined. Note that ∂
λ
i
/
∂
σ
k represents the variation of λ
i with respect to σ
k. With all ∂
λ
i
/
∂
σ
k, (
k
=
1,2
,
3
,
…
,
n
), the ith revised eigenvalue should be(20)
λ
i
=
λ
-
i
+
∑
∂
λ
i
∂
σ
k
Δ
σ
k
,
i
=
1,2
,
3
,
…
,
n
.Recall that, in this equation, λ
-
i represents the ith eigenvalue of -
Ω
-
+
f
x
x
d
,
u, where Ω
- contains the parameters σ
-
k, (k
=
1,2
,
3
,
…
,
n). Now since λ
i is the ith desired eigenvalue, (20) in fact becomes a set of n simultaneous equations from which Δ
σ
k can be solved. With σ
-
k being the initial guess, revisions for parameters can be made by(21)
σ
k
=
σ
-
k
+
Δ
σ
k
.Since λ
i
,
(
i
=
1,2
,
3
,
…
,
n
) are nonlinear function of σ
k
,
(
k
=
1,2
,
3
,
…
,
n
), iterative computations of (9) and (19)–(21) are required in order to obtain a set of convergent σ
k. To make it clear, the iteration procedures are summarized as follows.(1)
Let the desired eigenvalues beλ
i
(
i
=
1,2
,
3
,
…
,
n
) which are all distinct.
(2)
Guess a set of parametersσ
-
k
(
k
=
1,2
,
3
,
…
,
n
).
(3)
Determine the eigenvalues and eigenvectors of (9). Denote the ith eigenvalue and eigenvector as λ
-
i and e
-
i, respectively.
(4)
With eachσ
-
k
(
k
=
1,2
,
3
,
…
,
n
), determine ∂
λ
i
/
∂
σ
k
(
i
=
1,2
,
3
,
…
,
n
) from (19).
(5)
DetermineΔ
σ
k
(
k
=
1,2
,
3
,
…
,
n
) from (20). If Δ
σ
k is less than a preset small value, then stop; else continue with the next step.
(6)
Determineσ
k from (21).
(7)
Replaceσ
-
k with σ
k and go to step 3.At this point, it must be mentioned that the iterations will converge only if the initial guesses are very close to the true answers. It is also emphasized that the desired eigenvalues λ
i must be chosen to lie on the left-half s-plane and the resulting parameters in Ω must satisfy the stable model requirements in (3).


### 2.3. Robust Analysis

The nonlinear dynamic inversion control design developed so far is based on a nominal dynamical model. In the real world, however, there always exist some modeling uncertainties which cannot be determined in advance. In order to check if the control design is robust, let the actual dynamical model be represented by(22)
x
˙
=
f
′
x
,
u
.With this assumption, now (6) can be modified to(23)
x
˙
-
Ω
x
d
-
x
=
f
′
x
,
u
-
f
x
d
,
uwhich can be rewritten as follows:(24)
x
˙
+
Ω
x
=
Ω
x
d
+
f
x
x
d
,
u
x
-
x
d
+
Δ
f
n
x
,
x
d
,
u
+
f
′
x
,
u
-
f
x
,
u
,where Δ
f
n
x
,
x
d
,
u denotes the nonlinear part of f
x
,
u
-
f
x
d
,
u. The solution of the above equation can be represented as follows:(25)
x
t
=
e
-
Ω
t
x
0
+
I
-
e
-
Ω
t
x
d
+
∫
0
t
e
-
Ω
t
-
τ
f
x
x
d
,
u
τ
x
τ
-
x
d
+
Δ
f
n
x
τ
,
x
d
,
u
τ
+
f
′
x
τ
,
u
τ
-
f
x
τ
,
u
τ
d
τ
.Note that from (9), the eigenvalues of -
Ω
+
f
x
x
d
,
u are all in the left-half s-plane since Ω has been determined with that assumption. Also, at this point, it is not unreasonable to assume that both Δ
f
n
x
τ
,
x
d
,
u
τ and the modeling differences f
′
x
,
u
-
f
x
,
u are bounded. Therefore, if the model in (3) is stable enough, the integral part in (25) will vanish along with e
-
Ω
t. Accordingly, as the time t gets large enough, the state variables x will approach the desired values x
d as can be found from (25) even if there are some modeling uncertainties or disturbances.


## 2.1. Nonlinear Dynamic Inversion Control

In general, dynamic equations of motion with control inputs can be expressed by(1)
x
˙
=
f
x
,
u
,where x
∈
R
n is the state vector, u
∈
R
m (m
<
n) the control vector, and f the nonlinear function representing the model of dynamics with controls. By extending the concept of dynamic inversion (DI) [4], the control vector u can be assumed to be computed from(2)
x
˙
d
=
f
x
d
,
u
,where x
d
∈
R
n is a desired state vector with its rate of change being designated.In this paper, the desired dynamics is designated as a set of stable first-order differential equations:(3)
x
˙
d
=
Ω
x
d
-
x
,where Ω
∈
R
n
×
n represents a constant matrix with n independent parameters which can be chosen. Substituting (3) into (2) yields(4)
Ω
x
d
-
x
=
f
x
d
,
uwhich constitutes a set of n algebraic equations. Since m
<
n, obviously, u cannot satisfy (4) if all elements of x
d are to be designated. It means that only part of x
d can be designated. So let x
d be divided into two groups, say x
c
∈
R
m and x
r
∈
R
n
-
m, where x
c contains the state variables which are to be controlled or designated and x
r the residual ones. Both u and x
r constitute n unknown variables which are to be determined from (4). To solve a set of nonlinear algebraic equations, the Newton-Raphson iteration method can be employed as follows:(5)
u
x
r
k
+
1
=
u
x
r
k
-
∂
f
e
r
r
∂
u
∂
f
e
r
r
∂
x
r
-
1
f
e
r
r
,where f
e
r
r
≜
f
x
d
,
u
-
Ω
x
d
-
x and k
=
0,1
,
2
,
… represents the iteration number.


## 2.2. Parameter Determination

Now, a question arises whether the state vectorx will asymptotically follow the desired vector x
d if u is determined from (4) and substituted in (1). In order to answer it, let (1) and (4) be examined more carefully as follows.Subtracting (4) from (1) yields(6)
x
˙
-
Ω
x
d
-
x
=
f
x
,
u
-
f
x
d
,
u
.If x is very close to x
d, then the above equation can be linearized as follows:(7)
x
˙
-
Ω
x
d
-
x
≈
f
x
x
d
,
u
x
-
x
d
.With the concept that the desired variables x
d are near constant, say x
˙
d
≈
0, (7) can be rewritten as(8)
x
˙
-
x
˙
d
≈
-
Ω
+
f
x
x
d
,
u
x
-
x
d
.Defining an error vector e
=
x
-
x
d and replacing the approximate sign with the equal sign lead (8) to(9)
e
˙
=
-
Ω
+
f
x
x
d
,
u
e
.Equation (9) is a set of error dynamics in which the error vector e will vanish eventually if all the real parts of the eigenvalues of -
Ω
+
f
x
x
d
,
u are negative. This is possible if Ω is chosen appropriately. It means that the state vector x can approach the desired vector x
d once e approaches 0. Recall that x
d contains x
c, a state vector to be tracked.Notice that the eigenvalueλ in (9) can be determined by(10)
λ
I
+
Ω
-
f
x
x
d
,
u
=
0
,where I is an identity matrix. For simplicity, if Ω is diagonal and all elements are the same, say, Ω
=
σ
I, then (10) can be rewritten as(11)
λ
+
σ
I
-
f
x
x
d
,
u
=
0
.In order to make e bounded, σ can be so chosen that(12)
real
⁡
λ
i
+
σ
<
0for all λ
i
(
i
=
1,2
,
…
,
n
). Although this method is simple, the resulting σ may be unnecessarily large.For more general cases,Ω contains a set of parameters, σ
k
,
(
k
=
1,2
,
3
,
…
,
n
), which may be chosen to determine the eigenvalues and eigenvectors of -
Ω
+
f
x
x
d
,
u. Recall that (3) must be a stable model and, therefore, the simplest way of constructing Ω is to let all its elements vanish except those at the diagonal, which are assumed to be Ω
k
k
=
σ
k
>
0. In fact, some off-diagonal elements can also be allowed to exist. For example, let Ω
k
k
=
Ω
l
l
=
σ
k
>
0, and Ω
k
l
2
=
Ω
l
k
2
=
σ
l
2
<
σ
k
2. It is trivial to prove that the latter case is also a stable model.Now, assume that theith (
i
=
1,2
,
3
,
…
,
n
) eigenvalue and eigenvector are λ
-
i and e
-
i, respectively. Accordingly,(13)
-
Ω
-
+
f
x
x
d
,
u
e
-
i
=
λ
-
i
e
-
i
.In general, the eigenvector e
-
i can be normalized so that(14)
e
-
i
T
e
-
i
=
1
.However, if the eigenvalue λ
i is desired rather than λ
-
i, then it is necessary to adjust the parameters in Ω. Assume that an increment ∂
Ω
/
∂
σ
k
Δ
σ
k to Ω
- is required. Then(15)
-
Ω
-
-
∂
Ω
∂
σ
k
Δ
σ
k
+
f
x
x
d
,
u
e
i
=
λ
i
e
i
,where λ
i and e
i are assumed to be approximated to λ
-
i
+
∂
λ
i
/
∂
σ
k
Δ
σ
k and e
-
i
+
∂
e
i
/
∂
σ
k
Δ
σ
k, respectively, and e
i is also normalized so that(16)
e
-
i
+
∂
e
i
∂
σ
k
Δ
σ
k
T
e
-
i
+
∂
e
i
∂
σ
k
Δ
σ
k
=
1
.Accordingly, the derivative of (13) with respect to σ
k results in(17)
-
∂
Ω
∂
σ
k
e
-
i
+
-
Ω
-
+
f
x
x
d
,
u
∂
e
i
∂
σ
k
=
∂
λ
i
∂
σ
k
e
-
i
+
λ
-
i
∂
e
i
∂
σ
kand the derivative of (14) becomes(18)
e
-
i
T
∂
e
i
∂
σ
k
=
0
.Equations (17) and (18) can be rearranged to(19)
e
-
i
λ
-
i
I
+
Ω
-
-
f
x
x
d
,
u
0
e
-
i
T
∂
λ
i
∂
σ
k
∂
e
i
∂
σ
k
=
-
∂
Ω
∂
σ
k
e
-
i
0from which ∂
λ
i
/
∂
σ
k can be determined. Note that ∂
λ
i
/
∂
σ
k represents the variation of λ
i with respect to σ
k. With all ∂
λ
i
/
∂
σ
k, (
k
=
1,2
,
3
,
…
,
n
), the ith revised eigenvalue should be(20)
λ
i
=
λ
-
i
+
∑
∂
λ
i
∂
σ
k
Δ
σ
k
,
i
=
1,2
,
3
,
…
,
n
.Recall that, in this equation, λ
-
i represents the ith eigenvalue of -
Ω
-
+
f
x
x
d
,
u, where Ω
- contains the parameters σ
-
k, (k
=
1,2
,
3
,
…
,
n). Now since λ
i is the ith desired eigenvalue, (20) in fact becomes a set of n simultaneous equations from which Δ
σ
k can be solved. With σ
-
k being the initial guess, revisions for parameters can be made by(21)
σ
k
=
σ
-
k
+
Δ
σ
k
.Since λ
i
,
(
i
=
1,2
,
3
,
…
,
n
) are nonlinear function of σ
k
,
(
k
=
1,2
,
3
,
…
,
n
), iterative computations of (9) and (19)–(21) are required in order to obtain a set of convergent σ
k. To make it clear, the iteration procedures are summarized as follows.(1)
Let the desired eigenvalues beλ
i
(
i
=
1,2
,
3
,
…
,
n
) which are all distinct.
(2)
Guess a set of parametersσ
-
k
(
k
=
1,2
,
3
,
…
,
n
).
(3)
Determine the eigenvalues and eigenvectors of (9). Denote the ith eigenvalue and eigenvector as λ
-
i and e
-
i, respectively.
(4)
With eachσ
-
k
(
k
=
1,2
,
3
,
…
,
n
), determine ∂
λ
i
/
∂
σ
k
(
i
=
1,2
,
3
,
…
,
n
) from (19).
(5)
DetermineΔ
σ
k
(
k
=
1,2
,
3
,
…
,
n
) from (20). If Δ
σ
k is less than a preset small value, then stop; else continue with the next step.
(6)
Determineσ
k from (21).
(7)
Replaceσ
-
k with σ
k and go to step 3.At this point, it must be mentioned that the iterations will converge only if the initial guesses are very close to the true answers. It is also emphasized that the desired eigenvalues λ
i must be chosen to lie on the left-half s-plane and the resulting parameters in Ω must satisfy the stable model requirements in (3).


## 2.3. Robust Analysis

The nonlinear dynamic inversion control design developed so far is based on a nominal dynamical model. In the real world, however, there always exist some modeling uncertainties which cannot be determined in advance. In order to check if the control design is robust, let the actual dynamical model be represented by(22)
x
˙
=
f
′
x
,
u
.With this assumption, now (6) can be modified to(23)
x
˙
-
Ω
x
d
-
x
=
f
′
x
,
u
-
f
x
d
,
uwhich can be rewritten as follows:(24)
x
˙
+
Ω
x
=
Ω
x
d
+
f
x
x
d
,
u
x
-
x
d
+
Δ
f
n
x
,
x
d
,
u
+
f
′
x
,
u
-
f
x
,
u
,where Δ
f
n
x
,
x
d
,
u denotes the nonlinear part of f
x
,
u
-
f
x
d
,
u. The solution of the above equation can be represented as follows:(25)
x
t
=
e
-
Ω
t
x
0
+
I
-
e
-
Ω
t
x
d
+
∫
0
t
e
-
Ω
t
-
τ
f
x
x
d
,
u
τ
x
τ
-
x
d
+
Δ
f
n
x
τ
,
x
d
,
u
τ
+
f
′
x
τ
,
u
τ
-
f
x
τ
,
u
τ
d
τ
.Note that from (9), the eigenvalues of -
Ω
+
f
x
x
d
,
u are all in the left-half s-plane since Ω has been determined with that assumption. Also, at this point, it is not unreasonable to assume that both Δ
f
n
x
τ
,
x
d
,
u
τ and the modeling differences f
′
x
,
u
-
f
x
,
u are bounded. Therefore, if the model in (3) is stable enough, the integral part in (25) will vanish along with e
-
Ω
t. Accordingly, as the time t gets large enough, the state variables x will approach the desired values x
d as can be found from (25) even if there are some modeling uncertainties or disturbances.


## 3. Nonlinear Dynamic Inversion Flight Control System Design for a UAV: An Example

### 3.1. Flight Dynamics Equations of Motion

To illustrate the theory, a design of flight control system with the method developed is presented. Before the flight control design proceeds, a set of flight dynamics equations of motion must be formulated. Note that all aerodynamic forces and moments result from the relative motions between aircraft and the air. The aircraft is assumed to have a ground velocity:(26)
V
=
u
i
+
v
j
+
w
k
=
X
˙
E
I
+
Y
˙
E
J
-
h
˙
K
,where (
i
,
j
,
k
) are unit vectors in the aircraft body moving frame and (
I
,
J
,
K
) the unit vectors in a fixed flat earth frame. The air is assumed to have a velocity:(27)
V
w
=
u
g
i
+
v
g
j
+
w
g
k
=
V
w
x
I
+
V
w
y
J
-
V
h
Kwhich is also known as the wind velocity. Then, the velocity of the aircraft relative to the air can be represented by(28)
V
a
=
V
-
V
w
=
u
-
u
g
i
+
v
-
v
g
j
+
w
-
w
g
kfrom which, the aircraft total speed relative to the air, the angle of attack, and the sideslip angle can be determined, respectively, by the following equations:(29)
V
a
=
u
-
u
g
2
+
v
-
v
g
2
+
w
-
w
g
2
α
=
tan
-
1
⁡
w
-
w
g
u
-
u
g
β
=
tan
-
1
⁡
v
-
v
g
u
-
u
g
2
+
w
-
w
g
2
.With the assumptions of fixed flat earth and winds being present, the motions of aircraft with six degrees of freedom can be represented by a set of nonlinear first-order differential equations as follows:(30)
X
˙
E
≜
f
X
=
V
a
C
11
cos
⁡
α
cos
⁡
β
+
C
21
sin
⁡
β
+
C
31
sin
⁡
α
cos
⁡
β
+
V
w
x
(31)
Y
˙
E
≜
f
Y
=
V
a
C
12
cos
⁡
α
cos
⁡
β
+
C
22
sin
⁡
β
+
C
32
sin
⁡
α
cos
⁡
β
+
V
w
y
(32)
h
˙
≜
f
h
=
-
V
a
C
13
cos
⁡
α
cos
⁡
β
+
C
23
sin
⁡
β
+
C
33
sin
⁡
α
cos
⁡
β
+
V
h
(33)
ψ
˙
≜
f
ψ
=
q
sec
⁡
θ
sin
⁡
ϕ
+
r
sec
⁡
θ
cos
⁡
ϕ
(34)
θ
˙
≜
f
θ
=
q
cos
⁡
ϕ
-
r
sin
⁡
ϕ
(35)
ϕ
˙
≜
f
ϕ
=
p
+
q
tan
⁡
θ
sin
⁡
ϕ
+
r
tan
⁡
θ
cos
⁡
ϕ
(36)
V
˙
a
≜
f
V
=
F
x
m
cos
⁡
α
cos
⁡
β
+
F
y
m
sin
⁡
β
+
F
z
m
sin
⁡
α
cos
⁡
β
+
f
V
w
(37)
α
˙
≜
f
α
=
-
p
cos
⁡
α
tan
⁡
β
+
q
-
r
sin
⁡
α
tan
⁡
β
-
F
x
m
V
a
sin
⁡
α
sec
⁡
β
+
F
z
m
V
a
cos
⁡
α
sec
⁡
β
+
f
α
w
(38)
β
˙
≜
f
β
=
p
sin
⁡
α
-
r
cos
⁡
α
-
F
x
m
V
a
cos
⁡
α
sin
⁡
β
+
F
y
m
V
a
cos
⁡
β
-
F
z
m
V
a
sin
⁡
α
sin
⁡
β
+
f
β
w
(39)
p
˙
≜
f
p
=
I
x
z
-
I
x
x
+
I
y
y
-
I
z
z
I
x
x
I
z
z
-
I
x
z
I
x
z
p
q
+
I
y
y
I
z
z
-
I
z
z
I
z
z
-
I
x
z
I
x
z
I
x
x
I
z
z
-
I
x
z
I
x
z
q
r
+
I
z
z
I
x
x
I
z
z
-
I
x
z
I
x
z
M
x
-
I
x
z
I
x
x
I
z
z
-
I
x
z
I
x
z
M
z
(40)
q
˙
≜
f
q
=
I
x
z
I
y
y
p
2
-
r
2
+
I
z
z
-
I
x
x
I
y
y
p
r
+
1
I
y
y
M
y
(41)
r
˙
≜
f
r
=
I
x
x
I
x
x
-
I
x
x
I
y
y
+
I
x
z
I
x
z
I
x
x
I
z
z
-
I
x
z
I
x
z
p
q
+
I
x
z
I
x
x
-
I
y
y
+
I
z
z
I
x
x
I
z
z
-
I
x
z
I
x
z
q
r
-
I
x
z
I
x
x
I
z
z
-
I
x
z
I
x
z
M
x
+
I
x
x
I
x
x
I
z
z
-
I
x
z
I
x
z
M
z
,where (
X
E
,
Y
E
) represents the position in a fixed flat earth frame, h the altitude, C
i
j the elements of direction cosine matrix for transferring a fixed flat earth frame to the aircraft body moving frame, α the angle of attack, β the sideslip angle, ψ the heading angle, θ the pitch angle, ϕ the bank angle, p the roll rate, q the pitch rate, and r the yaw rate. Also, F
x, F
y, and F
z represent, respectively, three components of the total force in an aircraft body moving frame. The three force components are composed of the thrust T, the lift L, the drag D, the side force Y, and the weight m
g (m is the aircraft mass and g the gravity acceleration) by the following equations:(42)
F
x
=
T
cos
⁡
θ
T
+
L
sin
⁡
α
-
D
cos
⁡
α
-
m
g
sin
⁡
θ
F
y
=
Y
+
m
g
sin
⁡
ϕ
cos
⁡
θ
F
z
=
-
T
sin
⁡
θ
T
-
L
cos
⁡
α
-
D
sin
⁡
α
+
m
g
cos
⁡
ϕ
cos
⁡
θ
,where θ
T is the angle between the thrust and the longitudinal axis. Moreover, M
x, M
y, and M
z represent the roll moment, the pitch moment, and the yaw moment, respectively, about the center of gravity, and I
x
x, I
y
y, I
z
z, and I
x
z are the components of the moment-of-inertia tensor. Furthermore, the wind disturbances on V
˙
a, α
˙, and β
˙ are, respectively, represented by f
V
w, f
α
w, and f
β
w which are expressed as follows:(43)
f
V
w
=
F
x
w
m
cos
⁡
α
cos
⁡
β
+
F
y
w
m
sin
⁡
β
+
F
z
w
m
sin
⁡
α
cos
⁡
β
f
α
w
=
-
F
x
w
m
V
a
sin
⁡
α
sec
⁡
β
+
F
z
w
m
V
a
cos
⁡
α
sec
⁡
β
f
β
w
=
-
F
x
w
m
V
a
cos
⁡
α
sin
⁡
β
+
F
y
w
m
V
a
cos
⁡
β
-
F
z
w
m
V
a
sin
⁡
α
sin
⁡
β
,where(44)
F
x
w
=
m
r
v
g
-
q
w
g
-
u
˙
g
F
y
w
=
m
-
r
u
g
+
p
w
g
-
v
˙
g
F
z
w
=
m
q
u
g
-
p
v
g
-
w
˙
g
.In (44), F
x
w, F
y
w, and F
z
w are the three components of the force exerted by winds. At this point, it is worthy to mention that although there is no explicit term relating wind disturbances to p
˙, q
˙, and r
˙ in (39)–(41), winds do have effects on V
a, α, and β through which M
x, M
y, and M
z are affected. Also, winds not only have explicit effects on V
˙
a, α
˙, and β
˙ in (36)–(38), but also have implicit effects on them through L, Y, and D which obviously depend on V
a, α, and β.To validate the method developed in this paper, a UAV as shown in Figure1 is introduced. The parameters used for analysis are listed in Table 1.Table 1
The geometric data, weight, and moment of inertia of the UAV.
Parameter
Value
Units
Reference wing area
S
75
ft2
Wing span
b
15
ft
Mean aerodynamic chord
c
-
5.66
ft
Weight
W
2562.5
lb
x-axis inertia
I
x
x
296.75
slug/ft2
y-axis inertia
I
y
y
1744.1875
slug/ft2
z-axis inertia
I
z
z
1971.875
slug/ft2
x-z product of inertia
I
x
z
30.6875
slug/ft2Figure 1
The configuration of the UAV.The aerodynamic forces and moments of the UAV are computed by the following equations:(45)
L
=
q
-
S
C
L
M
,
α
,
δ
e
+
C
L
q
M
,
α
c
-
q
2
V
a
D
=
q
-
S
C
D
M
,
α
,
δ
e
Y
=
q
-
S
C
Y
β
M
,
α
β
+
C
Y
p
M
,
α
b
p
2
V
a
+
C
Y
r
M
,
α
b
r
2
V
a
+
Δ
C
Y
M
,
α
,
δ
a
+
Δ
C
Y
M
,
α
,
δ
r
M
x
=
q
-
S
b
C
l
β
M
,
α
β
+
C
l
p
M
,
α
b
p
2
V
a
+
C
l
r
M
,
α
b
r
2
V
a
+
Δ
C
l
M
,
α
,
δ
a
+
Δ
C
l
M
,
α
,
δ
r
M
y
=
q
-
S
c
-
C
m
M
,
α
,
δ
e
+
C
m
q
M
,
α
c
-
q
2
V
a
M
z
=
q
-
S
b
C
n
β
M
,
α
β
+
C
n
p
M
,
α
b
p
2
V
a
+
C
n
r
M
,
α
b
r
2
V
a
+
Δ
C
n
M
,
α
,
δ
a
+
Δ
C
n
M
,
α
,
δ
r
,where q
- is the dynamic pressure and M the Mach number. The aerodynamic coefficients and stability derivatives are determined by a computer code dubbed “VORSTAB” [17]. The computation results are in the form of discrete data which are then interpolated as continuous functions in the computer program. These functions are represented by C
L
M
,
α
,
δ
e, C
L
q
M
,
α, and so forth, in (45). The nonlinear database includes α ranging from -
20
° to 40
°, δ
e from -
24
° to 24
°, δ
a from -
25
° to 25
°, and δ
r from -
30
° to 30
°. Also, the flight conditions include altitude h ranging from sea level to 40,000ft and Mach number M from 0.1 to 0.9 and from 1.1 to 1.9. The thrust model of the UAV is assumed as(46)
T
˙
=
0.1
T
c
-
T
,where T
c is the thrust command. The thrust is assumed to have a limitation, say 0
≤
T
≤
1.1
W.


### 3.2. Flight Control Design

In flight control design, only nominal flight dynamics models are considered since modeling uncertainties or wind disturbances, and so forth, cannot be determined in advance. Hence, all terms related to the wind disturbances in the flight dynamics equations of motion are neglected in this stage. For this flight dynamics model, the parameters involved in the control analysis are further elaborated as follows.(1)
In the flight dynamics equations of motion, onlyh, θ, ϕ, V, α, β, p, q, and r are coupled. Therefore x
≜
h
,
θ
,
ϕ
,
V
,
α
,
β
,
p
,
q
,
r
T, f
=
f
h
,
f
θ
,
f
ϕ
,
f
V
,
f
α
,
f
β
,
f
p
,
f
q
,
f
r
T, and n
=
9.
(2)
The control input vectoru
≜
δ
e
,
δ
a
,
δ
r
,
T
c
T and m
=
4. According to the theory derived in the above, only 4 state variables can be controlled or designated. In this paper, the state vector to be controlled x
c
≜
h
c
,
V
c
,
ϕ
c
,
β
c
T is chosen and the residual vector is x
r
≜
θ
r
,
α
r
,
p
r
,
q
r
,
r
r
T.
(3)
Determine the state variables and control variables for some trim conditions, on which the control design is based.
(4)
The model matrix in (3) is constructed as(47)
Ω
=
σ
1
0
0
0
0
0
0
0
0
0
σ
2
0
σ
4
0
0
0
0
0
0
0
σ
3
0
0
0
0
0
0
0
σ
4
0
σ
2
0
0
0
0
0
0
0
0
0
σ
5
0
0
σ
8
0
0
0
0
0
0
σ
6
0
0
σ
9
0
0
0
0
0
0
σ
7
0
0
0
0
0
0
σ
8
0
0
σ
5
0
0
0
0
0
0
σ
9
0
0
σ
6
,where the arrangement for the pair (σ
2
,
σ
4) is enlightened from the phugoid mode in which θ and V are coupled. Similar arrangements are for the pairs (σ
5
,
σ
8) and (σ
6
,
σ
9). In order to make the desired dynamics stable, the elements on the diagonal, σ
1, σ
2, σ
3, σ
5, σ
6, and σ
7, must be positive. Also, those off-diagonal elements must satisfy the conditions, σ
4
2
<
σ
2
2, σ
8
2
<
σ
5
2, and σ
9
2
<
σ
6
2. For convenience of later usage, these necessary conditions are dubbed “the stable model requirements.”In the design, the initial cruise flight conditions areh
=
600ft and M
=
0.25. The trim conditions are α
t
r
i
m
=
θ
t
r
i
m
=
6.045
°, T
t
r
i
m
=
341.6lb, and δ
e
t
r
i
m
=
-
0.872
°. Based on these data, the open-loop eigenvalues of f
x
(
x
t
r
i
m
,
u
t
r
i
m
) are determined and listed in Table 2.Table 2
The open-loop eigenvalues off
x
x
trim
,
u
trim.
λ
1,2
−1.124786 ± 1.566875i
λ
6
-
3.129774
λ
3,4
−0.008857 ± 0.158058i
λ
7,8
−0.297956 ± 0.410071i
λ
5
-
0.000133
λ
9
-
0.023006Longitudinal dynamics analysis reveals thatλ
1,2 and λ
3,4 are closely associated with the short period mode and the phugoid mode, respectively, and lateral dynamics analysis shows that λ
6, λ
7,8, and λ
9 are closely associated with the pure roll mode, the Dutch roll mode, and the spiral mode, respectively. Finally, the rest one, λ
5, can be inferred to be associated with the altitude dynamics.Intuitively,σ
i
i
=
1,2
,
3
,
…
,
9
∈
σ may just be arbitrarily assigned as long as they satisfy the stable model requirements. However, the resulted error dynamics in (9) may not just be stable. As will be shown in simulations, the desired eigenvalues must not only be able to make the system follow the desired dynamics, but also be able to generate a good set of σ which makes e
-
Ω
t decrease so quickly that the system is also robust if modeling uncertainties exist or wind disturbances are encountered. The choice of these eigenvalues is not only just based on the control analysis alone but must also be simultaneously based on simulations.Note that, in this case, longitudinal dynamics analysis also reveals a right-half zeroz
=
10.85 in the h
-
δ
e transfer function. It can therefore be identified as a nonminimum phase problem. For this kind of problem, the desired eigenvalues must be very carefully assigned lest the resulting σ does not satisfy the stable model requirements. To determine σ, a two-way approach (σ
⇌
λ) is proposed as follows.(1)
If the desired eigenvalues are equal to the open-loop eigenvalues as listed in Table2, obviously, σ
=
0. Choose small σ
(
≠
0
) which satisfies the stable model requirements.
(2)
Determine the eigenvaluesλ
i
i
=
1,2
,
3
,
…
,
9
∈
λ of the error dynamics in (9).
(3)
If the error dynamics is not stable, modifyλ so that it makes the error dynamics stable. Use the modified values as the desired eigenvalues.
(4)
Determineσ by following the 7 steps of iteration procedure described in Section 2.2.
(5)
Ifσ does not satisfy the stable model requirements, modify σ to make it satisfy. Go to step 2.
(6)
In each step of simulations, determinex
r and u from (4) by using the iteration method described in (5). In this step, only the nominal dynamics model is used.
(7)
Useu in (1) for simulations. The equations of motion may include modeling uncertainties or wind disturbances.
(8)
If the simulation results are not satisfactory, modifyσ and then go to step 2, or modify λ and then go to step 4.By using these 8 steps of procedure, the designated eigenvalues and the resulting parameters for desired dynamics are obtained and listed in Tables3 and 4, respectively.Table 3
The designated eigenvalues.
λ
1,2
−4.22180 ± 1.57684i
λ
7
-
1.75509
λ
3,4
−0.112728 ± 0.0902199i
λ
8
-
0.920091
λ
5
-
0.0204641
λ
9
-
0.146833
λ
6
-
4.73168Table 4
The resulting parameters for desired dynamics.
σ
1
0.0201234
σ
6
1.04305
σ
2
0.101328
σ
7
1.60654
σ
3
0.112353
σ
8
0.00193758
σ
4
-
0.0000325585
σ
9
0.295443
σ
5
3.09966These parameters indeed satisfy the stable model requirements. Note that initiallyσ
5 is small and λ
1,2 and λ
3,4 are not so deviated from their counterparts of open-loop eigenvalues. However, the simulation results show that the performance in tracking h
c is not good enough. Enlarging σ
5 does improve the performance but makes λ
1,2 and λ
3,4 deviate their counterparts of open-loop eigenvalues a lot. Also note that the large deviations of λ
6, λ
7,8, and λ
9 from their counterparts of open-loop eigenvalues are due to the iteration procedure in determining a set of σ
6, σ
7, σ
8, and σ
9 for satisfying the stable model requirements.


### 3.3. Flight Simulations

In simulations, the initial conditions are(
X
E
,
Y
E
)
=
(
0
,
0
), h
=
600ft, ψ
=
0
°, θ
=
θ
t
r
i
m
=
6.04
°, ϕ
=
0
°, M
=
0.25, α
=
α
t
r
i
m
=
6.04
°, β
=
0
°, p
=
0, q
=
0, and r
=
0. The states to be tracked are h
c
=
2,000ft, ψ
c
=
90
°, M
c
=
0.40, and β
c
=
0
°. Although the heading angle is to be tracked, for pilots, it would make more sense to regulate the bank angle rather than the heading angle by assuming(48)
ϕ
c
=
k
ϕ
ψ
c
-
ψ
,where k
ϕ is a constant parameter. In this study, k
ϕ
=
0.12.At each instant,h
c, V
c, ϕ
c, and β
c are given, the controls δ
e, δ
a, δ
r, and T
c along with the residual state variables α
r, θ
r, p
r, q
r, and r
r can be determined from (4) by using the iteration method described in (5). Note that, in the computation of the control, only nominal flight dynamics equations of motion are used since modeling uncertainties or wind disturbances are not known. Also note that, in using the Newton-Raphson method, the convergence can be guaranteed if the initial guess is sufficiently close to the solution. Since the solution changes only very little in each step of integration, numerical practices show that it takes only 4 iterations to converge to within 0.0001
% of the correct value if the solution in the previous step of integration is taken as an initial guess. In the first step of integration, it may be necessary to take a few more iterations, say, 10 iterations, since the solution is different from that in the trim conditions, which is usually taken as an initial guess.For illustration, simulation block diagrams includingx, x
d, x
c, x
r, and u are shown in Figure 2.Figure 2
Simulation block diagrams for the NDI flight control system.After the controlsδ
e, δ
a, δ
r, and T
c are determined, (30)–(41) are used for simulations. In order to ascertain if the control system is robust or not, a fictitious wind model is put enroute as follows:(49)
V
w
x
=
-
30
×
6080
3600
sin
⁡
2
π
X
E
-
X
w
o
-
X
s
c
a
l
e
2
X
s
c
a
l
e
if
X
E
-
X
w
o
<
X
s
c
a
l
e
,
0
else
.
(50)
V
w
y
=
30
×
6080
3600
sin
⁡
2
π
Y
E
-
Y
w
o
-
Y
s
c
a
l
e
2
Y
s
c
a
l
e
if
Y
E
-
Y
w
o
<
Y
s
c
a
l
e
,
0
else
.
(51)
V
h
=
10
×
6080
3600
sin
⁡
2
π
h
-
H
w
o
-
H
s
c
a
l
e
2
H
s
c
a
l
e
if
h
-
H
w
o
<
H
s
c
a
l
e
,
0
else
.
(52)
V
˙
w
x
=
-
30
×
6080
3600
×
2
π
X
˙
E
2
X
s
c
a
l
e
cos
⁡
2
π
X
E
-
X
w
o
-
X
s
c
a
l
e
2
X
s
c
a
l
e
if
X
E
-
X
w
o
<
X
s
c
a
l
e
,
0
else
.
(53)
V
˙
w
y
=
30
×
6080
3600
×
2
π
Y
˙
E
2
Y
s
c
a
l
e
cos
⁡
2
π
Y
E
-
Y
w
o
-
Y
s
c
a
l
e
2
Y
s
c
a
l
e
if
Y
E
-
Y
w
o
<
Y
s
c
a
l
e
,
0
else
.
(54)
V
˙
h
=
10
×
6080
3600
×
2
π
h
˙
2
H
s
c
a
l
e
cos
⁡
2
π
h
-
H
w
o
-
H
s
c
a
l
e
2
H
s
c
a
l
e
if
h
-
H
w
o
<
H
s
c
a
l
e
,
0
else
,where X
w
o
=
29
,
727ft, Y
w
o
=
10
,
180ft, and H
w
o
=
1
,
000ft represent the center position of wind zone, and X
s
c
a
l
e
=
5
,
000ft, Y
s
c
a
l
e
=
5
,
000ft, and H
s
c
a
l
e
=
1
,
000ft represent the maximum range of wind zone from its center. From (49)–(51), it is also known that the three maximum wind velocity components in the earth fixed frame are 30 knots, 30 knots, and 10 knots, respectively.In this study, two cases of simulations are conducted; one is without wind disturbances and another with wind disturbances. The results are presented in Figures3–6.Figure 3
(a) Trajectory inX-Y plane, (b) wind speed, and (c) wind acceleration.
(a)
(b)
(c)Figure 4
The responses of (a) altitude, (b) Mach number, (c) heading angle, (d) bank angle, and (e) sideslip angle.
(a)
(b)
(c)
(d)
(e)Figure 5
The responses of (a) control surface deflections and (b) thrust.
(a)
(b)Figure 6
The responses of (a) angle of attack, (b) pitch angle, (c) pitch rate, (d) roll rate, and (e) yaw rate.
(a)
(b)
(c)
(d)
(e)As shown in Figure3(a), there is no much difference in ground trajectories whether wind disturbances are present or not. While flying, the UAV first suffers vertical wind disturbances in between t
=
15.10sec and 32.04sec and then horizontal wind disturbances in between t
=
50.60sec and 90.00sec as shown in Figures 3(b) and 3(c).The state variables to be tracked are all plotted in Figures4(a)–4(e). From Figure 4(a), it is found that h decreases initially when h
c
=
2
,
000ft is commanded. A careful study reveals that the elevator angle computed is δ
e
=
-
0.733
°, which is not enough to hold the UAV level as δ
e
t
r
i
m
=
-
0.872
° is required. The lowest altitude reached is 578ft with σ
5
=
3.09966 in Table 4 being used. Increasing σ
5 will increase the lowest altitude reached but at the expense of increasing overshoot. Fortunately, after descending to 578ft, the UAV begins to climb. During its climb, the UAV encounters an ascending wind and then a descending wind. But the influences on the altitude are almost negligible. In contrary, the horizontal wind seems to have more influences as it makes the altitude fluctuate. Whether wind disturbances are present or not, the UAV can always approach the commanded altitude asymptotically.As shown in Figure4(b), the Mach number does not seem to be so affected by the vertical wind as does by the horizontal wind. In this case, the horizontal wind causes the Mach number to fluctuate up to ±
10
%. The sharp angles shown in the figure are due to the discontinuity of wind acceleration at the edges of wind zone as shown in Figure 3(c).When the UAV is commanded to make a90
° turn, its heading angle gradually increases and asymptotically approaches the command as shown in Figure 4(c) whether winds are present or not. Figure 4(d) shows that the bank angle is closely associated with the heading angle as the bank angle is computed based on (48). Although the wind disturbances do not explicitly affect ϕ
˙ in (35), they do implicitly affect it through α and β which in turn affect the aerodynamic moments. Fluctuations in the bank angle and the heading angle are remarkable when the UAV passes through the wind zone but are still tolerable.Sideslip angle seems to be affected by winds more heavily as shown in Figure4(e). Fortunately, the control system is robust enough as being able to keep it small.As shown in Figure5(a), all control surface deflection angles remain small. The elevator needs only to deflect slightly to rotate the UAV to climb. As mentioned early, δ
e
=
-
0.733
° at the instant when h
c and M
c are given. When the UAV reaches the commanded state conditions, the elevator trim angle approaches -
0.510
°. Small deflection in aileron angle is enough to make the UAV bank turn and the rudder just keeps very small as does the sideslip angle. All control surface deflection angles are not remarkably affected as the UAV passes through the wind zone.As shown in Figure5(b), the thrust command increases very sharply as a demand to increase both the altitude and the Mach number simultaneously is given. Also, when the UAV passes the wind zone, the thrust command fluctuates very violently. Obviously, the thrust is very closely interacted with the Mach number. A low pass filter with time constant 10sec alleviates the sharp and violent responses a lot for the actual thrust at the expense of delaying its reaction time.In Figures6(a)–6(e), responses of all residual state variables along with the computed commands for them are presented. In fact, these commands are generated in the inner system, not input from outer designations. It is interesting to observe how the actual state variables track the commands.In Figure6(a), it is observed that, in comparison to α
t
r
i
m
=
6.04
°, the computed command α
r
=
4.85
° at the instant when h
c and M
c are given. Then the angle of attack α follows α
r closely without apparent short-period mode oscillations when the designated short-period mode eigenvalues are as high as λ
1,2
=
-
4.22180
±
1
.
57684
i. Winds do have remark effect on both α
r and α. Finally, both approach a new value, α
t
r
i
m
=
2.88
°, for flight conditions at higher altitude and faster speed.In Figures6(b) and 6(c), it is observed that although both θ
r and q
r increase in response to the requirement for increasing the cruise altitude, θ and q decrease remarkably. As revealed in explaining h response, the elevator angle computed is not enough to hold q
˙
=
0 and therefore q
˙
<
0 which in turn makes θ
˙
<
0. The minimum pitch rate is q
=
-
0.717deg./s at t
=
0.78sec. In this case, the vertical wind does have more remarkable effects on θ and q than the horizontal wind. Finally, both approach θ
t
r
i
m
=
α
t
r
i
m
=
2.88
° and q
t
r
i
m
=
0, respectively, for the new flight conditions.It is very interesting to find that, in Figures6(d) and 6(e), the shapes of p and r responses resemble, respectively, those of ϕ and β responses as shown in Figures 4(d) and 4(e). The effects of horizontal wind on both responses are more remarkable than those of vertical wind in this case. The decay of r seems to be very slow.


## 3.1. Flight Dynamics Equations of Motion

To illustrate the theory, a design of flight control system with the method developed is presented. Before the flight control design proceeds, a set of flight dynamics equations of motion must be formulated. Note that all aerodynamic forces and moments result from the relative motions between aircraft and the air. The aircraft is assumed to have a ground velocity:(26)
V
=
u
i
+
v
j
+
w
k
=
X
˙
E
I
+
Y
˙
E
J
-
h
˙
K
,where (
i
,
j
,
k
) are unit vectors in the aircraft body moving frame and (
I
,
J
,
K
) the unit vectors in a fixed flat earth frame. The air is assumed to have a velocity:(27)
V
w
=
u
g
i
+
v
g
j
+
w
g
k
=
V
w
x
I
+
V
w
y
J
-
V
h
Kwhich is also known as the wind velocity. Then, the velocity of the aircraft relative to the air can be represented by(28)
V
a
=
V
-
V
w
=
u
-
u
g
i
+
v
-
v
g
j
+
w
-
w
g
kfrom which, the aircraft total speed relative to the air, the angle of attack, and the sideslip angle can be determined, respectively, by the following equations:(29)
V
a
=
u
-
u
g
2
+
v
-
v
g
2
+
w
-
w
g
2
α
=
tan
-
1
⁡
w
-
w
g
u
-
u
g
β
=
tan
-
1
⁡
v
-
v
g
u
-
u
g
2
+
w
-
w
g
2
.With the assumptions of fixed flat earth and winds being present, the motions of aircraft with six degrees of freedom can be represented by a set of nonlinear first-order differential equations as follows:(30)
X
˙
E
≜
f
X
=
V
a
C
11
cos
⁡
α
cos
⁡
β
+
C
21
sin
⁡
β
+
C
31
sin
⁡
α
cos
⁡
β
+
V
w
x
(31)
Y
˙
E
≜
f
Y
=
V
a
C
12
cos
⁡
α
cos
⁡
β
+
C
22
sin
⁡
β
+
C
32
sin
⁡
α
cos
⁡
β
+
V
w
y
(32)
h
˙
≜
f
h
=
-
V
a
C
13
cos
⁡
α
cos
⁡
β
+
C
23
sin
⁡
β
+
C
33
sin
⁡
α
cos
⁡
β
+
V
h
(33)
ψ
˙
≜
f
ψ
=
q
sec
⁡
θ
sin
⁡
ϕ
+
r
sec
⁡
θ
cos
⁡
ϕ
(34)
θ
˙
≜
f
θ
=
q
cos
⁡
ϕ
-
r
sin
⁡
ϕ
(35)
ϕ
˙
≜
f
ϕ
=
p
+
q
tan
⁡
θ
sin
⁡
ϕ
+
r
tan
⁡
θ
cos
⁡
ϕ
(36)
V
˙
a
≜
f
V
=
F
x
m
cos
⁡
α
cos
⁡
β
+
F
y
m
sin
⁡
β
+
F
z
m
sin
⁡
α
cos
⁡
β
+
f
V
w
(37)
α
˙
≜
f
α
=
-
p
cos
⁡
α
tan
⁡
β
+
q
-
r
sin
⁡
α
tan
⁡
β
-
F
x
m
V
a
sin
⁡
α
sec
⁡
β
+
F
z
m
V
a
cos
⁡
α
sec
⁡
β
+
f
α
w
(38)
β
˙
≜
f
β
=
p
sin
⁡
α
-
r
cos
⁡
α
-
F
x
m
V
a
cos
⁡
α
sin
⁡
β
+
F
y
m
V
a
cos
⁡
β
-
F
z
m
V
a
sin
⁡
α
sin
⁡
β
+
f
β
w
(39)
p
˙
≜
f
p
=
I
x
z
-
I
x
x
+
I
y
y
-
I
z
z
I
x
x
I
z
z
-
I
x
z
I
x
z
p
q
+
I
y
y
I
z
z
-
I
z
z
I
z
z
-
I
x
z
I
x
z
I
x
x
I
z
z
-
I
x
z
I
x
z
q
r
+
I
z
z
I
x
x
I
z
z
-
I
x
z
I
x
z
M
x
-
I
x
z
I
x
x
I
z
z
-
I
x
z
I
x
z
M
z
(40)
q
˙
≜
f
q
=
I
x
z
I
y
y
p
2
-
r
2
+
I
z
z
-
I
x
x
I
y
y
p
r
+
1
I
y
y
M
y
(41)
r
˙
≜
f
r
=
I
x
x
I
x
x
-
I
x
x
I
y
y
+
I
x
z
I
x
z
I
x
x
I
z
z
-
I
x
z
I
x
z
p
q
+
I
x
z
I
x
x
-
I
y
y
+
I
z
z
I
x
x
I
z
z
-
I
x
z
I
x
z
q
r
-
I
x
z
I
x
x
I
z
z
-
I
x
z
I
x
z
M
x
+
I
x
x
I
x
x
I
z
z
-
I
x
z
I
x
z
M
z
,where (
X
E
,
Y
E
) represents the position in a fixed flat earth frame, h the altitude, C
i
j the elements of direction cosine matrix for transferring a fixed flat earth frame to the aircraft body moving frame, α the angle of attack, β the sideslip angle, ψ the heading angle, θ the pitch angle, ϕ the bank angle, p the roll rate, q the pitch rate, and r the yaw rate. Also, F
x, F
y, and F
z represent, respectively, three components of the total force in an aircraft body moving frame. The three force components are composed of the thrust T, the lift L, the drag D, the side force Y, and the weight m
g (m is the aircraft mass and g the gravity acceleration) by the following equations:(42)
F
x
=
T
cos
⁡
θ
T
+
L
sin
⁡
α
-
D
cos
⁡
α
-
m
g
sin
⁡
θ
F
y
=
Y
+
m
g
sin
⁡
ϕ
cos
⁡
θ
F
z
=
-
T
sin
⁡
θ
T
-
L
cos
⁡
α
-
D
sin
⁡
α
+
m
g
cos
⁡
ϕ
cos
⁡
θ
,where θ
T is the angle between the thrust and the longitudinal axis. Moreover, M
x, M
y, and M
z represent the roll moment, the pitch moment, and the yaw moment, respectively, about the center of gravity, and I
x
x, I
y
y, I
z
z, and I
x
z are the components of the moment-of-inertia tensor. Furthermore, the wind disturbances on V
˙
a, α
˙, and β
˙ are, respectively, represented by f
V
w, f
α
w, and f
β
w which are expressed as follows:(43)
f
V
w
=
F
x
w
m
cos
⁡
α
cos
⁡
β
+
F
y
w
m
sin
⁡
β
+
F
z
w
m
sin
⁡
α
cos
⁡
β
f
α
w
=
-
F
x
w
m
V
a
sin
⁡
α
sec
⁡
β
+
F
z
w
m
V
a
cos
⁡
α
sec
⁡
β
f
β
w
=
-
F
x
w
m
V
a
cos
⁡
α
sin
⁡
β
+
F
y
w
m
V
a
cos
⁡
β
-
F
z
w
m
V
a
sin
⁡
α
sin
⁡
β
,where(44)
F
x
w
=
m
r
v
g
-
q
w
g
-
u
˙
g
F
y
w
=
m
-
r
u
g
+
p
w
g
-
v
˙
g
F
z
w
=
m
q
u
g
-
p
v
g
-
w
˙
g
.In (44), F
x
w, F
y
w, and F
z
w are the three components of the force exerted by winds. At this point, it is worthy to mention that although there is no explicit term relating wind disturbances to p
˙, q
˙, and r
˙ in (39)–(41), winds do have effects on V
a, α, and β through which M
x, M
y, and M
z are affected. Also, winds not only have explicit effects on V
˙
a, α
˙, and β
˙ in (36)–(38), but also have implicit effects on them through L, Y, and D which obviously depend on V
a, α, and β.To validate the method developed in this paper, a UAV as shown in Figure1 is introduced. The parameters used for analysis are listed in Table 1.Table 1
The geometric data, weight, and moment of inertia of the UAV.
Parameter
Value
Units
Reference wing area
S
75
ft2
Wing span
b
15
ft
Mean aerodynamic chord
c
-
5.66
ft
Weight
W
2562.5
lb
x-axis inertia
I
x
x
296.75
slug/ft2
y-axis inertia
I
y
y
1744.1875
slug/ft2
z-axis inertia
I
z
z
1971.875
slug/ft2
x-z product of inertia
I
x
z
30.6875
slug/ft2Figure 1
The configuration of the UAV.The aerodynamic forces and moments of the UAV are computed by the following equations:(45)
L
=
q
-
S
C
L
M
,
α
,
δ
e
+
C
L
q
M
,
α
c
-
q
2
V
a
D
=
q
-
S
C
D
M
,
α
,
δ
e
Y
=
q
-
S
C
Y
β
M
,
α
β
+
C
Y
p
M
,
α
b
p
2
V
a
+
C
Y
r
M
,
α
b
r
2
V
a
+
Δ
C
Y
M
,
α
,
δ
a
+
Δ
C
Y
M
,
α
,
δ
r
M
x
=
q
-
S
b
C
l
β
M
,
α
β
+
C
l
p
M
,
α
b
p
2
V
a
+
C
l
r
M
,
α
b
r
2
V
a
+
Δ
C
l
M
,
α
,
δ
a
+
Δ
C
l
M
,
α
,
δ
r
M
y
=
q
-
S
c
-
C
m
M
,
α
,
δ
e
+
C
m
q
M
,
α
c
-
q
2
V
a
M
z
=
q
-
S
b
C
n
β
M
,
α
β
+
C
n
p
M
,
α
b
p
2
V
a
+
C
n
r
M
,
α
b
r
2
V
a
+
Δ
C
n
M
,
α
,
δ
a
+
Δ
C
n
M
,
α
,
δ
r
,where q
- is the dynamic pressure and M the Mach number. The aerodynamic coefficients and stability derivatives are determined by a computer code dubbed “VORSTAB” [17]. The computation results are in the form of discrete data which are then interpolated as continuous functions in the computer program. These functions are represented by C
L
M
,
α
,
δ
e, C
L
q
M
,
α, and so forth, in (45). The nonlinear database includes α ranging from -
20
° to 40
°, δ
e from -
24
° to 24
°, δ
a from -
25
° to 25
°, and δ
r from -
30
° to 30
°. Also, the flight conditions include altitude h ranging from sea level to 40,000ft and Mach number M from 0.1 to 0.9 and from 1.1 to 1.9. The thrust model of the UAV is assumed as(46)
T
˙
=
0.1
T
c
-
T
,where T
c is the thrust command. The thrust is assumed to have a limitation, say 0
≤
T
≤
1.1
W.


## 3.2. Flight Control Design

In flight control design, only nominal flight dynamics models are considered since modeling uncertainties or wind disturbances, and so forth, cannot be determined in advance. Hence, all terms related to the wind disturbances in the flight dynamics equations of motion are neglected in this stage. For this flight dynamics model, the parameters involved in the control analysis are further elaborated as follows.(1)
In the flight dynamics equations of motion, onlyh, θ, ϕ, V, α, β, p, q, and r are coupled. Therefore x
≜
h
,
θ
,
ϕ
,
V
,
α
,
β
,
p
,
q
,
r
T, f
=
f
h
,
f
θ
,
f
ϕ
,
f
V
,
f
α
,
f
β
,
f
p
,
f
q
,
f
r
T, and n
=
9.
(2)
The control input vectoru
≜
δ
e
,
δ
a
,
δ
r
,
T
c
T and m
=
4. According to the theory derived in the above, only 4 state variables can be controlled or designated. In this paper, the state vector to be controlled x
c
≜
h
c
,
V
c
,
ϕ
c
,
β
c
T is chosen and the residual vector is x
r
≜
θ
r
,
α
r
,
p
r
,
q
r
,
r
r
T.
(3)
Determine the state variables and control variables for some trim conditions, on which the control design is based.
(4)
The model matrix in (3) is constructed as(47)
Ω
=
σ
1
0
0
0
0
0
0
0
0
0
σ
2
0
σ
4
0
0
0
0
0
0
0
σ
3
0
0
0
0
0
0
0
σ
4
0
σ
2
0
0
0
0
0
0
0
0
0
σ
5
0
0
σ
8
0
0
0
0
0
0
σ
6
0
0
σ
9
0
0
0
0
0
0
σ
7
0
0
0
0
0
0
σ
8
0
0
σ
5
0
0
0
0
0
0
σ
9
0
0
σ
6
,where the arrangement for the pair (σ
2
,
σ
4) is enlightened from the phugoid mode in which θ and V are coupled. Similar arrangements are for the pairs (σ
5
,
σ
8) and (σ
6
,
σ
9). In order to make the desired dynamics stable, the elements on the diagonal, σ
1, σ
2, σ
3, σ
5, σ
6, and σ
7, must be positive. Also, those off-diagonal elements must satisfy the conditions, σ
4
2
<
σ
2
2, σ
8
2
<
σ
5
2, and σ
9
2
<
σ
6
2. For convenience of later usage, these necessary conditions are dubbed “the stable model requirements.”In the design, the initial cruise flight conditions areh
=
600ft and M
=
0.25. The trim conditions are α
t
r
i
m
=
θ
t
r
i
m
=
6.045
°, T
t
r
i
m
=
341.6lb, and δ
e
t
r
i
m
=
-
0.872
°. Based on these data, the open-loop eigenvalues of f
x
(
x
t
r
i
m
,
u
t
r
i
m
) are determined and listed in Table 2.Table 2
The open-loop eigenvalues off
x
x
trim
,
u
trim.
λ
1,2
−1.124786 ± 1.566875i
λ
6
-
3.129774
λ
3,4
−0.008857 ± 0.158058i
λ
7,8
−0.297956 ± 0.410071i
λ
5
-
0.000133
λ
9
-
0.023006Longitudinal dynamics analysis reveals thatλ
1,2 and λ
3,4 are closely associated with the short period mode and the phugoid mode, respectively, and lateral dynamics analysis shows that λ
6, λ
7,8, and λ
9 are closely associated with the pure roll mode, the Dutch roll mode, and the spiral mode, respectively. Finally, the rest one, λ
5, can be inferred to be associated with the altitude dynamics.Intuitively,σ
i
i
=
1,2
,
3
,
…
,
9
∈
σ may just be arbitrarily assigned as long as they satisfy the stable model requirements. However, the resulted error dynamics in (9) may not just be stable. As will be shown in simulations, the desired eigenvalues must not only be able to make the system follow the desired dynamics, but also be able to generate a good set of σ which makes e
-
Ω
t decrease so quickly that the system is also robust if modeling uncertainties exist or wind disturbances are encountered. The choice of these eigenvalues is not only just based on the control analysis alone but must also be simultaneously based on simulations.Note that, in this case, longitudinal dynamics analysis also reveals a right-half zeroz
=
10.85 in the h
-
δ
e transfer function. It can therefore be identified as a nonminimum phase problem. For this kind of problem, the desired eigenvalues must be very carefully assigned lest the resulting σ does not satisfy the stable model requirements. To determine σ, a two-way approach (σ
⇌
λ) is proposed as follows.(1)
If the desired eigenvalues are equal to the open-loop eigenvalues as listed in Table2, obviously, σ
=
0. Choose small σ
(
≠
0
) which satisfies the stable model requirements.
(2)
Determine the eigenvaluesλ
i
i
=
1,2
,
3
,
…
,
9
∈
λ of the error dynamics in (9).
(3)
If the error dynamics is not stable, modifyλ so that it makes the error dynamics stable. Use the modified values as the desired eigenvalues.
(4)
Determineσ by following the 7 steps of iteration procedure described in Section 2.2.
(5)
Ifσ does not satisfy the stable model requirements, modify σ to make it satisfy. Go to step 2.
(6)
In each step of simulations, determinex
r and u from (4) by using the iteration method described in (5). In this step, only the nominal dynamics model is used.
(7)
Useu in (1) for simulations. The equations of motion may include modeling uncertainties or wind disturbances.
(8)
If the simulation results are not satisfactory, modifyσ and then go to step 2, or modify λ and then go to step 4.By using these 8 steps of procedure, the designated eigenvalues and the resulting parameters for desired dynamics are obtained and listed in Tables3 and 4, respectively.Table 3
The designated eigenvalues.
λ
1,2
−4.22180 ± 1.57684i
λ
7
-
1.75509
λ
3,4
−0.112728 ± 0.0902199i
λ
8
-
0.920091
λ
5
-
0.0204641
λ
9
-
0.146833
λ
6
-
4.73168Table 4
The resulting parameters for desired dynamics.
σ
1
0.0201234
σ
6
1.04305
σ
2
0.101328
σ
7
1.60654
σ
3
0.112353
σ
8
0.00193758
σ
4
-
0.0000325585
σ
9
0.295443
σ
5
3.09966These parameters indeed satisfy the stable model requirements. Note that initiallyσ
5 is small and λ
1,2 and λ
3,4 are not so deviated from their counterparts of open-loop eigenvalues. However, the simulation results show that the performance in tracking h
c is not good enough. Enlarging σ
5 does improve the performance but makes λ
1,2 and λ
3,4 deviate their counterparts of open-loop eigenvalues a lot. Also note that the large deviations of λ
6, λ
7,8, and λ
9 from their counterparts of open-loop eigenvalues are due to the iteration procedure in determining a set of σ
6, σ
7, σ
8, and σ
9 for satisfying the stable model requirements.


## 3.3. Flight Simulations

In simulations, the initial conditions are(
X
E
,
Y
E
)
=
(
0
,
0
), h
=
600ft, ψ
=
0
°, θ
=
θ
t
r
i
m
=
6.04
°, ϕ
=
0
°, M
=
0.25, α
=
α
t
r
i
m
=
6.04
°, β
=
0
°, p
=
0, q
=
0, and r
=
0. The states to be tracked are h
c
=
2,000ft, ψ
c
=
90
°, M
c
=
0.40, and β
c
=
0
°. Although the heading angle is to be tracked, for pilots, it would make more sense to regulate the bank angle rather than the heading angle by assuming(48)
ϕ
c
=
k
ϕ
ψ
c
-
ψ
,where k
ϕ is a constant parameter. In this study, k
ϕ
=
0.12.At each instant,h
c, V
c, ϕ
c, and β
c are given, the controls δ
e, δ
a, δ
r, and T
c along with the residual state variables α
r, θ
r, p
r, q
r, and r
r can be determined from (4) by using the iteration method described in (5). Note that, in the computation of the control, only nominal flight dynamics equations of motion are used since modeling uncertainties or wind disturbances are not known. Also note that, in using the Newton-Raphson method, the convergence can be guaranteed if the initial guess is sufficiently close to the solution. Since the solution changes only very little in each step of integration, numerical practices show that it takes only 4 iterations to converge to within 0.0001
% of the correct value if the solution in the previous step of integration is taken as an initial guess. In the first step of integration, it may be necessary to take a few more iterations, say, 10 iterations, since the solution is different from that in the trim conditions, which is usually taken as an initial guess.For illustration, simulation block diagrams includingx, x
d, x
c, x
r, and u are shown in Figure 2.Figure 2
Simulation block diagrams for the NDI flight control system.After the controlsδ
e, δ
a, δ
r, and T
c are determined, (30)–(41) are used for simulations. In order to ascertain if the control system is robust or not, a fictitious wind model is put enroute as follows:(49)
V
w
x
=
-
30
×
6080
3600
sin
⁡
2
π
X
E
-
X
w
o
-
X
s
c
a
l
e
2
X
s
c
a
l
e
if
X
E
-
X
w
o
<
X
s
c
a
l
e
,
0
else
.
(50)
V
w
y
=
30
×
6080
3600
sin
⁡
2
π
Y
E
-
Y
w
o
-
Y
s
c
a
l
e
2
Y
s
c
a
l
e
if
Y
E
-
Y
w
o
<
Y
s
c
a
l
e
,
0
else
.
(51)
V
h
=
10
×
6080
3600
sin
⁡
2
π
h
-
H
w
o
-
H
s
c
a
l
e
2
H
s
c
a
l
e
if
h
-
H
w
o
<
H
s
c
a
l
e
,
0
else
.
(52)
V
˙
w
x
=
-
30
×
6080
3600
×
2
π
X
˙
E
2
X
s
c
a
l
e
cos
⁡
2
π
X
E
-
X
w
o
-
X
s
c
a
l
e
2
X
s
c
a
l
e
if
X
E
-
X
w
o
<
X
s
c
a
l
e
,
0
else
.
(53)
V
˙
w
y
=
30
×
6080
3600
×
2
π
Y
˙
E
2
Y
s
c
a
l
e
cos
⁡
2
π
Y
E
-
Y
w
o
-
Y
s
c
a
l
e
2
Y
s
c
a
l
e
if
Y
E
-
Y
w
o
<
Y
s
c
a
l
e
,
0
else
.
(54)
V
˙
h
=
10
×
6080
3600
×
2
π
h
˙
2
H
s
c
a
l
e
cos
⁡
2
π
h
-
H
w
o
-
H
s
c
a
l
e
2
H
s
c
a
l
e
if
h
-
H
w
o
<
H
s
c
a
l
e
,
0
else
,where X
w
o
=
29
,
727ft, Y
w
o
=
10
,
180ft, and H
w
o
=
1
,
000ft represent the center position of wind zone, and X
s
c
a
l
e
=
5
,
000ft, Y
s
c
a
l
e
=
5
,
000ft, and H
s
c
a
l
e
=
1
,
000ft represent the maximum range of wind zone from its center. From (49)–(51), it is also known that the three maximum wind velocity components in the earth fixed frame are 30 knots, 30 knots, and 10 knots, respectively.In this study, two cases of simulations are conducted; one is without wind disturbances and another with wind disturbances. The results are presented in Figures3–6.Figure 3
(a) Trajectory inX-Y plane, (b) wind speed, and (c) wind acceleration.
(a)
(b)
(c)Figure 4
The responses of (a) altitude, (b) Mach number, (c) heading angle, (d) bank angle, and (e) sideslip angle.
(a)
(b)
(c)
(d)
(e)Figure 5
The responses of (a) control surface deflections and (b) thrust.
(a)
(b)Figure 6
The responses of (a) angle of attack, (b) pitch angle, (c) pitch rate, (d) roll rate, and (e) yaw rate.
(a)
(b)
(c)
(d)
(e)As shown in Figure3(a), there is no much difference in ground trajectories whether wind disturbances are present or not. While flying, the UAV first suffers vertical wind disturbances in between t
=
15.10sec and 32.04sec and then horizontal wind disturbances in between t
=
50.60sec and 90.00sec as shown in Figures 3(b) and 3(c).The state variables to be tracked are all plotted in Figures4(a)–4(e). From Figure 4(a), it is found that h decreases initially when h
c
=
2
,
000ft is commanded. A careful study reveals that the elevator angle computed is δ
e
=
-
0.733
°, which is not enough to hold the UAV level as δ
e
t
r
i
m
=
-
0.872
° is required. The lowest altitude reached is 578ft with σ
5
=
3.09966 in Table 4 being used. Increasing σ
5 will increase the lowest altitude reached but at the expense of increasing overshoot. Fortunately, after descending to 578ft, the UAV begins to climb. During its climb, the UAV encounters an ascending wind and then a descending wind. But the influences on the altitude are almost negligible. In contrary, the horizontal wind seems to have more influences as it makes the altitude fluctuate. Whether wind disturbances are present or not, the UAV can always approach the commanded altitude asymptotically.As shown in Figure4(b), the Mach number does not seem to be so affected by the vertical wind as does by the horizontal wind. In this case, the horizontal wind causes the Mach number to fluctuate up to ±
10
%. The sharp angles shown in the figure are due to the discontinuity of wind acceleration at the edges of wind zone as shown in Figure 3(c).When the UAV is commanded to make a90
° turn, its heading angle gradually increases and asymptotically approaches the command as shown in Figure 4(c) whether winds are present or not. Figure 4(d) shows that the bank angle is closely associated with the heading angle as the bank angle is computed based on (48). Although the wind disturbances do not explicitly affect ϕ
˙ in (35), they do implicitly affect it through α and β which in turn affect the aerodynamic moments. Fluctuations in the bank angle and the heading angle are remarkable when the UAV passes through the wind zone but are still tolerable.Sideslip angle seems to be affected by winds more heavily as shown in Figure4(e). Fortunately, the control system is robust enough as being able to keep it small.As shown in Figure5(a), all control surface deflection angles remain small. The elevator needs only to deflect slightly to rotate the UAV to climb. As mentioned early, δ
e
=
-
0.733
° at the instant when h
c and M
c are given. When the UAV reaches the commanded state conditions, the elevator trim angle approaches -
0.510
°. Small deflection in aileron angle is enough to make the UAV bank turn and the rudder just keeps very small as does the sideslip angle. All control surface deflection angles are not remarkably affected as the UAV passes through the wind zone.As shown in Figure5(b), the thrust command increases very sharply as a demand to increase both the altitude and the Mach number simultaneously is given. Also, when the UAV passes the wind zone, the thrust command fluctuates very violently. Obviously, the thrust is very closely interacted with the Mach number. A low pass filter with time constant 10sec alleviates the sharp and violent responses a lot for the actual thrust at the expense of delaying its reaction time.In Figures6(a)–6(e), responses of all residual state variables along with the computed commands for them are presented. In fact, these commands are generated in the inner system, not input from outer designations. It is interesting to observe how the actual state variables track the commands.In Figure6(a), it is observed that, in comparison to α
t
r
i
m
=
6.04
°, the computed command α
r
=
4.85
° at the instant when h
c and M
c are given. Then the angle of attack α follows α
r closely without apparent short-period mode oscillations when the designated short-period mode eigenvalues are as high as λ
1,2
=
-
4.22180
±
1
.
57684
i. Winds do have remark effect on both α
r and α. Finally, both approach a new value, α
t
r
i
m
=
2.88
°, for flight conditions at higher altitude and faster speed.In Figures6(b) and 6(c), it is observed that although both θ
r and q
r increase in response to the requirement for increasing the cruise altitude, θ and q decrease remarkably. As revealed in explaining h response, the elevator angle computed is not enough to hold q
˙
=
0 and therefore q
˙
<
0 which in turn makes θ
˙
<
0. The minimum pitch rate is q
=
-
0.717deg./s at t
=
0.78sec. In this case, the vertical wind does have more remarkable effects on θ and q than the horizontal wind. Finally, both approach θ
t
r
i
m
=
α
t
r
i
m
=
2.88
° and q
t
r
i
m
=
0, respectively, for the new flight conditions.It is very interesting to find that, in Figures6(d) and 6(e), the shapes of p and r responses resemble, respectively, those of ϕ and β responses as shown in Figures 4(d) and 4(e). The effects of horizontal wind on both responses are more remarkable than those of vertical wind in this case. The decay of r seems to be very slow.


## 4. Conclusions

In this paper, a unified nonlinear dynamic inversion control system is successfully developed. With this method, the parameters for desired dynamics can be determined with a set of assigned eigenvalues for error dynamics. The control system has been proved to be robust to modeling uncertainties or wind disturbances. A nonaffine nonlinear flight dynamics system with right-half zero has been used as an example for the control design. In the design process, it is not necessary to use the feedback linearization to transform the nonlinear equations into linear equations. Numerical simulations of the control system show that the desired state variables can be successfully tracked whether winds are present or not.



---

*Source: 548050-2015-12-07.xml*