System Identification of General Aviation Aircraft
Using The Filter Error Technique
Except where reference is made to the work of others, the work described in this thesis is
my own or was done in collaboration with my advisory committee. This thesis does not
include proprietary or classi ed information.
Dakshesh Patel
Certi cate of Approval:
John E. Cochran, Jr.
Professor
Aerospace Engineering
Gilbert L. Crouse, Chair
Associate Professor
Aerospace Engineering
David A. Cicci
Professor
Aerospace Engineering
Joe F. Pittman
Interim Dean
Graduate School
System Identification of General Aviation Aircraft
Using The Filter Error Technique
Dakshesh Patel
A Thesis
Submitted to
the Graduate Faculty of
Auburn University
in Partial Ful llment of the
Requirements for the
Degree of
Master of Science
Auburn, Alabama
December 17, 2007
System Identification of General Aviation Aircraft
Using The Filter Error Technique
Dakshesh Patel
Permission is granted to Auburn University to make copies of this thesis at its
discretion, upon the request of individuals or institutions and at
their expense. The author reserves all publication rights.
Signature of Author
Date of Graduation
iii
Vita
Dakshesh Patel was born to Mahesh and Geeta, on March 15, 1983 in Ahmedabad.
He attended Nirma Institute of Technology, Ahmedabad and graduated in July 2004 with
Bachelor of Engineering in Instrumentation and Control. He joined Aerospace Engineering
Program at Auburn University as a Masters student in August 2004.
iv
Thesis Abstract
System Identification of General Aviation Aircraft
Using The Filter Error Technique
Dakshesh Patel
Master of Science, December 17, 2007
(B.E., Nirma Institute of Technology, Gujarat University, 2004)
123 Typed Pages
Directed by Gilbert L. Crouse
An overview of some past and modern applications of system identi cation techniques
to aircraft is presented in this thesis, which traces progress and accomplishments in the
eld of system identi cation - especially application to aircraft. The thesis also includes
the de nition, meaning, importance, applications and necessity of the same. The funda-
mental parts of system identi cation, including model postulation, experimental design,
data analysis, parameter estimation, and model validation are explained. The lter error
technique of data analysis is detailed, for both linear and nonlinear systems. Techniques to
predict aircraft stability and control derivatives theoretically and analytically are explained
in some details. Two successfully investigated experiments of system identi cation using
ight data with results and comparison of the estimated parameters with published data
are demonstrated with details. The thesis is completed by some remarks made on possible
future advancement and concluding notes.
v
Acknowledgments
The author would like to thank his research advisor, Dr. Gilbert L. Crouse, for render-
ing constant guidance and motivation throughout his research. He also thanks Dr. Crouse
for showing complete trust in his abilities and providing him with an opportunity to con-
duct research in the exciting eld of system identi cation. Dr. Crouse is a great mentor
and working with him was a wonderful and memorable experience for Dakshesh. The au-
thor would like to express his gratitude to Dr. John E. Cochran for teaching the courses
which greatly helped him in understanding the research problem and also for being part
of his graduate committee. He would like to thank Dr. David A. Cicci, for teaching the
fundamentals of estimation theory through orbit determination class and also for serving
as a graduate committee member.
The author is thankful to Ran Dai and Daroe Lee, for their friendship and help in
understanding the dynamics and control subjects. He would like to acknowledge his friends
at Auburn, with whom he spent memorable years. Particularly Palak, Vinita and Dhrumil?s
support and encouragement have always made his stay at Auburn, pleasant and enthusiastic.
The author wishes to dedicate this work to his parents and teachers for their enduring
love, immense moral support and encouragement in the journey of life.
vi
Style manual or journal used AIAA Journal of Aircraft
Computer software used LATEX and the departmental style- le \aums.sty"
vii
Table of Contents
List of Figures x
List of Tables xii
Nomenclature xiii
1 Introduction 1
1.1 System Identi cation - A De nition . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Signi cance and Necessity of System Identi cation . . . . . . . . . . 2
1.1.2 Objectives of Research . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Problem Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Progress of Aircraft System identi cation . . . . . . . . . . . . . . . . . . . 7
1.3.1 System Identi cation in General . . . . . . . . . . . . . . . . . . . . 7
1.3.2 System Identi cation of Aircraft . . . . . . . . . . . . . . . . . . . . 8
1.3.3 Deterministic and Non-Deterministic Analysis . . . . . . . . . . . . . 10
1.3.4 Latest Techniques of System Identi cation . . . . . . . . . . . . . . . 13
2 Background and Theoretical Development 16
2.1 Dynamical Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2 Cost Function (Performance Index) J . . . . . . . . . . . . . . . . . . . . . . 17
2.3 Essentials of Cost Function Minimization . . . . . . . . . . . . . . . . . . . 17
2.4 Modi ed Newton-Raphson Algorithm . . . . . . . . . . . . . . . . . . . . . 20
3 System Identification of General Aviation Aircraft 23
3.1 Filter Error Technique for Linear Systems . . . . . . . . . . . . . . . . . . . 23
3.1.1 Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1.2 Solution of Riccati Equation . . . . . . . . . . . . . . . . . . . . . . 26
3.1.3 Formulations for Process Noise . . . . . . . . . . . . . . . . . . . . . 28
3.1.4 The Filter Error Algorithm for Linear Systems . . . . . . . . . . . . 30
3.1.5 Parameter Update . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2 Filter Error Technique for Non-linear Systems . . . . . . . . . . . . . . . . . 36
3.2.1 Steady State Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.2.2 Parameter Update . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
viii
4 Prediction of Aircraft Stability and Control Derivatives 42
4.1 Steady State Coe cients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.2 Stability Derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2.1 Aircraft Speed Derivatives . . . . . . . . . . . . . . . . . . . . . . . . 46
4.2.2 Angle of Attack Derivatives . . . . . . . . . . . . . . . . . . . . . . . 47
4.2.3 Angle of Sideslip Derivatives . . . . . . . . . . . . . . . . . . . . . . 48
4.2.4 Roll Rate Derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.2.5 Pitch Rate Derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.2.6 Yaw Rate Derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.3 Control Derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.3.1 Aileron Control Derivatives . . . . . . . . . . . . . . . . . . . . . . . 56
4.3.2 Elevator Control Derivatives . . . . . . . . . . . . . . . . . . . . . . 57
4.3.3 Rudder Control Derivatives . . . . . . . . . . . . . . . . . . . . . . . 57
5 Investigated Examples of System Identification Using Flight Test Data 59
5.1 Estimation of Lateral Motion Derivatives using the Simulated Flight Data . 59
5.1.1 Comparison of the Results with the Original Experiment . . . . . . 61
5.1.2 Responses and Results . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.2 Estimation of Longitudinal Motion Derivatives of HFB-320 Aircraft . . . . 67
5.2.1 Comparison of the Results with the Original Experiment . . . . . . 70
5.2.2 Responses and Results . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.3 System Identi cation Applied to Model Aircraft Trainer . . . . . . . . . . . 73
5.3.1 The Investigated Aircraft Trainer - Hanger 9 Extra Easy . . . . . . 73
5.3.2 Responses and Results . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.4 Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6 Summary 88
6.1 Concluding Notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.2 Expected Future Applications of System Identi cation to Aircraft . . . . . . 89
Bibliography 91
A Discretization of continuous-time state equation 99
A.1 Need of Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
A.2 Determine and . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
B Numerical Integration using Runge-Kutta method, Order 4 102
ix
List of Figures
1.1 Block Diagram for System Identi cation to aircraft [1] . . . . . . . . . . . . 14
3.1 Block Diagram for Filter Error Technique [68] . . . . . . . . . . . . . . . . . 24
3.2 Algorithm of Filter Error Technique for Linear Systems [1] . . . . . . . . . . 31
3.3 Algorithm of Filter Error Technique for Nonlinear Systems [1] . . . . . . . . 37
4.1 Determination of Drag Coe cient from a known Lift Coe cient . . . . . . 44
5.1 The DLR de Havilland DHC II Research Aircraft [84] . . . . . . . . . . 60
5.2 Comparison of Estimated(green) and Measured(blue) Responses for the Lat-
eral Directional Motion of Aircraft . . . . . . . . . . . . . . . . . . . . . . . 64
5.3 Convergence Plots of Dimensional Derivatives for Lateral Directional Motion
of Aircraft . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.4 Cost Function Reduction for Lateral Directional Motion of Aircraft . . . . . 66
5.5 The DLR HFB 320 Research Aircraft [87] . . . . . . . . . . . . . . . . . 67
5.6 Comparison of Estimated(green) and Measured(blue) Responses for the Lon-
gitudinal Directional Motion of Aircraft . . . . . . . . . . . . . . . . . . . . 71
5.7 Convergence Plots of Non-Dimensional Derivatives for Longitudinal Direc-
tional Motion of Aircraft . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.8 Inputs for Longitudinal Motion . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.9 Cost Function Reduction for Longitudinal Motion . . . . . . . . . . . . . . 73
5.10 Aircraft Trainer - Hanger 9 Extra Easy [91] . . . . . . . . . . . . . . . . . . 74
5.11 Comparison of Estimated(green) and Measured(blue) Responses for the Lat-
eral Directional Motion of Aircraft . . . . . . . . . . . . . . . . . . . . . . . 77
x
5.12 Convergence Plots of Dimensional Derivatives for Lateral Directional Motion
of Aircraft . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.13 Cost Function Reduction for Lateral Directional Motion of Aircraft . . . . . 79
5.14 Comparison of Estimated(green) and Measured(blue) Responses for the Lat-
eral Directional Motion of Aircraft . . . . . . . . . . . . . . . . . . . . . . . 80
5.15 Convergence Plots of Dimensional Derivatives for Lateral Directional Motion
of Aircraft . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.16 Cost Function Reduction for Lateral Directional Motion of Aircraft . . . . . 82
5.17 Comparison of Estimated(green) and Measured(blue) Responses for the Lon-
gitudinal Directional Motion of Aircraft . . . . . . . . . . . . . . . . . . . . 84
5.18 Convergence Plots of Non-Dimensional Derivatives for Longitudinal Direc-
tional Motion of Aircraft . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.19 Input for Longitudinal Motion . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.20 Cost Function Reduction for Longitudinal Motion . . . . . . . . . . . . . . 86
xi
List of Tables
5.1 Estimated dimensional derivatives of lateral motion of aircraft . . . . . . . . 62
5.2 Estimated non-dimensional derivatives of longitudinal motion of aircraft . . 69
5.3 Aerodynamic, geometric and mass data for the aircraft trainer . . . . . . . 75
5.4 Comparison between the predicted values and the estimated values of the
dimensional derivatives for the lateral motion of the aircraft trainer . . . . . 76
5.5 Comparison between the predicted values and the estimated values of the
dimensional derivatives for the lateral motion of the aircraft trainer . . . . . 82
5.6 Comparison between the predicted values and the estimated values of the
non-dimensional derivatives for the longitudinal motion of aircraft trainer . 83
xii
Nomenclature
1. English Symbols
a:c: aerodynamic center
A state matrix of a linear system, size n n
Ah = bh2 = Sh horizontal tail aspect ratio
Av = bv2 = Sv vertical tail aspect ratio
Aw = b2 = S wing aspect ratio
ax, ay, az accelerations along x, y, z body axes, ft=s2
b wing span
bh horizontal tail span
bv vertical tail span
bx, by bias parameters of state and observation equations
c chord, ft
c mean geometric chord, ft
cr root chord, ft
ct tip chord, ft
c:g: aircraft center of gravity
C observation matrix, size m n
CD aircraft drag coe cient
CD0 zero-lift drag coe cient
xiii
Cl aerodynamic rolling moment coe cient
CL aircraft lift coe cient
CL0 lift coe cient at zero angle of attack
Cm aerodynamic pitching moment coe cient
Cn aerodynamic yawing moment coe cient
CX coe cient of longitudinal force
CY coe cient of lateral(side) force
CZ coe cient of vertical force
D Drag
e span e ciency factor
ej jth unit vector
F state noise matrix, size n n
also: aerodynamic force, N
FG magnitude of force due to gravity, N
FT magnitude of force due to thrust, N
Fe net thrust magnitude, N
f[ ] vector of system state functions
G measurement noise matrix, size m m
also: moment, Nm
g acceleration due to gravity, ft=s2
g[ ] system observation function
xiv
I inertia matrix
Ixx rolling moment of inertia in body axes, slug-ft2
Ixz XZ product of inertia in body axes, slug-ft2
Iyy pitching moment of inertia in body axes, slug-ft2
Izz yawing moment of inertia in body axes, slug-ft2
J cost function
K lter-gain matrix, size n m
k discrete-time index
lf fuselage length
lp horizontal distance between vertical tail a.c. and aircraft wing a.c.
lv horizontal distance between vertical tail a.c. and aircraft c.g.
L lift
m aircraft mass, lb
also: number of observation variables
m:g:c: mean geometric chord, ft
M free stream mach number
N number of data points
n number of state variables
P covariance matrix of the state error, size n n
p, q, r roll, pitch, and yaw rates, rad=s
_p, _q, _r roll, pitch, and yaw accelerations, rad=s2
xv
q dynamic pressure, N=ft2
rk kth diagonal element of R
R covariance matrix of residuals (innovations), size m m
RN reynold?s number
S wing area, ft2
Sbfuse fuselage base area, ft2
Sfuse maximum fuselage cross section area, ft2
Sh horizontal tail area, ft2
Splf aircraft planform area, ft2
Sv vertical tail area, ft2
Swet total wetted area of aircraft, ft2
s number of input variables
t time, s
t=c thickness ratio at c
u control input vector, size s 1
u, v, w velocity components along x, y, and z body axes, ft=s
v measurement noise vector, size m 1
V airspeed, ft=s
Vh horizontal tail volume coe cient
w state noise vector, size n 1
W aircraft weight, lb
xvi
x state vector, size n 1
Xac longitudinal position of a.c. on wing m.g.c., ft
Xach longitudinal position of horizontal tail a.c. on wing m.g.c., ft
Xw longitudinal distance from wing quarter chord m.g.c. to c.g., ft
y observation vector, size m 1
z measurement vector, size m 1
zf vertical height of fuselage at wing root chord, ft
zp vertical distance between vertical tail a.c. and aircraft wing a.c., ft
zv vertical distance between vertical tail a.c. and aircraft c.g., ft
zw wing distance to fuselage centerline, ft
2. Stability and Control Derivatives
CDV non-dimensional derivative of drag with respect to speed
CLV non-dimensional derivative of lift with respect to speed
CmV non-dimensional derivative of pitching moment with respect to speed
CD non-dimensional derivative of drag with respect to angle of attack
CL non-dimensional derivative of lift with respect to angle of attack
Cm non-dimensional derivative of pitching moment with respect to angle of attack
Cmq derivative of pitching moment with respect to pitch rate
Cm e derivative of pitching moment with respect to elevator de ection
Lp derivative of rolling moment with respect to roll rate
xvii
Lr derivative of rolling moment with respect to yaw rate
L a derivative of rolling moment with respect to aileron de ection
L r derivative of rolling moment with respect to rudder de ection
Lv derivative of rolling moment with respect to lateral(side) velocity
L derivative of rolling moment with respect to sideslip angle
Np derivative of yawing moment with respect to roll rate
Nr derivative of yawing moment with respect to yaw rate
N a derivative of yawing moment with respect to aileron de ection
N r derivative of yawing moment with respect to rudder de ection
Nv derivative of yawing moment with respect to lateral(side) velocity
N derivative of yawing moment with respect to sideslip angle
Yp derivative of side-force with respect to roll rate
Yr derivative of side-force with respect to yaw rate
Y a derivative of side-force with respect to aileron de ection
Y r derivative of side-force with respect to rudder de ection
Yv derivative of side-force with respect to lateral(side) velocity
Y derivative of side-force with respect to sideslip angle
3. Greek Symbols
angle of attack, deg
sideslip angle, deg
xviii
also: (1 M2)1=2
ight path angle, deg
dihedral, deg
a, e, r aileron, elevator, and rudder de ections, deg
x perturbation in x
perturbation in
t sampling time, s
vector of parameter update
downwash angle at the horizontal tail, deg
t wing twist angle, deg
span fraction
h ratio of horizontal tail to wing dynamic pressure
= ct=cr taper tatio
c=2 semi-chord sweep angle, deg
c=4 quarter-chord sweep angle, deg
LE leading edge sweep angle, deg
coe cient of viscosity for air
! angular velocity, rad=s
vector of unknown parameters
density of air, lb=ft3
side-wash angle, deg
xix
T tilt angle of the engines, deg
pitch angle, deg
state transition matrix
# vector of process noise distribution matrix, F
vector of system parameters
4. Superscripts
T transpose
1 inverse
predicted estimates
^ corrected estimates
5. Subscripts
a aileron
b base
cg center of gravity
cs cross sectional
e elevator
fus fusalage
h horizontal tail
i, j general indices
xx
k discrete-time index
LE leading edge
m measured variables
max maximum
M (at a given) mach number
p perturbed variables
plf planform
r rudder
ref reference, usually the wing, or a point on the wing
v vertical tail
w wing
wet wetted
0 initial condition
xxi
Chapter 1
Introduction
The thesis describes the system identi cation technique and its applications related
only to aircraft, but the system identi cation has been successfully applied to many diverse
elds. The system identi cation is basically related to modeling from experimental data and
can be applied to many di erent elds such as: economics, biology, materials, engineering,
automobiles, medicine and aerospace engineering.
1.1 System Identi cation - A De nition
System identi cation is a discipline of science which attempts to answer the inverse
problem. It tries to characterize a system in a suitable form based on observations of
the system?s behavior. The process of system identi cation involves certain fundamental
assumptions [1]:
1. The true state of the dynamic system is deterministic.
2. Physical principles underlying the dynamic process can be modeled.
3. It is possible to carry out speci c experiments.
4. Measurements of system inputs and outputs are available.
In the year 1962, Zadeh [2] gave a technical de nition for system identi cation as: \System
identi cation is the determination, on the basis of observation of input and output, of a
1
system within a speci ed class of systems to which the system under test is equivalent."
Several important features of system identi cation are highlighted by this description: the
necessity of obtaining input and output data, model characterization and selecting the best
model from the selected class. Ili [3] provided a philosophical de nition in contrast with
the above technical de nition as: \Given the answer, what are the questions, that is, look
at the results and try to gure out what situation caused those results." This de nition
points out the prime approach behind the model building process, namely that system
identi cation is an inverse problem.
1.1.1 Signi cance and Necessity of System Identi cation
The system identi cation tries to determine the best parameters for an assumed math-
ematical model, which is made up of di erential equations. The unknown parameters are
determined indirectly from measured data. Considering the aircraft system identi cation
problem, assume is an unknown parameter vector, which contains unknown stability and
control parameters of the dynamics equations, which are to be estimated. As a rst part
of system identi cation, a parameter estimation technique is used to estimate the unknown
parameter vector , such that the system response y ts the measured system response
z adequately. In aircraft system identi cation the measured system response is typically
ight test data. After the parameter estimation is performed, a second and more important
part, known as model validation needs to follow. It evaluates model reliability. If it seems
that the identi ed model does not agree with the speci c system standards and also the
2
statistical properties of the estimated parameters do not agree with the standards, then the
postulated model structure needs to be changed. Generally both steps must be repeated
again and again until the model is validated.
If the experimental model is already validated, then only the rst part, parameter
estimation, is dealt with. If the model structure is not known and adequate information
about the system structure is not available, the system identi cation technique is dealt with
as a whole, which includes both parameter estimation and model validation.
After mentioning the de nition, meaning and signi cance of system identi cation and
parameter estimation, the next question is, what is the necessity of system identi cation.
Why and where is system identi cation needed? This question must be answered due to
the two following examples of questions that have been raised about the utility of system
identi cation.
1. During one of the early American Institute of Aeronautics and Astronautics (AIAA)
conferences, someone raised a question: \When the aircraft is past the design and pro-
duction stage, and already ying, why and where do you need system identi cation?"[4]
2. An anonymous paper reviewers comment: \After perhaps 30 years of journal articles
and conferences on system identi cation, hardly any results have been of engineering
utility. Or if they are in use, the fact is not widely published."[4]
In ight vehicle development, system identi cation is a necessary step because it leads to
adequately accurate and validated mathematical models of ight vehicles, which are required
to [1],
3
1. understand the cause-e ect relationship that underlines a physical phenomenon,
2. investigate system performance and characteristics,
3. verify wind-tunnel and analytical predictions,
4. develop high- delity aerodynamic databases for ight simulators meeting FAA delity
requirements,
5. support ight envelope expansion during prototype testing,
6. derive high- delity and high-bandwidth models for in- ight simulators,
7. design ight control laws including stability augmentation systems,
8. reconstruct the ight path trajectory, including wind estimation and incidence anal-
ysis,
9. perform fault-diagnosis and adaptive control or recon guration, and
10. analyze handling qualities speci cation compliance.
Due to the advantages and the uses given above the system identi cation is made an essential
part of aerospace system design and development.
1.1.2 Objectives of Research
The aim of this thesis and research work, is to contribute to the eld of system iden-
ti cation of aircraft which has advanced over the recent decades, concentrating speci cally
4
on nonlinear systems and time domain methodologies. E cient and e ective system iden-
ti cation of aircraft is only possible if the applied approach is well coordinated and follows
speci c research guidelines. The thesis uses the lter error approach for all three examples
and tries to accommodate full details of such an e cient approach, summarizing the gen-
eral important concepts, techniques, computational procedures and three selected examples,
which were investigated based on the method.
1.2 Problem Description
The objective of system identi cation is to estimate the values of some unknown pa-
rameters in the system equations, which best represent the actual aircraft response, given a
set of ight time histories of an aircraft?s response variables. Here, the unknown parameters
are stability and control parameters.
The typical mathematical approach to this problem is to minimize the di erence be-
tween the ight response and the response computed from the system equations. This
di erence can be de ned for each response variable as the integral of the error squared.
Then the signal errors can be multiplied by weighting factors and summed to obtain the
response error which de nes an integral squared error criterion.
A mathematically more precise formulation can be made in probabilistic terms. A
probability that the estimated aircraft response time histories take on the values actually
observed can be de ned for each possible estimate of unknown parameters. These parameter
5
estimates should be chosen so that this probability is maximized. This is known as a
maximum likelihood formulation of the problem.
To describe the maximum likelihood estimator mathematically, it is necessary to de ne
the equations of motion for the aircraft system. The following steps should be followed.
1. Derive the dynamical equations of motion.
2. Predict the unknown stability and control parameters.
3. De ne the control inputs.
4. Express integral squared error criterion (to nd a vector of unknowns), which mini-
mizes the cost function (performance index) J.
5. Estimate the state error covariance matrix to construct the likelihood function.
6. Minimize the cost function. The modi ed Newton-Raphson method seems most suit-
able for aircraft derivative determination, both in terms of computer time and con-
vergence properties.
7. Include prior information. Information from wind tunnel tests, previous real-time
ight tests and simulated ight test results are often available with the values of
some of the aircraft derivatives. It may be desirable to include this information in
the algorithm. The use of this information is particularly important when there is a
linear dependence between the response and the unknown parameters.
6
1.3 Progress of Aircraft System identi cation
A brief survey of the various contributions to the system identi cation eld is presented
in this section, starting with applications in many elds. Applications to aircraft will follow.
Finally, the most modern techniques are briefed.
1.3.1 System Identi cation in General
The evolution from older system identi cation methods to statistically more accurate
approaches has been gradual. References can be found going back to the 18th century,
when Gauss [5] and Bernoulli [6] approached the system identi cation problem. The oldest
reference was found in the year 1777, when Bernoulli arrived at a solution by di erentiating
the likelihood function [6]. Bernoulli used the concepts of maximum likelihood and least-
squares solution, though he did not introduce these terms. In the year 1809, Gauss [5]
used the maximum likelihood principle. In the reference mentioned, he discussed the least-
squares method for orbit determination of the earth from astronomical measurements.
Following the fundamental methods of Bernoulli and Gauss during the 18th century,
the maximum likelihood estimator was rst discussed by Fisher [7] in the year 1912, as a
general statistical parameter estimator. Douglas [8] discussed \Inverse Problem Theorems"
in 1940. Feldbaum [9] considered the identi cation and control of the system as a single
problem in the \dual-control" theory; which was a little di erent than others, but still his
work was one of the most signi cant and also aimed at the present direction of investigation.
7
Research work in system identi cation theory prior to 1970, has been briefed in outstanding
research survey papers [10-13].
There are two kinds of system identi cation problems: deterministic (without state
noise) and non-deterministic (with state noise). Until the early forties, work on system
identi cation focus on the deterministic problem. But, succeeding the initial work in 1941-
1942 of Kolmogorov [14] and Wiener [15] the focus gradually shifted from deterministic
estimation to stochastic estimation. It was in the year 1960, when Kalman [16] followed basic
theories of Wiener and formulated a recursive solution to a ltering problem. This recursive
solution was directly executable through digital computations. Therefore, the Kalman lter
became popular quickly. Today it is the most widely used method for stochastic estimation.
The rst experimentation of the maximum likelihood estimator on a digital computer and
application to parameters estimation in an industrial plant was done by Astrom and Bohlin
in 1965 [17]. This was the beginning of the modern era of system identi cation methods.
1.3.2 System Identi cation of Aircraft
The aerodynamic modeling, which obtains relationships between the three forces X,
Y and Z along the three cartesian coordinates, and the moments L, M and N about these
axes as functions of linear translational motion variables u, v, w and rotational rates p, q,
r was introduced in the early 20th century [18, 19]. This was the initiation of the evolution
of aircraft system identi cation. After the pioneering work of Bryan [18], numerical, sta-
tistical and experimental research was started in the eld of real-time dynamics stability
8
and control. Glauert [20] initiated analyzing the phugoid motion of an aircraft in 1919 and
Norton [21, 22] presented papers in 1923 on estimation of stability and control derivatives.
The developments over the last nine decades have led to three di erent, but complementary
techniques for determining aerodynamic coe cients: 1) analytical methods, 2) wind-tunnel
methods, and 3) ight test methods [19].
The analytical methods and the wind-tunnel methods are performed to generate basic
information about the aerodynamic ight coe cients. Research in the computational uid
dynamics eld in recent years has positively a ected the analytical approaches by providing
numerical solutions of complete con gurations via sophisticated and advanced Euler and
Navier-Stokes ow solvers [19, 23, 24]. Experimental methods are important to validate the
analytical estimations. Wind-tunnel techniques, have in the past provided a huge amount
of data on innumerable ight vehicle con gurations and are, as a rule, a basis for any new
ight vehicle design [19]. These techniques are, however, often associated with certain lim-
itations of validity arising out of, for example, model scaling, Reynolds number, dynamic
derivatives, cross coupling, and aero servo elasticity e ects [19]. Determination of aerody-
namic derivatives from ight measurements is, therefore, important and necessary to reduce
limitations and uncertainties of the aforementioned two methods [19] and [25].
While detailing the chronological survey of research which led to the development and
universal acceptance of the maximum-likelihood estimation technique for aircraft stability
and control derivatives estimation, the more straightforward deterministic analysis will be
discussed rst, which will be followed by non-deterministic analysis. References [26] and [27]
9
discuss some of the investigations in estimation of unknown stability and control derivatives
of an aircraft from aircraft dynamic response data.
1.3.3 Deterministic and Non-Deterministic Analysis
1. Deterministic Analysis: System identi cation methods have become more ele-
gant and complex since Kalman?s contribution. The steady-state oscillator excitation
analysis by Breuhaus [28] and Milliken [29, 30] and also the pulse input method uti-
lizing Fourier analysis [31] are examples of the frequency response method, which
increased in popularity for aircraft analysis during the 1940s and 1950s. These meth-
ods produced the frequency response of the vehicle, but not the stability and control
derivatives of the di erential equations. Greenberg [32] and Shinbrot [33] attempted
to extract those parameters by de ning the values of the aircraft parameters that re-
sulted in the best t of the frequency response. Weighted least-square [33] and linear
least-squares [34] techniques were also applied to ight data in the 1960s. But these
techniques give poor results in the presence of measurement noise and yield biased
estimates. Also, they are time consuming and singularities are a serious problem with
them. The response curve- tting technique [35] was formulated in 1951, which is
equivalent to the output error technique. But due to the lack of e cient and faster
computational means it did not work at that time.
In the 1950s and 1960s, time-vector methods [36-40] were commonly applied graph-
ical procedures to determine aircraft parameters from ight data. However, those
10
methods give an incomplete set of parameters and only the responses of fairly simple
motion can be analyzed. In the early 1960s, before the invention of digital computers,
analog matching techniques [40, 41] which are time consuming and somewhat tedious,
were applied to ight data. Resulting estimates di ered depending on the skill and
knowledge of the person. Discussion of these past techniques concluded that a more
complete method of identi cation was needed, which can be faster and gives better
results.
Two separate articles [42, 43] on output error methods for obtaining aerodynamic
coe cients were published in the 1968. Taylor, Lawrence and Ili [42] discussed the
maximum-likelihood estimator for obtaining a complete set of aerodynamic parame-
ters from ight data. Larson and Fleck [43] described a quasi-linearization method for
parameter identi cation of an aircraft. This research had produced an excellent model
which su ciently described the resulting motion of the aircraft; which is the reason
for the success of these two methods. Due to these two studies on aircraft identi -
cation using nonlinear minimization techniques, the interest in analysis of ight data
was increased. After only a year in the 1969, Taylor and Ili [44] modi ed these two
techniques to include prior information. Minimization of this modi ed cost function
does not result in a maximum-likelihood estimator because it was based on the joint
probability distribution rather than the conditional probability [45]. Some excellent
and very successful computer algorithms on system identi cation have been published
11
[46-50]. The maximum likelihood estimator method was found very useful for ight
dynamic response [47-55].
Mehra [51] applied the Kalman lter to estimate aircraft aerodynamic parameters,
which produced poor results, because the state estimation and parameter estimation
were both biased. Hence, they did not converge to good results. Taylor, Powell and
Mehra [56] obtained better results after addition of the derivative of the state. In
the Netherlands, two successful applications of the Kalman lter, providing the state
estimation and aircraft parameters identi cation were obtained [57, 58].
2. Non-Deterministic Analysis: There are two types of techniques applied for pa-
rameter estimation with measurement and state noise: the Kalman lter technique
[51, 56-62] and maximum-likelihood estimator technique [52, 53, 63-65]. In the case of
non-deterministic analysis, the maximum likelihood estimator technique is normally
known as the lter error method. It was after 1965 that Astrom [59] and Kashyap
[60] described general applications of the extended Kalman lter. During 1970, for
the discrete-time system Taylor, Powell and Mehra [56] applied the extended Kalman
lter to simulated aircraft data with a state noise input. Chen and Eulrich [61] ex-
perimented with an application to that of Taylor, Powell and Mehra in 1971, but
it produced inconclusive results, as the state noise input was small, and also due to
nonlinearity of the system. Yazawa [62] produced very good results, when he applied
a simpli ed extended Kalman lter.
12
Ili [65] applied the maximum-likelihood technique to output data of an aircraft ying
in atmospheric turbulence. The results were similar for the same aircraft ying in
smooth air, which is without state noise.
1.3.4 Latest Techniques of System Identi cation
The maximum likelihood estimator technique is used in most of the latest techniques
of system identi cation. One of these methods, the lter error technique is completely
detailed in Chapter 3. The digital computers of the modern era have the speed to imple-
ment self-regulating and self-governing data processing capability, which has changed the
concentration of ight test data analysis utilizing methods based on frequency domain to
methods based on time domain. Using modern computers, relatively larger numbers of
aircraft stability and control derivatives can be estimated from only a single aircraft test.
A coordinated approach based on ight test instrumentation, ight test techniques, and
methods of data analysis gradually evolved for system identi cation as applied to aircraft
[19, 46]:
1. Instrumentation and lters, which cover the entire ight data acquisition process
including adequate instrumentation and airborne or ground-based digital recording
equipment. E ects of all kinds of data quality should be accounted.
2. Flight test techniques, which are related to the selected ight vehicle maneuvering
procedures. The input signals should be optimized in their spectral composition to
excite all response modes from which parameters are to be estimated.
13
Figure 1.1: Block Diagram for System Identi cation to aircraft [1]
3. Analysis of ight data, which includes the mathematical model of the ight vehicle
and an estimation criterion that devises a suitable computational algorithm to adjust
starting values or a priori estimates of the unknown parameters until a set of best
parameter estimates is obtained that minimizes the response error.
Synchronizing these interdependent subjects, four important aspects (Fig. 1.1) of the art
and science of system identi cation have to be carefully treated:
1. design of the control input shape to excite all modes of the vehicle dynamic motion
2. selection of instrumentation and lters for high accuracy measurements
3. type of ight vehicle under investigation to de ne the structure of a possible mathe-
matical model
14
4. the quality of data analysis by selecting the most suitable time or frequency domain
identi cation method
The above four points must be followed carefully while investigating each aircraft. If they
are followed, modern system identi cation methods for aircraft can generally be expected
to provide good results. The modern techniques of data analysis for parameter estimation
are classi ed as: 1) Equation error technique, 2) Output error technique, 3) Filter error
technique. The last technique has been chosen for this study and is discussed in Chapter 3.
System identi cation is a multidisciplinary technique, which covers the elds of control
theory, numerical techniques, theories of statistics, sensors and instrumentation, ight test
techniques, signal processing, and ight dynamics. Basic knowledge of each of these elds
enables the user to generate e cient results.
15
Chapter 2
Background and Theoretical Development
Dynamical systems can generally be characterized by di erential equations, whose or-
der is dependent on the complexity of the whole process. The system identi cation process
starts with postulating a model, so that it becomes possible to estimate the parameters of
state through simulation. Usually this is done by solving an initial value problem utilizing
numerical integration techniques. For system identi cation of aircraft, state space dynam-
ical models characterizing the aircraft motion are usually used. These models are derived
from Newtonian mechanics; which gives well formulated kinematic equations of aircraft
motion with translational and rotational degrees of freedom [66].
2.1 Dynamical Model Description
As stated earlier, the system identi cation process for aircraft starts with postulating
a dynamical model and governing equations of aircraft motion. The equations of motion
also include the output and measurement equations. The complete set of equations is then
[66]:
_x=f[x(t);u(t); ]+F(#)w(t); x(t0)=x0
y=g[x(t);u(t); ]
z(tk)=y(tk)+Gv(tk) (2.1)
16
where x is the (n 1) state vector, u is the (s 1) input vector, y is the (m 1) observation
vector and z is the (m 1) measurement vector. The unknown parameter vector is given
as:
=[ T #T x0T]T (2.2)
In this speci c case of system identi cation to aircraft, the system parameter vector is
made up of aircraft stability and control derivatives. The above nonlinear dynamical model
can be linearized if required, using numerical approximation.
2.2 Cost Function (Performance Index) J
Using the measurement vector z(t), observation vector y(t), measurement noise covari-
ance matrix R (m m), for ny observation variables and N data points, a cost function, or
performance index, J can be de ned as follows [67].
J( ;R) = 12
Z N
1
[z(t) y(t)]T R 1 [z(t) y(t)]dt+ N2 ln [det(R)] + Nny2 ln det(2 ) (2.3)
where is a column vector of parameters, which are to be estimated.
2.3 Essentials of Cost Function Minimization
For any experiment, the number of observations ny and number of data points N shall
be xed. So, the last term in the above equation is a constant and hence can be neglected
17
without a ecting the results [1]. Now for the measurement noise covariance matrix R, there
are only two possibilities.
1. Known Measurement Noise Covariance Matrix R: Based on previous ight tests or
using wind tunnel data, the measurement noise covariance matrix R can be assumed
for this experiment. In this case, the second term in the above equation is constant.
Then the cost function is simpli ed as [68].
J( ) = 12
Z N
1
[z(t) y(t)]T R 1 [z(t) y(t)]dt (2.4)
Therefore, the cost function is easy to calculate and also quadratic in nature. It can
be reduced using modern optimization methods given in the next section.
2. Unknown Measurement Noise Covariance Matrix R: Jategaonkar [1] used a relaxation
strategy in which optimization of the likelihood function, Eq. 2.3 was carried out in
two steps. In the rst step, di erentiate Eq. 2.3 partially with respect to R, set it
equal to zero and derive
R = 1N
Z N
1
[z(t) y(t)] [z(t) y(t)]T dt (2.5)
Substituting the value of R, obtained in Eq. 2.5 in Eq. 2.3 will give the result:
J( ) = 12nyN + N2 ln [det(R)] + Nny2 ln det(2 ) (2.6)
18
The rst term and the last term on the right side of the above equation are con-
stant because N and ny are xed for any model. Therefore, with out a ecting the
minimization results the cost function reduces to [1]:
J( ) = det(R) (2.7)
Now, the value of parameter vector is be determined, which should minimize det(R), that
is the cost function. The elements of may be determined as per the requirement by a
variety of optimization techniques. A relaxation strategy will be followed in this thesis, as
outlined below.
1. Assume parameter vector with suitable initial values.
2. Determine system observation (or output) matrix Y. Compute the innovation vector
(Z - Y). Determine measurement noise covariance matrix R.
3. Calculate cost function J.
4. Apply any of the optimization techniques to reduce Cost function J.
5. Repeat the algorithm from step 2 to step 4 until the cost function reduces satisfactorily
and also, the estimates converge to a speci c and constant values.
19
Numerous optimization techniques have been developed for linear and nonlinear optimiza-
tion. Direct search techniques, accelerated gradient based techniques and modi ed Newton-
Raphson techniques [69] are few of them. In this thesis the modi ed Newton-Raphson
technique is utilized.
2.4 Modi ed Newton-Raphson Algorithm
The Newton-Raphson technique is an iterative method, which nds a zero of the gra-
dient of a cost function, that is
@J( )
@ = 0 (2.8)
@J
@ is expanded using two-term Taylor?s series about the i
th value of the parameter vector
:
@J( )
@
i+1
@J( )
@
i
+
@2J( )
@ 2
!
i
(2.9)
where = i+1 i is the parameter change and (@2J( )@ 2 )i is the second derivative of the
cost function with respect to at the ith iteration [68]. The change in can be written
=
"
@2J
@ 2
!
i
# 1
@J
@
i
(2.10)
The Newton-Raphson algorithm described above, is more e cient than the gradient method
because it attempts to predict where the local minimum point is quadratically convergent
[69, 70]. To compute the rst gradient of the cost function, Eq. 2.3 can be partially
20
di erentiated with respect to and the result is:
@J
@ =
Z N
1
@y(t)
@
T
R 1 [z(t) y(t)]dt (2.11)
Di erentiating the above equation with respect to , the second gradient matrix can be
formed as [68]:
@2J
@ 2 =
Z N
1
@y(t)
@
T
R 1
@y(t)
@ dt
+
Z N
1
"
@2y(t)
@ 2
#T
R 1 [z(t) y(t)]dt (2.12)
The rst gradient of the cost function can be calculated easily but computation of the
second gradient is complicated and time consuming. The prime reason for this complexity
is the computation of the second gradient of the response @2y@ 2 , which is found in the second
integral of Eq. 2.12. Balakrishnan [71] suggested a simpli cation as follows. The second
term of Eq. 2.12 goes to zero as the process converges. Hence, it can be neglected and the
second gradient can be approximated as:
@2J
@ 2 =
Z N
1
@y(t)
@
T
R 1
@y(t)
@ dt
(2.13)
This yields a more exible technique called the modi ed Newton-Raphson method. This
technique is also referred to as the Newton-Balakrishnan method or Gauss-Newton method.
21
This leads to a system of linear equations that can be summarized as follow [68]:
Denote P1 =
"
@2J
@ 2
!
i
# 1
and P2 =
@J
@
i
(2.14)
From Eq. 2.11 and Eq. 2.13:
P1 =
Z N
1
@y(t)
@
T
R 1
@y(t)
@ dt
(2.15)
P2 =
Z N
1
@y(t)
@
T
R 1 [z(t) y(t)]dt (2.16)
Then,
i+1 = i + ; where = P1 1P2 (2.17)
22
Chapter 3
System Identification of General Aviation Aircraft
The various aircraft parameter estimation techniques are:
1. The equation error technique
2. The output error technique
3. The lter error techniques
4. The frequency domain techniques
5. The recursive parameter estimation technique
6. The arti cial neural network based techniques
A lter error technique is used in this thesis. The lter error techniques represent the general
stochastic approach to aircraft system identi cation introduced by Balakrishnan [63]. Ili
[65] and Mehra [51, 56, 72] utilized these techniques in the 1970s. These techniques are
highly e cient while working with measurement and process noise.
3.1 Filter Error Technique for Linear Systems
A block diagram of the lter error technique is shown in Fig. 3.1. In this section, the
lter error technique for linear systems will be discussed with theoretical advancements and
computational aspects.
23
Figure 3.1: Block Diagram for Filter Error Technique [68]
A dynamic system can be described by the following stochastic linear mathematical
model [1].
_x(t) = Ax(t) +Bu(t) +Fw(t) +bx; x(t0) = 0 (3.1)
y(t) = Cx(t) +Du(t) +by (3.2)
z(tk) = y(tk) +Gvk; k = 1;2;:::N (3.3)
where x is a state vector (n 1), u is a control input vector (s 1), y is an output vector
(m 1), and z is a measurement vector (m 1) sampled at N discrete time points. Matrices A
(n n), and C (m n), contain unknown stability derivatives and matrices B (n s), and D
(m s), contain unknown control derivatives. Matrices F (n n), and G (m m), represent
24
the process and measurement noise distribution matrices, respectively. In the state and
observation equations, bx and by are bias terms. In the most general case all of the elements
of the matrices A, B, C, D and bias terms, bx and by form the parameter vector , which
is to be estimated.
From the Eqs. 2.3-2.4, output vector, y, is required now, to calculate the cost function
J and the measurement noise covariance matrix, R. A Kalman lter is used to estimate the
output vector given the observed data, Ztk. Parameter estimates for and R are obtained
by minimizing the cost function.
Though it can be assumed that the process and measurement noises are additive, the
estimation process is still complicated because the nature of the system is non-deterministic
[67]. Because of the process and measurement noise, it is not possible to integrate the state
equations. Instead, a standard Kalman lter is utilized as an optimal state estimator for
the linear system [73, 74].
3.1.1 Kalman Filter
The Kalman lter consists of two steps, a prediction and a correction step.
1. Prediction step: For the linear model postulated in Eqs. 3.1-3.2, the prediction step
is given by [1]:
ex(tk+1) = bx(tk) + Bu(tk) + bx (3.4)
y(tk) = Cex(tk) +Du(tk) +bx (3.5)
25
2. Correction step
bx(tk) = ex(tk) +K[z(tk) y(tk)] (3.6)
Here, ex is a predicted state vector, bx is a corrected state vector, u is an average of the
control input at two successive discrete time steps, ey is a predicted observation vector and
[z(tk) ey(tk)] is the residual (or innovation). and are the state transition matrix and
integral of state transition matrix (derived in appendix A). K is the Kalman gain (n m)
and is dependent on the observation matrix C, state prediction error covariance matrix P
(n n), and measurement error covariance matrix R (m m), and given as [67]:
K = PCTR 1 (3.7)
The measurement noise covariance matrix R can be found using Eq. 2.5 and the state
prediction error covariance matrix P is calculated by solving the Riccati equation.
3.1.2 Solution of Riccati Equation
The Riccati equation is solved to compute the state prediction error covariance matrix,
P. A steady-state form of the continuous-time Riccati equation is used here because it is
easier than the discrete-time equation [1, 67, 75]. The rst order projection of the Riccati
equation is given below.
AP +PAT 1 tPCTR 1CP +FFT = 0 (3.8)
26
Eq. 3.8 can be solved using Potter?s method [1, 76] based on eigenvector decomposition.
Only the necessary details are provided here. Solution of the Riccati equation follows three
steps and is very straightforward.
1. Start with formulation of the Hamiltonian matrix, H, as [67]:
H =
"
A j FFT
1 tCTR 1C j AT
#
(3.9)
2. Compute the eigenvalues and eigenvector of the Hamiltonian matrix, H. Divide the
eigenvector matrix into four equal size matrices such that the eigenvectors related to
eigenvalues with positive real parts remain on the left side. Controllability due to
process noise and observability ensure that exactly one half of the eigenvectors will
have positive real parts (i.e., unstable eigenvalues). Suppose M is the eigenvector
matrix of H, then it can be divided and arranged as [67]:
M =
M
11jM12
M21jM22
(3.10)
3. The solution of the steady-state, continuous-time Riccati equation is as follows [67]:
P = M11M21 1 (3.11)
The state prediction error covariance matrix, P, is a steady-state matrix and is not depen-
dent on time [1].
27
3.1.3 Formulations for Process Noise
Three principle formulations have been used to account for both the process and mea-
surement noises based on the way the measurement noise covariance matrix, R; Kalman
gain matrix, K; and the noise distribution matrices, F and G, are estimated [54, 67, 68].
1. Natural Formulation
The natural formulation [67, 68] utilizes the unknown parameter vector in the best
possible optimization method and minimizes the cost function. The elements of parameter
vector are all or some elements of matrices the A, B, C, D, F, G. This formulation is
very easy to de ne and theoretically a possible and perfect technique. But, if the problem
considered has process and measurement noises, this formulation is not a practical one.
Mainly because this formulation does not have an explicit solution to estimate the noises.
There are three primary disadvantages of this formulation, primarily due to the complica-
tions in the estimation of F and G matrices. The drawbacks are convergence, singularity
and computational burden. Therefore, this technique is not used extensively these days.
2. Innovation Formulation
The innovation formulation was suggested to overcome the di culties of the natural
formulation. It is very similar to the formulation which was followed in the output error
method [51, 62, 67]. As it was mentioned in the previous section, the prime di culty
of the natural formulation is estimation of the F and G matrices. From the de nitions
of the cost function, J, in Eq. 2.3 and measurement noise covariance matrix, R, in Eq.
2.5, it is clear that they are dependent on the F and G matrices indirectly through the
28
Kalman gain. Therefore, if the Kalman gain is directly estimated, then the complications
due to the natural formulation would be solved. In this case, the elements of the parameter
vector are all or some elements of the matrices A, B, C, D, K. This formulation
works well with process and measurement noises and also for linear and nonlinear systems.
The straightforward solution of the cost function, J, from Eq. 2.6 and measurement noise
covariance matrix, R, in Eq. 2.5 solves the convergence problems.
Because the innovation formulation estimates the Kalman gain K directly (in place
of deriving the same from Eq. 3.7, it solves the remaining parameter estimation problem
utilizing only three equations (3.4-3.6). The best part of this formulation is it excludes
the need to compute the most complex part of the algorithm, that is the solution of the
Riccati equation. Thus, this technique looks more appealing computationally, however, it
has a few disadvantages [67]. The size of matrix K is larger than that of matrix F. The
parameter vector is larger here, which takes more computational time. The elements
of K do not have direct physical meaning [68]. Moreover, the Kalman gain, K as well as
system parameters and measurement noise covariance matrix, R, are computed separately,
which may give unsuitable estimates for R and K.
3. Combined Formulation
Maine and Ili [55, 67] proposed a Combined Formulation in 1981, which takes the best
features from both of the previous formulations. The combined formulation estimates the
process noise matrix, F, and it calculates the measurement noise covariance matrix, R, from
Eq. 2.5. Hence, the parameter vector, , consists of all or some elements of the matrices A,
29
B, C, D, F, and a two-step relaxation strategy is applied. Here the rst step is to estimate
the measurement noise covariance matrix, R, which is relatively simpler. The Parameter
vector, , is estimated in the second step utilizing the modi ed Newton-Raphson technique.
The primary advantage of this method is that the process noise matrix, F, has a physical
meaning. To estimate the state prediction error covariance matrix, P, the solution of the
Riccati equation is necessary. Although this solution is tedious, it is worthwhile because
the results of this formulation are more practical with respect to convergence, singularities,
parameter estimates and computation burden [1].
3.1.4 The Filter Error Algorithm for Linear Systems
After comparing the three formulations, only the combined formulation is followed in
this thesis due to its inherent e ciency. The lter error algorithm for linear systems is
outlined in Fig. 3.2. In this technique the following tasks are important.
1. Select suitable initial values for the unknown parameter vector, , the Kalman gain,
K, the state prediction error covariance matrix, P, and states x. Load the available
ight data (from real time ight or from a simulator).
2. Apply the Kalman lter and calculate the observation vector, Y; innovation matrix,
initial measurement noise covariance matrix, R; and based on that the initial cost
function, J.
3. Solve the Riccati equation and nd the new Kalman gain. Calculate the new mea-
surement noise covariance matrix, R, and based on that the new cost function, J.
30
Figure 3.2: Algorithm of Filter Error Technique for Linear Systems [1]
31
4. Apply the modi ed Newton-Raphson technique to minimize the cost function.
5. Update the parameter vector .
6. Repeat the steps 3 - 5 until all parameters converge to speci c and constant values.
3.1.5 Parameter Update
The modi ed Newton-Raphson method is applied in this section to update the param-
eter vector, , such that it minimizes the cost function after each iteration. The following
equations were derived in section 2.4 [1].
i+1 = i + ; where = P1 1P2 (3.12)
P1 =
Z N
1
@y(t
k)
@
T
R 1
@y(t
k)
@ dt
(3.13)
P2 =
Z N
1
@y(t
k)
@
T
R 1 [z(tk) y(tk)]dt (3.14)
In the above equations, everything is readily available to use, except the sensitivity matrix
or the gradient matrix of response, given below [1]:
@y
@
ij
= @yi@
j
(3.15)
The sensitivity matrix (@y@ ) can be solved using several ways. However, the classical method
is used in this thesis by solving the sensitivity equations [69]. The nite di erence method
is also an excellent method, which is usually applied for nonlinear system. The required
32
sensitivity equations are obtained by di erentiating the system equations with respect to
elements of the parameter matrix, . Partial di erentiation of Eq. 3.4 results in the
following equation [1].
@ex(tk+1)
@ =
@bx(tk)
@ +
@
@ bx(tk) +
@B
@ u(tk) +
@
@ Bu(tk) +
@bx
@ +
@
@ bx (3.16)
The above equation is the exact gradient of Eq. 3.4 and is necessary to compute @ @ and
@
@ . Maine and Ili [67] eliminated this computation burden by approximating the gradient
computation signi cantly, explained below:
@ex(tk+1)
@
@bx(tk)
@ +
@A
@ x(tk) +
@B
@ u(tk) +
@bx
@ ;
@ex(1)
@ = 0 (3.17)
Two new terms are introduced in the above equations. An average of two consecutive
inputs is shown as u and x is an average of states at two successive time steps. Partial
di erentiation of Eq. 3.5 and Eq. 3.6 result in the following equation [67].
@y(tk)
@ = C
@ex(tk)
@ +
@C
@ ex(tk) +
@D
@ u(t) +
@by
@ (3.18)
@bx(tk)
@ =
@ex(tk)
@ K
@y(tk)
@ +
@K
@ [z(tk) y(tk)] (3.19)
Eqs. 3.17-3.19 represent a set of sensitivity equations. Except for the @K@ all terms are
available at this point. Therefore, the gradient of the Kalman gain matrix is computed to
complete this optimization problem. Di erentiate Eq. 3.7 to derive the following result
33
[67].
@K
@ =
@P
@ C
TR 1 +P@CT
@ R
1 (3.20)
For each speci c time step and at this stage in the algorithm measurement noise covariance
matrix R, system matrix C and state prediction error covariance matrix P are xed and
constant. Partial di erentiation of each element matrixC with respect to related parameters
is one and elsewhere it is zero. Thus, only @P@ needs to be calculated to compute the gradient
of K. Note that @P@ is a three dimensional matrix and hence @K@ is also a three dimensional
matrix. @P@ is derived by di erentiating the steady-state Riccati equation 3.8 [67].
A@P@ + @A@ P +P@A
T
@ +
@P
@ A
T 1
t
@P
@ C
TR 1CP 1
tP
@CT
@ R
1CP
1 tPCTR 1@C@ 1 tPCTR 1C@P@ +F@F
T
@ +
@F
@ F
T = 0 (3.21)
Arrange the left-hand side terms and simplify mathematically using matrix properties like
(A+B)T = AT +BT and (AB)T = BTAT. The result is shown here [67].
A@P@ + @P@ AT = C +CT (3.22)
where
A = A 1 tPCTR 1C = A 1 tKC (3.23)
and
C = @A@ P + 1 tPCTR 1@C@ P @F@ FT (3.24)
34
There is one equation like Eq. 3.22 for each element of the system parameter vector .
This leads to a set of Lyapunov equations of the form AX + XAT = B [1, 67]. Since the
A matrix remains similar for the entire set, the Lyapunov equations are solved e ciently
utilizing the following transformation [67].
A0 = T 1AT (3.25)
and
(C +CT)0 = T 1(C +CT)T 1 (3.26)
Here, T is the matrix of the eigenvectors of A. Due to this alteration A0 is now diagonal
and thus Eq. 3.22 is left as [67]:
A0
@P
@
0
+
@P
@
0
A0 = (C +CT)0 (3.27)
As the A0 is diagonal, the above equation can be solved for
@P
@
0
easily. Generally for a
diagonal matrix A, AX +XAT = B can be solved by [67]
Xij = Bij(A
ii +Ajj)
(3.28)
The gradient of the state prediction error covariance matrix P is calculated using back
transformation [67].
@P
@ = T
@P
@
0
T (3.29)
35
Now, Eqs. 3.12-3.14 can be solved easily. As stated earlier, this technique leads to an
unconditional minimum of the cost function. All system parameters are treated separately.
The results obtained using this technique will be provided and discussed in the next chapter.
3.2 Filter Error Technique for Non-linear Systems
The general mathematical model of a non-linear system is given below:
_x(t) = f[x(t);u(t); ] +Fw(t); x(t0) = x0 (3.30)
y(t) = g[x(t);u(t); ] (3.31)
z(tk) = y(t) +Gv(tk) (3.32)
Here f and g are nonlinear functions, x(t) is a state vector and u(t) is an input vector. F
and G are time-invariant noise distribution matrices.
The lter error technique for nonlinear systems is illustrated in Fig. 3.3. For non-linear
systems, response gradients are calculated based on nite di erence approximations. An
extended Kalman lter [77, 78] dependent on the rst-order approximation of the dynamical
equations is used for non-linear ltering.
36
Figure 3.3: Algorithm of Filter Error Technique for Nonlinear Systems [1]
37
3.2.1 Steady State Filter
The combined formulation is followed, for the nonlinear case just as was done for the
linear case. The two-step relation strategy, discussed in section 2.4 is applied to estimate
R. While doing this, the following tasks are important.
1. Select some suitable initial values for the unknown parameter vector . Load the
available ight data (from real time ight or from simulator).
2. Apply the extended Kalman lter and calculate the observation vector, Y; innovation
matrix; initial measurement noise covariance matrix, R; and the initial cost function,
J. Numerically integrate the state equations. For numerical integration, a fourth order
Runge-Kutta method was used. A description of this method is given in appendix B.
3. Solve the Riccati equation and nd the new Kalman gain. Calculate a new measure-
ment noise covariance matrix, R, and new cost function, J.
4. Apply the modi ed Newton-Raphson technique to minimize the cost function.
5. Update the parameter vector, .
6. Repeat steps 3 - 5 until all parameters converge to constant values.
The optimal lters are practically impossible for non-linear systems. Hence, an extended
Kalman lter is used here. A two-step, non-linear, constant-gain lter is formulated as [1].
38
Prediction step
ex(tk+1) = bx(tk) +
Z tk+1
tk
f[x(t);u(t); ]dt; bx(t0) = x0 (3.33)
y(tk) = g[ex(t);u(t); ] (3.34)
Correction step
bx(tk) = ex(tk) +K[z(tk) y(tk)] (3.35)
The Kalman gain K is dependent on the measurement noise covariance matrix, R, state
prediction error covariance matrix, P; and system matrix, C. It is expressed as [1]:
K = PCTR 1 (3.36)
where
C =
@g[x(t);u(t); ]
@x
t=t0
(3.37)
A steady-state, continuous-time, Riccati equation is solved here to calculate P. The lin-
earized matrices A and C are used, and solved using a nite di erence method which is
explained below.
Finite Di erence Technique:
39
If each of the state variables is perturbed by xj, the elements of the matrices A and
C can be approximated using a central di erence formula as [1]:
Aij
fi[x+ xj;u; ] fi[x xj;u; ]
2 xj
!
x=x0
(3.38)
Cij
gi[x+ xj;u; ] gi[x xj;u; ]
2 xj
!
x=x0
(3.39)
For the non-linear system, the response gradients are approximated using the nite di er-
ence approximation. For a small perturbation j in each variable of the unknown param-
eter vector , the perturbed response variable yp for each of the unperturbed variable yi is
computed. The sensitivity response gradient is approximated as
@y(t
k)
@
= yp(tk) yi(tk) (3.40)
For each element of the unknown parameter vector, , the state and observation vari-
ables are calculated utilizing a two step non-linear steady-state lter as given below [1].
Prediction step
exp j(tk+1) = bxp j(tk) +
Z tk+1
tk
f[xp j(tk+1);u(tk+1); + jej]dt (3.41)
yp j(tk) = g[exp j(tk+1);u(tk+1); + jej] (3.42)
40
Correction step
bxp j(tk) = exp j(tk) +Kp j[z(tk) yp j(tk)] (3.43)
The predicted state variables ex in Eq. 3.41 are computed by numerical integration. To
compute Eq. 3.43 the perturbed Kalman gain matrix Kp j needs to be calculated, which
is given by [1]
Kp j = Pp jCTp jR 1 (3.44)
here, the perturbed system matrices Ap j and Cp j need to be calculated to use them to
solve the Riccati equation, which are approximated using a central di erence method.
3.2.2 Parameter Update
Referring to section 3.1.5 and Eqs. 3.13-3.14, all quantities can be derived for non-
linear systems. The parameter vector will be updated using Eq. 3.12. This technique
also leads to an unconditional minimum of the cost function. The results obtained using
this technique will be provided and discussed in the next chapter.
41
Chapter 4
Prediction of Aircraft Stability and Control Derivatives
Techniques to calculate the stability and control derivatives analytically, are presented
in this chapter. These derivatives are an important part of the experiment because they
are required as an input to determine dynamic stability and control behavior [79]. These
methods are referred from [79-83]. Only preliminary information is presented here. Details
should be found from [79-83]. Stability and control derivatives given in Eqs. 5.1-5.5 are
predicted using the following methods. Initial values derived using these methods were used
in the experimental aircraft Trainer.
The important notes while using these methods are [79]:
1. The methods presented in this chapter apply only to rigid aircraft.
2. The methods presented in this chapter apply only to subsonic speed regimes.
3. All derivatives are in rad 1.
4. All coe cients and derivatives are de ned in the stability axes system.
5. The drag prediction methods apply only in the ight cases where the boundary layer
is turbulent.
6. The drag prediction methods apply only to smooth surfaces.
42
4.1 Steady State Coe cients
The de nitions of steady state coe cients CD, CL and Cm0 are given below. They
were used to determine the lift, drag and pitching moment steady state coe cients [79].
CL is the aircraft?s steady-state lift coe cient, given as [79]:
CL = W=(qS) (4.1)
CD is the aircraft?s steady-state drag coe cient. It is calculated using the following
method. As shown in Fig. 4.1, it is related to a speci c value of steady-state lift coe cient
CL. The steady-state aircraft drag coe cient is given as follows [79]:
CD = CDw +CDfus (4.2)
The wing drag coe cient is predicted from [79]:
CDw = CD0w +CDLw
where CD0w = (Rwf)(RLS)(Rfw)f1 +L0(t=c) + 100(t=c)4g(Swet=S)
and CDLw = (CLw)2= Ae+ 2 CLw tv + 4 2( t)2w (4.3)
43
Figure 4.1: Determination of Drag Coe cient from a known Lift Coe cient
The fuselage drag coe cient is predicted from [79]:
CDfus = CD0fus +CDLfus
where CD0fus = (Rwf)(Cffus)f1 + 60(lf=df)3 + 0:0025(lf=df)g(Swet=S) +CDbfus
and CDLfus = 2 2(Sb=S) + cdc 3(Splffus=S) (4.4)
Cm0 is the aircraft zero lift pitching moment coe cient and can be predicted using the
following equation [79]:
Cm0 = Cm0wf +Cm0h (4.5)
44
Cm0wf is the zero lift pitching moment coe cient due to wing-fuselage combination. It is
predicted from [79]:
Cm0wf =f(Cm0w)+(Cm0f )gf(Cm0)M+(Cm0)M=0g
where Cm0w =f(Acos2 c=4)=(A+2cos c=4)g(cm0r+cm0t)=2
+( Cm0= t) t
and Cm0f =f(k2 k1)=36:5Scg[Xf(wfi2)(iw+ 0Lw+iclf) xig] (4.6)
Cm0h is the zero lift pitching moment coe cient because of the horizontal tail. It is predicted
from [79]:
Cm0h = (xach xref)CL0h (4.7)
4.2 Stability Derivatives
The techniques for calculating the following stability derivatives are addressed in this
section:
1. Aircraft speed derivatives.
2. Angle of attack derivatives.
3. Angle of sideslip derivatives.
4. Roll rate derivatives.
45
5. Pitch rate derivatives.
6. Yaw rate derivatives.
As it was stated before, only preliminary equations are addressed here. Complete analysis
can be found in [79-83].
4.2.1 Aircraft Speed Derivatives
The aerodynamic speed derivative CLV , lift with respect to aircraft speed is determined
using the following equation [79].
CLV =fM2(cos c=4)2CLg=f1 M2(cos c=4)2g (4.8)
The aerodynamic speed derivative CDV , drag with respect to aircraft speed is deter-
mined using the following equation [79].
CDV = M(@CD=@M) (4.9)
The aerodynamic speed derivative CmV , pitching moment with respect to aircraft speed
is determined using the following equation [79].
CmV = CL(@XacA=@M) (4.10)
46
4.2.2 Angle of Attack Derivatives
The derivative of the aerodynamic lift with respect to aircraft angle of attack (also
called as the aircraft lift curve slope) can be determined from [79]:
CL = CL wf +CL h h(Sh=S)(1 d =d ) (4.11)
CL wf is the wing-fuselage lift curve slope, determined as [79]:
CL wf = KwfCL w (4.12)
where CL w is the wing lift curve slope and is given by [79]:
CL w = 2 A=[2 +f(A2 2=k2)(1 + tan2 c=2= 2) + 4g1=2] (4.13)
CL h is the horizontal tail lift curve slope. It can be calculated using the above equation
with suitable substitution of horizontal tail parameters for the wing parameters. All other
quantities in Eqs. 4.11-4.13 are de ned in [79].
The derivative of the aerodynamic drag with respect to aircraft angle of attack can be
computed as follows [79].
CD = (@CD=@CL)CL (4.14)
47
The derivative of the aerodynamic pitching moment with respect to aircraft angle of
attack can be calculated as follows [79].
Cm = (@Cm=@CL)CL (4.15)
4.2.3 Angle of Sideslip Derivatives
The aircraft equations of motion contain the aerodynamic pitch velocity derivatives
Lv, Nv, Yv. They are derivatives of rolling moment, yawing moment and side force with
respect to pitch velocity. If the angle of sideslip is very small then it is given by:
= v=u (4.16)
From the above equation:
L = Lvu
N = Nvu
Y = Yvu (4.17)
48
L , N and Y are calculated rst to predict Lv, Nv and Yv and given as [83],
L = qSbCl I
x
N = qSbCn I
z
Y = qSCy m (4.18)
The theoretical methods to predict the coe cients given in the above equations are
described in this section.
The coe cient of the rolling moment derivative with respect to the sideslip, Cl (also
called the dihedral e ect) can be calculated from [79]:
Cl = Cl wf +Cl h +Cl v (4.19)
The wing-fuselage contribution is given by [79]:
Cl wf = 180 [CLwff(Cl =CL) c=2(KM )(Kf)+(Cl =CL)Ag
+ f(Cl = )KM +( Cl = )g+( Cl )zw
+( ttan c=4)f( Cl )=( ttan c=4)g] (4.20)
The horizontal and vertical tail contributions are given by [79]:
Cl h = (Cl hf )(Shbh=Sb)
49
Cl v = (Cy v)f(zvcos lvsin )=bg (4.21)
The coe cient of the yawing moment derivative with respect to the sideslip, Cn (also
called static dihedral stability) is computed from [79]:
Cn = Cn w +Cn f +Cn v (4.22)
The wing, fuselage and vertical tail contributions are given in order as [79]:
Cn w = 0
Cn f = 180pi KNKR1(Sfsbf=Sb)
Cn v = (Cy v)f(lvcos +zvsin )=bg (4.23)
The coe cient of the side force derivative with respect to the sideslip, Cy is calculated
from [79]:
Cy = Cy w +Cy f +Cy v (4.24)
The wing, fuselage and vertical tail contributions are given in order as [79]:
Cy w = 0:00573(j j)
Cy f = 2Ki(Scs=S)
50
Cy v = kv(CL v)(1+d =d ) v(Sv=S) (4.25)
4.2.4 Roll Rate Derivatives
The rolling moment with respect to the roll rate is [83],
Lp = qSb
2Cl
p
2Ixu0
The coe cient of the rolling moment with respect to the roll rate Clp, also called the
damping derivative, is found from [79]:
Clp = Clpw+Clph+Clpv (4.26)
The wing, horizontal and vertical tail contributions are given in order as [79]:
Clpw = ( Clp=k)CL=0(k= )f(CL w)CL=(CL w)CL=0g
f(Clp) =(Clp) =0g+ ( Clp)drag
Clph = 0:5(Clp)h(Sh=S)(bh=b)2
Clpv = 2(zv=b)2Y v (4.27)
The derivative of the yawing moment with respect to the roll rate is [83],
Np = qSb
2Cn
p
2Izu0
51
The derivative of the coe cient of the yawing moment with respect to the roll rate, Cnp, is
found from [79]:
Cnp = Cnpw+Cnpv (4.28)
The wing and vertical tail contributions are given in order as [79]:
Cnpw =f(Cnp=CL)CL=0gCL+(Cnp= t) t
+[ Cnp=f( f)( f)g]( f)( f)
Cnpv = (2=b2)(lvcos +zvsin )(zvcos lvsin zv)Cy v (4.29)
The derivative of the side force with respect to the roll rate is [83],
Yp = qSbCnp2mu
0
The derivative of the coe cient of the side force with respect to the roll rate Cyp can be
found from [79]:
Cyp = 2(Cy v)f(zvcos lvsin )=bg (4.30)
52
4.2.5 Pitch Rate Derivatives
The pitching moment with respect to pitch rate derivative, also called pitch damping
derivative is calculated from [79]:
Cmq = Cmqw +Cmqh (4.31)
The wing contribution is given by [79]:
Cmqw = Cmqw=atM=0f A
3tan2 c=4
AB+6cos c=4 +
3
Bg=f
A3tan2 c=4
A+6cos c=4 +3g
where Cmqw=atM=0 = kwCl wcos c=4[Af2(xw=c)
2+0:5(xw=c)g
(A+2cos c=4) +
A3tan2 c=4
24(A+6cos c=4)+
1
8]
The horizontal tail contribution is given by [79]:
Cmqh = 2(CL h) hVh(xach xcg) (4.32)
4.2.6 Yaw Rate Derivatives
The derivative of the rolling moment with respect to the yaw rate is [83],
Lr = qSb
2Cl
r
2Ixu0
53
The coe cient of the rolling moment with respect to the yaw rate Clr is determined
from [79]:
Clr = Clrw+Clrv (4.33)
The wing and vertical tail contributions are given in order as [79]:
Clrw = (CLw)(Clr=CL)CL=0+( Clr= ) +( Clr= t) t
+[ Clr=f( f)( f)g]( f)( f)
Clrv = (2=b2)(lvcos +zvsin )(zvcos lvsin )Cy v (4.34)
The yawing moment with respect to the yaw rate is [83],
Nr = qSb
2Cn
r
2Izu0
The coe cient of the yawing moment with respect to the yaw rate Cnr, also called the
yaw damping is found from [79]:
Cnr = Cnrw+Cnrv (4.35)
The wing and vertical tail contributions are given in order as [79]:
Cnrw = (Cnr=CL2)(CLw)2+(Cnr=CD0)CD0w
54
Cnrv = (2=b2)(lvcos +zvsin )2Cy v (4.36)
The side force with respect to the yaw rate is [83],
Yr = qSbCyr2mu
0
The coe cient of the side force with respect to the yaw rate Cyr, is found from [79]:
Cyr = 2(Cy v)f(lvcos +zvsin )=bg (4.37)
4.3 Control Derivatives
The techniques to calculate the following control derivatives are addressed in this sec-
tion:
1. Aileron control derivatives.
2. Elevator control derivatives.
3. Rudder control derivatives.
As it was stated before, only preliminary equations are addressed here. Complete analysis
can be found from [79-83].
55
4.3.1 Aileron Control Derivatives
The rolling moment with respect to the aileron de ection is given as [83]
L a = qSbCl aI
x
The coe cient of the rolling moment with respect to the aileron de ection, Cl a, also
called the roll control power, is calculated from [83]:
Cl a = 2CL w Sb
Z y2
y1
cy dy (4.38)
The yawing moment with respect to the aileron de ection is given as [83]
N a = qSbCn aI
z
The coe cient of the yawing moment with respect to the aileron de ection, Cn a, also
called the adverse aileron yaw, is predicted from [83]:
Cn a = KaCLwCl a (4.39)
The side force with respect to the aileron de ection is given as [83]
Y a = qSCy am
56
The coe cient of the side force with respect to the aileron de ection, Cy a is negligible
for most of the conventional aileron arrangements [80]:
Cy a = 0 (4.40)
4.3.2 Elevator Control Derivatives
The pitching moment with respect to elevator derivative, Cm e, also called the elevator
control power, can be found from [79]:
Cm e = ( e)Cmih
( e) = kbfcl =(cl )theoryg(cl )theory(k0=cl h)[( )CL=( )cl]
Cmih = (CL h) hVh (4.41)
4.3.3 Rudder Control Derivatives
The rolling moment with respect to the rudder de ection is given as [83]
L r = qSbCl rI
x
The coe cient of the rolling moment with respect to the aileron de ection, Cl r, also
called the roll control power, is calculated from [83]:
Cl r =f(zvcos lvsin )=bgCy v (4.42)
57
The yawing moment with respect to the rudder de ection is given as [83]
N r = qSbCn rI
z
The coe cient of the yawing moment with respect to the aileron de ection, Cn r, also
called the rudder control power, can be found from [79]:
Cn r = Cy vf(lvcos +zvsin )=bg (4.43)
The side force with respect to the rudder de ection is given as [83]
Y r = qSCy rm
The coe cient of the side force with respect to the rudder de ection, Cy r is given as
[79]:
Cy r = (CL v)(k0Kb)f( )CL=( )Clg( )Cl(Sv=S) (4.44)
58
Chapter 5
Investigated Examples of System Identification Using Flight Test Data
Two separate experiments were performed to investigate the System Identi cation the-
ory, which was detailed in Chapters 2 and 3. The lter error technique was used in both
experiments.
Sections 5.1-5.2 provide experimental details, which were originally carried out at
DFVLR, Institute of Flight Mechanics, Braunschweig, Federal Republic of Germany [68].
The rst experiment of this thesis utilized the same ight parameters and ight test data,
which were used in the original experiments. The results of the experiments and compar-
ison with the earlier experiments [1, 68] are also presented. The second experiment was
performed on a radio controlled, glow-powered aircraft trainer, the Hanger 9 Extra Easy.
The ight test data were acquired at Auburn University, Auburn, AL, USA.
5.1 Estimation of Lateral Motion Derivatives using the Simulated Flight Data
This experiment was performed to estimate the lateral motion derivatives of an aircraft
using a linear model and simulated aircraft response. The aircraft response was generated
incorporating moderate to high level of turbulence [1]. Nominal values of the lateral motion
aerodynamic derivatives used correspond to those obtained by parameter estimation from
ight data recorded during the tests in a steady atmosphere with the research aircraft de
Havilland DHC-2, shown in Fig. 5.1 [68]. Classical equations of aircraft motion including
59
Figure 5.1: The DLR de Havilland DHC II Research Aircraft [84]
additional state and measurement noise were used to generate aircraft response time histo-
ries. The rudder and aileron excitations are applied in the experiment to provide realistic
control inputs [68].
An aircraft model for lateral directional motion is given below, which was used to
estimate the derivatives [68, 85, 86].
State equations:
_p = Lpp+Lrr +L a a +L r r +Lvv +bx_p
_r = Npp+Nrr +N a a +N r r +Nvv +bx_r (5.1)
60
Observation equations:
_pm = Lpp+Lrr +L a a +L r r +Lvv +by_p
_rm = Npp+Nrr +N a a +N r r +Nvv +by_r
_aym = Ypp+Yrr +Y a a +Y r r +Yvv +byay
pm = p+byp
rm = r +byr (5.2)
The unknown parameter vector for the considered model has dimensional derivatives and
bias terms. It is given as [68]:
T = [Lp Lr L a L r Lv Np Nr N a N r Nv Yp Yr
Y a Y r Yv bx_p bx_r by_p by_r byay byp byr fpp frr] (5.3)
Each element of the parameter vector is de ned in the above equations except fpp and
frr. They form diagonal process noise distribution matrix F. The initial values were kept
the same as they were in the original experiment [68].
5.1.1 Comparison of the Results with the Original Experiment
The estimated results for the lateral motion derivatives, applying the system identi -
cation techniques discussed in Chapters 2 and 3, are summarized in Table 5.1. Initial and
nominal values of the derivatives [1] are also given in the table. The results obtained from
61
Table 5.1: Estimated dimensional derivatives of lateral motion of aircraft
Parameter Initial Value Nominal Value Estimated Value Estimated Value
from this from the original
experiment experiment [68]
Lp -6.700 -5.820 -5.818 (0.33) -5.718 (6.6)
Lr 1.830 1.782 1.781 (0.42) 1.720 (9.1)
L a -18.300 -16.434 -16.431 (0.53) -14.925 (11.1)
L r 0.430 0.434 0.444 (4.80) 0.200 (216)
Lv -0.114 -0.097 -0.097 (0.59) -0.0883 (12.8)
Np -0.906 -0.665 -0.667 (1.23) -0.621 (9.4)
Nr -0.665 -0.712 -0.711 (0.45) -0.722 (3.4)
N a -0.660 -0.428 -0.427 (8.61) -0.427 (58.9)
N r -2.820 -2.824 -2.825 (0.32) -2.828 (2.4)
Nv 0.0069 0.0084 0.0084 (2.88) 0.0092 (19.0)
Yp -0.640 -0.278 -0.272 (8.79) -0.297 (27.7)
Yr 1.300 1.410 1.408 (0.66) 1.415 (2.5)
Y a -1.400 -0.447 -0.449 (23.91) -0.514 (72.9)
Y r 2.790 2.657 2.659 (1.02) 2.688 (3.7)
Yv -0.193 -0.180 -0.180 (0.39) -0.180 (1.5)
Iterations 6 10
Cost Function 6.406 10 32 2.146 10 12
The values in parenthesis indicate percent standard deviation
this experiment are compared with the results of the original experiment [68]. The com-
parison of the estimated values shows that the results from the present work are as close
to the nominal values as the results from the original experiment. The lter error method
converged within 6 iterations. Neither singularity nor convergence problems were found to
occur.
After utilizing the statistical estimation techniques, the rst and the most question
is asked about the statistical accuracy of the estimates which are obtained in these ex-
periments. The maximum likelihood estimator is asymptotically more e cient [1] in the
sense of achieving smaller standard deviations which can be computed by i=P1 1. The
62
asymptotic e ciency is an important property with practical importance. It implies that
the maximum likelihood estimator makes e cient use of the available data and that the
smaller standard deviation indicates theoretically maximum achievable accuracy of the es-
timates [1]. Tables 5.1-5.5 indicate that the standard deviation of all estimated derivatives
are smaller, which means the estimates are very accurate.
5.1.2 Responses and Results
The comparison of the estimated model response with the measured data and the time
history plots of the control inputs are shown in Fig. 5.2. As it was stated before, the control
inputs and ight data used in this experiment are same as the original experiment [68]. The
estimated data are shown in green while the measured data are shown as blue. Fig. 5.3
shows the convergence plots of the estimated parameters, which suggest all parameters are
converging to a speci c constant value within 6 iterations.
Reduction of the cost function appears in Fig. 5.4, which shows how rapidly the cost
function reduces. This is due to the modi ed Newton-Raphson technique, which nds zero
for the gradient of the cost function at the local minimum point. Also, the theoretically
predicted parameters are very close to the estimated parameters using the algorithm. This
is the prime reason the cost function is very small. Recalling the formulations for process
noise discussed in Section 3.1.2, the combined formulation was used in this experiment.
From Table 5.1 and Figs. 5.2-5.4 the following points are observed.
1. The measured response and the estimated response are in agreement with each other.
63
Figure 5.2: Comparison of Estimated(green) and Measured(blue) Responses for the Lateral
Directional Motion of Aircraft
64
Figure 5.3: Convergence Plots of Dimensional Derivatives for Lateral Directional Motion of
Aircraft
65
Figure 5.4: Cost Function Reduction for Lateral Directional Motion of Aircraft
2. All derivatives converged smoothly within only 4 iterations in comparison with the
original experiment [68], where they converged in 10 iterations. Also, the estimated
derivatives of this experiment are very close to the nominal values in comparison with
the original one [68].
3. The cost function was reduced to a very low value, which was approximately zero.
The observations showed the combined formulation works well with process and measure-
ment noise.
66
Figure 5.5: The DLR HFB 320 Research Aircraft [87]
5.2 Estimation of Longitudinal Motion Derivatives of HFB-320 Aircraft
This experiment was performed to estimate the non dimensional derivatives of another
aircraft. A non-linear model of a xed wing research aircraft, the HFB-320, shown in Fig.
5.5, was used and the drag, lift and pitching moment coe cients of this research aircraft
were estimated [87]. The longitudinal motion of the aircraft was excited through ight tests
[88]. A multi-step elevator input was used in the ight test, which resulted in a short period
motion, as well as a pulse input, that lead to a phugoid motion [89]. The inputs are shown
in Fig. 5.8.
A non-linear model of longitudinal motion of aircraft was utilized in the experiment
because the non-linear model was found to provide excellent estimation results [88-90]. The
postulated non-linear model to estimate the non-dimensional aerodynamic derivatives is
67
described below [1]. This model is in terms of non-dimensional derivatives (the drag, lift
and pitching moment coe cients), which are functions of variables in the wind axes (V, ,
q, etc.) [68].
State equations:
_V = qS
mCD +gsin ( ) +
Fe
msin ( + T)
_ = qSmV CL +q + gV cos ( ) FemV sin ( + T)
_ = q
_q = qScI
y
Cm + FeI
y
(ltxsin T +ltzcos T) (5.4)
where the lift, drag and pitching moment coe cients are given as:
CD = CD0 +CDV VV
0
+CD
CL = CL0 +CLV VV
0
+CL
Cm = Cm0 +CmV VV
0
+Cm +Cmq qc2V
0
+Cm e e (5.5)
Observation equations:
Vm = V
m =
m =
68
Table 5.2: Estimated non-dimensional derivatives of longitudinal motion of aircraft
Parameter Initial Value Estimated Value from Estimated Value from the
this experiment original experiment [68]
CD0 0.0057 0.12374 (9.37) 0.1227 (2.4)
CDV 0.0077 -0.0654 (9.88) -0.0645 (3.9)
CD 0.6742 0.31952 (11.41) 0.3201 (2.2)
CL0 -0.3183 -0.09548 (13.18) -0.0952 (20)
CLV 0.2603 0.15600 (8.78) 0.1500 (10)
CL 5.3758 4.27838 (4.92) 4.3386 (1.1)
Cm0 0.0498 0.09319 (11.47) 0.1141 (3.3)
CmV 0.0189 0.01488 (14.07) 0.0022 (152)
Cm -0.4986 -0.89514 (3.24) -0.9705 (1.1)
Cmq -25.844 -38.24428 (5.93) -34.021 (2.3)
Cm e -0.9907 -1.49040 (3.66) -1.5228 (1.3)
V0 104.67 104.36 106.024
0 0.1971 0.1085 0.1117
0 0.1317 0.1543 0.1049
q0 0.0469 -0.00468 -0.00333
Iterations 6 12
Cost Function 8.514 10 30 3.1256 10 31
Values in parenthesis indicate percent standard deviation
qm = q
_qm = qScI
y
Cm + FeI
y
(ltxsin T +ltzcos T)
axm = qSmCX + Femcos T
azm = qSmCZ Femsin T (5.6)
where the longitudinal and vertical force coe cients CX and CZ are modeled as
CX = CLsin CDcos
CZ = CLcos CDsin (5.7)
69
The above mentioned non-linear aircraft model of Eqs. 5.4-5.6, contains non-linearities due
to trigonometric terms, q and also on account of the transformation in Eq. 5.7.
The unknown parameter vector for the considered model has non-dimensional deriva-
tives, and is given below as [1]:
T = [CD0 CDV CD CL0 CLV CL Cm0 CmV Cm Cmq Cm e] (5.8)
5.2.1 Comparison of the Results with the Original Experiment
The estimated results from the longitudinal motion derivatives, applying the system
identi cation technique discussed in Chapters 2 and 3, are summarized in the Table 5.2.
Initial values of the non-dimensional derivatives [1] are also given in the table. The results
obtained from this experiment are compared with the results of the original experiment
[68]. Neither singularity nor convergence problems were found to occur in this experiment.
5.2.2 Responses and Results
The comparison of the estimated model response with the measured data and the time
history plots of the control inputs are shown in Fig. 5.6. The estimated data are shown
in green while the measured data are shown in blue. Fig. 5.7 shows the convergence plots
of the estimated parameters, which suggest all derivatives converged smoothly within only
4 iterations in comparison with the original experiment [68], where they converged in 12
70
Figure 5.6: Comparison of Estimated(green) and Measured(blue) Responses for the Longi-
tudinal Directional Motion of Aircraft
iterations. Also, the estimated derivatives of this experiment are very close to the nominal
values [68].
Reduction of the cost function appears in Fig. 5.9, which shows how rapidly the cost
function was reduced.
71
Figure 5.7: Convergence Plots of Non-Dimensional Derivatives for Longitudinal Directional
Motion of Aircraft
Figure 5.8: Inputs for Longitudinal Motion
72
Figure 5.9: Cost Function Reduction for Longitudinal Motion
5.3 System Identi cation Applied to Model Aircraft Trainer
The experiments discussed in sections 5.1-5.2 were successfully completed. As a nal
part of the exercise, the system identi cation technique was applied to model aircraft trainer.
The lter error technique for linear and nonlinear systems discussed in Chapter 3 were
applied here, as well as the modi ed Newton-Raphson algorithm to reduce the cost function.
The experiment and results are discussed in the following subsections.
5.3.1 The Investigated Aircraft Trainer - Hanger 9 Extra Easy
An unmanned aerial aircraft trainer, Hanger 9 Extra Easy which is made by Hobby
Outlet (shown in Fig. 5.10), was investigated in the third experiment. The Hanger 9 Extra
73
Figure 5.10: Aircraft Trainer - Hanger 9 Extra Easy [91]
Easy is a radio controlled glow powered aircraft. Both, lateral and longitudinal motion
derivatives were estimated for this aircraft. For lateral directional motion, a linear model
was used and the state and observation equations of Eqs. 5.1-5.2 were utilized. The unknown
parameter vector of Eq. 5.3 was estimated. For longitudinal directional motion, a non-
linear model was used and the state and observation equations of Eqs. 5.4-5.7 were utilized.
The unknown parameter vector of Eq. 5.8 was estimated. The aerodynamic, geometric
and mass details of the trainer aircraft [91] are given in the following Table 5.3. Geometric
details were provided by the manufacturer and the remaining details were measured from
the aircraft.
5.3.2 Responses and Results
The stability and control derivatives of the aircraft trainer, given in Eq. 5.3 and
Eq. 5.8 were estimated in this part of experiment. The initial value of each derivative
74
Table 5.3: Aerodynamic, geometric and mass data for the aircraft trainer
Parameter Detail Data Parameter Detail Data
Wing Area, S, ft2 4.8472 Horizontal tail Area, St, ft2 1.0938
Mach number, M 0.2 Vertical tail Area, St, ft2 0.9375
Aircraft side Area, Sside, ft2 1.4637 Aircraft planform Area, Splf, ft2 1.0664
Aircraft cross sectional Area, 0.167 Aircraft wetted Area, Swet, ft2 4.3005
Scs, ft2
Wing Span, ft 5.167 Horizontal tail Span, bt, ft
Wing mean aerodynamic 0.9375 Horizontal tail mean 0.6035
chord, c, ft aerodynamic chord, ct, ft
Wing aspect ratio, A 5.5072 Horizontal tail aspect ratio, At 0.9142
Location of the wing 1=4 root 0.2614 Aircraft center of 0.2919
chord on the fuselage, fraction arm, lt, gravity, Xcg,
of the fuselage length, lf Back from the leading
edge of the wing, percent
Wing lift curve slop 4.8 Horizontal tail lift curve 3.5
CL w, rad 1 slope, CL t, rad 1
Span e ciency factor, e 0.75 Ambient air density, , slug=ft3 0.002377
Fuselage length, lf, ft 4.2917 Fuselage width, wf, ft 0.2817
Horizontal distance between 2.33 Horizontal distance between 2.52
vertical tail a.c. and vertical tail a.c. and
aircraft wing a.c., lp, ft aircraft c.g. lv, ft
Vertical distance between 0.145 Vertical distance between 0.4292
vertical tail a.c. and vertical tail a.c. and
aircraft wing a.c., zp, ft aircraft c.g., zv, ft
Aircraft weight, W, lbs 8.64 Aircraft speed V0, ft=s 58.67
Center of gravity location, 31 Compressible Sweep 1
percent of c, Correction Factor, B
measured from leading edge
Position of a.c. on wing 0.1781 Position of horizontal tail a.c. 2.7552
m.g.c., Xac, ft on wing m.g.c., Xach, ft
Aircraft mass moment of 0.2712 Sweep angle at leading 0
inertia, Iyy, slug ft2 edge, LE, deg
Angle of attack, , deg 5 sideslip angle, , deg 0
Dynamic pressure, q, lb=ft2 4.1 Taper ratio, 1
Wing dihedral angle, , deg 3 Wing twist angle, t, deg 0
Semi chord Sweep angle, c=2 0 Quarter chord Sweep angle, c=2 0
Coe cient of viscosity 0.374e-6 Thickness ratio at mean 0.12
of air, , slug=ft3, deg geometric chord, t=c, deg
Reynolds number of wing, 0.807e6 Reynolds number of fuselage, 3.7e6
RNwing RNfuse
75
Table 5.4: Comparison between the predicted values and the estimated values of the
dimensional derivatives for the lateral motion of the aircraft trainer
Parameter Predicted Values Estimated Values
Lp -8.9091 -9.5691 (5.68)
Lr 3.8293 |||
L a 25.3741 30.6561 (5.28)
L r -3.8108 |||
Lv 0.6043 2.6734 (8.65)
Np 1.0091 1.3091 (8.12)
Nr -0.4112 |||
N a -2.0027 -2.9756 (11.88)
N r -5.7743 |||
Nv 0.1953 0.3232 (18.9)
Yp -0.1964 -0.2967 (6.11)
Yr -0.0564 |||
Y a 0.0000 |||
Y r -0.3438 |||
Yv 0.5251 0.1203 (10.09)
Values in parenthesis indicate percent standard deviation
was determined using the techniques detailed in Chapter 4. The comparison between the
predicted initial values and the estimated values, of the lateral dimensional derivatives
are summarized in the Tables 5.4-5.5. The comparison between the predicted values and
the estimated values of the longitudinal non-dimensional derivatives are summarized in
the Table 5.6. The gures showing responses use the metric units for the length related
parameters and the degree units for the angle related parameters.
For the lateral motion of the aircraft trainer, two di erent data sections from the same
ight data were analyzed. The aileron and the rudder control inputs were applied simul-
taneously in both the data sections. First consider that data section in which the aileron
control input was applied. While analyzing this particular data section, the derivatives with
respect to the yaw velocity and the derivatives with respect to the rudder de ection, were
76
Figure 5.11: Comparison of Estimated(green) and Measured(blue) Responses for the Lateral
Directional Motion of Aircraft
77
Figure 5.12: Convergence Plots of Dimensional Derivatives for Lateral Directional Motion
of Aircraft
78
Figure 5.13: Cost Function Reduction for Lateral Directional Motion of Aircraft
kept constant. Table 5.4 shows comparison between the predicted values and the estimated
values of the accounted derivatives. The predicted and estimated values are very close as
that can be concluded from the Table 5.4.
The Fig. 5.11 shows comparison between the estimated and the measured parameters.
It also shows the inputs utilized for this experiment. The estimated parameters are very
close to the measured parameters. Fig. 5.12 shows the convergence of the derivatives.
Now, consider that data section in which the rudder control input was applied. While
analyzing this data section the derivatives with respect to the roll velocity and the deriva-
tives with respect to the aileron de ection, were kept constant. Table 5.5 shows comparison
79
Figure 5.14: Comparison of Estimated(green) and Measured(blue) Responses for the Lateral
Directional Motion of Aircraft
80
Figure 5.15: Convergence Plots of Dimensional Derivatives for Lateral Directional Motion
of Aircraft
81
Table 5.5: Comparison between the predicted values and the estimated values of the
dimensional derivatives for the lateral motion of the aircraft trainer
Parameter Predicted Values Estimated Values
Lp -8.9091 |||
Lr 3.8293 4.5890 (8.46)
L a 25.3741 |||
L r -3.8108 -4.7749 (19.71)
Lv 0.6043 -0.0435 (9.43)
Np 1.0091 |||
Nr -0.4112 -0.6200 (17.31)
N a -2.0027 |||
N r -5.7743 -6.5220 (7.3)
Nv 0.1953 0.0135 (25.36)
Yp -0.1964 |||
Yr -0.0564 -0.0614 (27.9)
Y a 0.0000 |||
Y r -0.3438 -0.4526 (25.26)
Yv 0.5251 -0.1016 (3.39)
Values in parenthesis indicate percent standard deviation
Figure 5.16: Cost Function Reduction for Lateral Directional Motion of Aircraft
82
Table 5.6: Comparison between the predicted values and the estimated values of the non-
dimensional derivatives for the longitudinal motion of aircraft trainer
Parameter Predicted Values Estimated Values
CD0 0.0232 |||
CDV 0.0000 |||
CD 0.0021 0.0893 (9.25)
CL0 0.2243 -0.0202 (4.55)
CLV 0.0000 |||
CL 0.1543 0.1274 (7.38)
Cm0 0.0106 -0.0194 (19.72)
CmV 0.0000 |||
Cm -0.4286 -1.2563 (10.97)
Cmq -23.7257 -30.2730 (6.36)
Cm e 2.9213 3.3677 (5.36)
Values in parenthesis indicate percent standard deviation
between the predicted values and the estimated values of the accounted derivatives. The
predicted and estimated values are very close as that can be concluded from the Table 5.5.
The Fig. 5.14 shows comparison between the estimated and the measured parameters.
It also shows the inputs utilized for this experiment. The estimated parameters are very
close to the measured parameters. Fig. 5.15 shows the convergence of the derivatives.
To estimate the longitudinal motion of the aircraft trainer, consider that data section
where the elevator control input was given. The aircraft drag due to zero lift was found
unobservable in this experiment. Hence it was not estimated and its predicted value was
used as a constant value. Table 5.6 shows comparison between the predicted values and
the estimated values of the accounted derivatives. From Table 5.6 it can be said that the
predicted and the estimated values are very close.
The Fig. 5.17 shows comparison between the estimated and the measured parameters.
Fig. 5.19 shows the inputs utilized for this experiment.
83
Figure 5.17: Comparison of Estimated(green) and Measured(blue) Responses for the Lon-
gitudinal Directional Motion of Aircraft
84
Figure 5.18: Convergence Plots of Non-Dimensional Derivatives for Longitudinal Directional
Motion of Aircraft
85
Figure 5.19: Input for Longitudinal Motion
Figure 5.20: Cost Function Reduction for Longitudinal Motion
86
5.4 Model Validation
Model validation is the nal step in the system identi cation process and must be
performed regardless of the accuracy of the parameter estimation techniques. The prime
reason for above statement is: estimates are not facts. The identi ed model must prove
that its parameters have physically acceptable values with reasonable accuracy and the best
prediction capability. Due to these arguments, the estimated parameters are compared with
the available aircraft aerodynamics data, which are obtained from wind tunnel experiments,
computational uid dynamics or theoretical predictions. While comparing these parameters,
accuracy and precision of the estimates must be taken in to account. Model validation is
very important to gain con dence in the identi ed model.
87
Chapter 6
Summary
6.1 Concluding Notes
The system identi cation technique applied to aircraft is an e cient method for aircraft
aerodynamic modeling, based on ight test data. This thesis traced several breakthroughs
and achievements in the eld of system identi cation, especially application to aircraft.
The thesis started with the de nition, meaning and signi cance of system identi cation
and also included the necessity for it. Minimization of the maximum likelihood function
was explained in detail.
Presently, this technique involves ve fundamental steps: model postulation, exper-
imental design, data analysis, parameter estimation, and model validation. All of them
were explained in detail in this thesis. The lter error technique for data analysis was dis-
cussed in detail, for both linear and nonlinear dynamical systems. Theoretical techniques
for predicting aircraft stability and control derivatives were explained in some details.
The nal part of the thesis included successfully investigated examples of system iden-
ti cation using ight test data with results and comparison of the estimated parameters
with published data. Using the e cient results derived and presented in this thesis, and
also high quality results of many other applications published recently [92-94], a conclusion
is be made that many of the application elds which were thought as growing areas few
88
years before are now more or less, well established. The system identi cation methods are
now highly matured, e cient and sophisticated tools; not only for the research purposes,
but also to support and to satisfy the requirements of the aeronautical industry.
6.2 Expected Future Applications of System Identi cation to Aircraft
There are many interesting and necessary applications of system identi cation tech-
niques to aircraft which can be studied in the future. Correlation between computational
uid dynamics (CFD) practitioners and system identi cation researchers should increase
considerably. In addition to the use of the system identi cation technique to validate the
wind tunnel and the CFD data, it should also be utilized in coordinated approaches with
the CFD, to gain advantages of the strengths of both the techniques or if one approach can
not be used e ciently somewhere, another technique can be used to ll in the gaps.
The test data of each individual parts of any aircraft can be measured using small
and inexpensive sensors. Using the modern computation capabilities, the aircraft can be
identi ed as distributed parameter models for the aircraft aerodynamics. E ects of each
individual sections can be considered as a results of that. This might not be an easy task, as
distributed models need more test data of each individual part, improved dynamical models
and the best wind tunnel and CFD data. The system identi cation researchers should
then progress towards controlled automated process and solve the particular problems if
encountered.
89
A lot of research work can be done with the nonlinear systems and unsteady aero-
dynamic models. It includes learning the basic techniques to design related experiments.
These should also ensure the necessary and critical experiments are being performed at safe
environment and at e ective cost; providing the state estimates are accurate and validated
models are appropriate. System identi cation techniques for the nonlinear systems and un-
steady aerodynamic modeling are very complex, so the experiments need to be performed
with high accuracy and should be repeated.
For few applications, the theoretical and analytical problems will arise in the system
identi cation which can not be ignored. Every new ying vehicle should be tested using the
system identi cation techniques to validate the model. The prime practice, known as the
system identi cation will continue to be an important part in the aeronautical engineering,
which involves identi cation of the mathematical model based on the measured data from
the experiment.
90
Bibliography
[1] Jategaonkar, R. V., \Flight Vehicle System Identi cation: A Time Domain Methodol-
ogy", edited by Lu, F. K., Vol. 216, Progress in Aeronautics and Astronautics, AIAA,
August, 2006, Reston, VA.
[2] Zadeh, L. A., \From Circuit Theory to System Theory," Proceedings of the IRE;
Vol. 50, May 1962, pp. 856-865.
[3] Anon., \He Has a Question for Any Answer," Interview with K. Ili by M. McCall,
Antelope Valley Press, CA, Nov. 1, 1994.
[4] Jategaonkar, R. V., \Flight Vehicle System Identi cation - Engineering Utility,"
Journal of Aircraft; Vol. 42, No. 1, 2005, pp. 11.
[5] Gauss, K. F., \Theory of the Motion of the Heavenly Bodies Moving About the Sun
in Conic Section," (Translated by Charles, H. D., 1857), Dover Publications Inc., New
York, 1847. (Translated from Theoria Motus, 1809)
[6] Bernoulli, D., \The Most Probable Choice Between Several Discrepant Observations
and the Formation Therefrom of the Most Likely Induction," Acta: Acad: Petrop:;
1777, pp. 3-33; English Translation by Kendall, M. G., \Daniel Bernoulli on Maximum
Likelihood," Biometrika; Vol. 48, No. 1, 1961, pp. 1-18.
[7] Fisher, R. A., \On an Absolute Criterion for Fitting Frequency Curves," Messenger
of Mathematics; Vol. 41, Macmillan, London, 1912, pp. 155-160.
[8] Douglas, J., \Theorems in the Inverse Problem in the Calculus of Variations,"
Proceedings of the National Academy of Science; Vol. 16, No. 3, 1940.
[9] Feldbaum, A. A., \Dual Control Theory," Automn: Remote Control; Vol. 22, 1961,
pp. 1-12 and 109-121.
[10] Cuenod, M. and Sage, A., \Comparison of Some Methods Used for Process Identi -
cation," Automatica; Vol. 4, 1968, pp. 235-269.
[11] Balakrishnan, A. V. and Peterka, V., \Identi cation in Automatic Control Systems,"
Automatica; Vol. 5, 1969, pp. 817-829.
[12] Astrom, K. J. and Eykho , P., \System Identi cation - A Survey," Automatica; Vol.
7, No. 2, 1971, pp. 123-162.
91
[13] Eykho , P., \System Identi cation, Parameter and State Estimation," John Wiley
and Sons, London, 1974.
[14] Kolmogorov, A. N., \Interpolation und Extrapolation von Stationaren Zufalligen Fol-
gen," Bulletin of the Academy of Sciences of the USSR Series Mathematics;
Vol. 5, 1941, pp. 3-14.
[15] Wiener, N., \The Extrapolation, Interpolation and Smoothing of Stationary Time Se-
ries," OSRD 370, Report to the Services, Research Project DIC-6037, Massachusetts
Inst. of Technology, Cambridge, MA, Feb. 1942.
[16] Kalman, R. E., \A New Approach to Linear Filtering and Prediction Problems,"
Transactions of the American Society of Mechanical Engineers, Series D, Journal of
Basic Engineering; Vol. 82, March 1960, pp. 35-45.
[17] Astrom, K. J., and Bohlin, T., \Numerical Identi cation of Linear Dynamic Systems
from Normal Operating Records," Proceedings of the 2nd IFAC Symposium on
Theory of Self Adaptive Control Systems (England, UK), edited by Hammond,
P. H., Plenum, New York, 1966, pp. 96-111.
[18] Bryan, G. H., \Stability in Aviation," Macmillan, London, 1911.
[19] Hamel, P. G. and Jategaonkar, R. V., \Evolution of Flight Vehicle System Identi ca-
tion," Journal of Aircraft; Vol. 33, No. 1, 1996, pp. 9-28.
[20] Glauert, H., \Analysis of Phugoids Obtained by a Recording Airspeed Indicator,"
Aeronautical Research Council R and M 576; Jan. 1919.
[21] Norton, F. H., \The Measurement of the Damping in Roll on a JN4h in Flight,"
NACA Rept. 167, 1923.
[22] Norton, F. H., \A Study of Longitudinal Dynamic Stability in Flight," NACA Rept.
170, 1923.
[23] Kroll, N., Radespiel, R., and Rossow, C., \Accurate and E cient Flow Solvers for 3D
Applications on Structured Meshes," VKI Lecture Series 1994-05, Brussels, Belgium,
March 1994.
[24] Sloo , J. W., and Schmidt, W. (eds.), \Computational Aerodynamics Based on the
Euler Equations," AG-325, AGARD, Sept. 1994.
[25] Hamel, P. G., \Determination of Aircraft Dynamic Stability and Control Parameters
from Flight Testing," LS-114, AGARD, May 1981 (Paper 10).
[26] \Parameter Estimation Techniques and Applications in Aircraft Flight Testing,"
NASA TN-D-7647, 1974.
92
[27] \Methods for Aircraft State and Parameter Identi cation," AGARD-CP-172, May
1975.
[28] Breuhaus, W. O., \Summary of Dynamic Stability and Control Flight Research Con-
ducted Utilizing a B-25J Airplane," Cornell Aeronautical Lab., Bu alo, NY, Rept.
TB-405-F-10, May 1948.
[29] Milliken, W. F. Jr., \Progress in Dynamic Stability and Control Research," Journal
of the Aeronautical Sciences; Vol. 14, No. 9, 1947, pp. 493-519.
[30] Milliken, W. F., Jr., \Dynamic Stability and Control Research," Proceedingsof the
3rdAnglo AmericanAeronauticalConference; Brighton, England, UK, 1951, pp.
447-524.
[31] Seamans, R. C. Jr., Blasingame, B. P., and Clementson, G. C., \The Pulse Method for
the Determination of Aircraft Dynamic Performance," Journal of the Aeronautical
Science; Vol. 17, No. 1, Jan. 1950, pp. 22-38.
[32] Greenberg, H., \A Survev of Methods for Determining Stability Parameters of an
Airplane From Dynamic Flight Measurements," NACA TN-2340, April 1951.
[33] Shinbrot, M., \On the Analysis of Linear and Nonlinear Dynamical Systems From
Transient-Response Data," NACA TN-3288, Dec. 1954.
[34] Howard, J., \The Determination of Lateral Stability and Control Derivatives From
Flight Data," Canadian Aeronautics Space Journal; Vol. 13, No. 3, March 1967,
pp. 126-134.
[35] Shinbrot, M., \A Least-Squares Curve Fitting Method with Applications to the Cal-
culation of Stability Coe cients from Transient Response Data," NACA TN 2341,
April 1951.
[36] Wolowicz, C. H., and Holleman, E. C., \Stability Derivative Determination from
Flight Data," Rept. 224, AGARD, Oct. 1958.
[37] Mueller, R. K., \The Graphical Solution of Stability Problems," Journal of the
Aeronautical Sciences; Vol. 4, No. 8, 1937, pp. 324-331.
[38] Doetsch, K. H., \The Time Vector Method for Stability Investigations," Royal
Aircraft Establishment Rept:; Aero 2495, Aug. 1953.
[39] Stern eld, L., \A Vector Method Approach to the Dynamic Lateral Stability of Air-
craft," Journal of the Aeronautical Sciences; Vol. 21, No. 4, 1954, pp. 251-256.
[40] Wolowicz, C. H., \Considerations in the Determination of Stability and Control
Derivatives and Dynamic Characteristics From Flight Data," AGARD-AR-549, Pt.
1, 1966.
93
[41] Rampy, J. M., and Berry, D. T., \Determination of Stability Derivatives From Flight
Test Data by Means of High Speed Repetitive Operation Analog Matching," FTC-
TDR-64-8, Edwards, CA, May 1964.
[42] Taylor, L. W. Jr., and Ili , K. W., \A Modi ed Newton-Raphson Method for Deter-
mining Stability Derivatives From Flight Data," Second International Conference
on Computing Methods in Optimization Problems; Academic Press, New York,
1969, pp. 353-364.
[43] Larson, D. B., and Fleck, J. T., \Identi cation of Parameters by the Method of
Quasilinearization," Cornell Aeronautical Lab., Bu alo, New York, Rept. 164, May
1968.
[44] Taylor, L. W. Jr., Ili , K. W., and Powers, B. G., \A Comparison of Newton-Raphson
and Other Methods for Determining Stability Derivatives From Flight Data," AIAA
Paper, 69-315, March 1969.
[45] Ili , K. W., \Aircraft Parameter Estimation" AIAA Dryden Lecture in Research for
1987.
[46] Hamel, P. G., \Flight Vehicle System Identi cation Status and Prospects," DFVLR-
Mitt. 87-22, Nov. 1987, pp. 51-90.
[47] Grove, R. D., Bowles, R. L., and Mayhew, S. C., \A Procedure for Estimating Sta-
bility and Control Parameters From Flight Test Data by Using Maximum Likelihood
Methods Employing a Real- Time Digital System," NASA TN-D-6735, 1972.
[48] Ross, A. J., and Foster, G. W., \FORTRAN Programs for the Determination of Aero-
dynamic Derivatives From Transient Longitudinal or Lateral Responses of Aircraft,"
Royal Aircraft Establishment; TR-75090, Sept. 1975.
[49] Maine, R. E., and Ili , K. W., \A FORTRAN Program for Determining Aircraft
Stability and Control Derivatives From Flight Data," NASA TN-D-7831, 1975.
[50] Nagy, C. J., \A New Method for Test and Analysis of Dynamic Stability and Control,"
Air Force Flight Test Center, Edwards, CA, AFFTC-TD-75-4, May 1976.
[51] Mehra, R. K., \Maximum Likelihood Identi cation of Aircraft Parameters,"
Proceedings of the Eleventh Joint Automatic Control Conference of the
American Automatic Control Council; Paper 18-C, Atlanta, GA, June 1970, pp.
442-444.
[52] Ili , K. W., and Maine, R. E., \Practical Aspects of Using a Maximum Likelihood
Estimatoion method to extract stability and control derivatives from ight data,"
NASA TN D-8209, April 1976, pp. 1-34
94
[53] Ili , K. W., and Maine, R. E., \Further Observations on Maximum Likelihood Esti-
mates of Stability and Control Characteristics Obtained From Flight Data," AIAA
Paper 77-1133, Aug. 1977, pp. 100-112.
[54] Jategaonkar, R. V., and Plaetschke, E., \Maximum Likelihood Estimation of Param-
eters in Linear Systems with Process and Measurement Noise," DFVLR-FB 87-20,
June 1987.
[55] Maine, R. E., and Ili , K. W., \User?s Manual for MMLE3, a General FORTRAN Pro-
gram for Maximum Likelihood Parameter Estimation," NASA TP-1563, Nov. 1980.
[56] Taylor, J. S., Powell, J. D., and Mehra, R. K., \The Use of Smoothing and Other Ad-
vanced Techniques for VTOL Aircraft Parameter Identi cation," Final Rept., Naval
Air Systems Command Contract N00019-69-C-0534, Systems Control, Inc., Palo Alto,
CA, June 1970.
[57] Gerlach, O. H., \The Determination of Stability Derivatives and Performance Char-
acteristics From Dynamic Maneuvers," Delft Univ. of Technology, Dept. of Aerospace
Engineering, Delft, The Netherlands, Rept. VTH-163, March 1971.
[58] Jonkers, H. L., \Application of the Kalman Filter to Flight Path Reconstruction From
Flight Test Data Including Estimation of Instrumental Bias Error Corrections," Delft
Univ. of Technology, Dept. of Aerospace Engineering, Delft, The Netherlands, Rept.
VTH-162, Feb. 1976.
[59] Astrom, K. J., \Control Problems in Papermaking," Proceedings of the IBM
ScientificComputingSymposiumonControlTheoryandApplications; IBM Data
Processing Division, White Plains, New York, 1966, pp. 135-167.
[60] Kashyap, R. L., \A New Method of Recursive Estimation in Discrete Linear Systems,"
IEEE TransactionsonAutomaticControl; Vol. AC-15, No. 1, Feb. 1970, pp. 18-24
[61] Chen, Robert, T. N., and Eulrich, B. J., \Parameter and Model Identi cation of
Nonlinear Dynamical Systems Using a Suboptimal Fixed-point Smoothing Algo-
rithm," Proceedings of the Twelfth Joint Automatic Control Conference of the
American Automatic Control Council; Paper 7-E2, Aug. 1971, pp. 731-740.
[62] Yazawa, K., \Identi cation of Aircraft Stability and Control Derivatives in the Pres-
ence of Turbulence," AIAA Paper 77-1134, Aug. 1977.
[63] Balakrishnan, A. V, \Stochastic System Identi cation Techniques," Stochastic
OptimizationandControl; edited by H. F. Karreman, John Wiley and sons, London,
1968.
[64] Balakrishnan, A. V., \Modelling and Identi cation Theory: A Flight Control Appli-
cation," Theory and Applications of Variable Structure Systems; edited by R. B.
Mohler and A. Ruberti, Academic, New York, 1972.
95
[65] Ili , K. W., \Identi cation and Stochastic Control with Application to Flight Control
in Turbulence," Ph.D. Dissertation, Univ. of California, Los Angeles, CA, May 1973.
[66] Etkin, B., \Dynamics of Atmospheric Flight," John Wiley and Sons, New York, 1972.
[67] Maine, R. E., and Ili , K. W., \Formulation and Implementation of a Practical al-
gorithm for Parameter Estimation with Process and Measurement Noise," SIAM J:
Applied Mathematics; Dec. 1981, pp. 558-579.
[68] Jategaonkar, R. V., and Plaetschket, E., \Algorithms for Aircraft Parameter Estima-
tion Accounting for Process and Measurement Noise," Journalof Aircraft; Vol. 26,
No 4, April 1989, pp. 360-372.
[69] Ili , K. W., and Taylor, L. W., \Determination of Stability Derivatives from Flight
Data Using a Newton-Raphson Minimization Technique," NASA TN D-6579, March
1972.
[70] Ortega, J. M., and Rheinboldt, W.C., \Iterative Solution of Nonlinear Equations in
Several Variables," Academic Press, New York, 1970.
[71] Balakrishnan, A. V., \Communication Theory," 1968, McGraw-Hill, NY.
[72] Mehra, R. K., \Identi cation of Stochastic Linear Dynamic Systems Using Kalman
Filter Representation," AIAA Journal, Vol. 9, Jan. 1971, pp. 28-31.
[73] Gelb, A., \Applied Optimal Estimation," The MIT Press, Cambridge, MA, 1974.
[74] Grewal, M. S., and Andrews, A. P., \Kalman Filtering Theory and Practice," Prentice
Hall, Upper Saddle River, NJ, 1993.
[75] Vaughan, D. R.,\A Non-recursive Algebraic Solution for the Discrete Riccati Equa-
tion," IEEE Transactions on Automatic Control; Vol. AC-15, Oct. 1970, pp. 597-
599.
[76] Potter, J. E., \Matrix Quadratic Solutions," SIAM Journal of Applied
Mathematics; Vol. 14, May 1966, pp. 496-501.
[77] Welch, G., and Bishop, G., \An Introduction to the Kalman Filter," SIGGRAPH,
Los Angeles, CA, Aug. 2001.
[78] Maybeck, P., \Stochastic models, estimation and control," Vol. 1, Academic Press,
New York, 1979.
[79] Roskam, J., \Airplane Design Part VI: Preliminary Calculation of Aerodynamic
Thrust and Power Characteristics", Vol. 2, 1990, Roskam Aviation and Engineering
Corporation, Ottawa, KS.
96
[80] Roskam, J., \Airplane Flight Dynamics and Automatic Flight Controls Part I", Vol.
1, 1995, DARcorporation, Lawrence, KS.
[81] Roskam, J., \Airplane Flight Dynamics and Automatic Flight Controls Part II", Vol.
1, 1998, DARcorporation, Lawrence, KS.
[82] Raymer, D. P., \Aircraft Design: A Conceptual Approach", Vol. 3, 1999, AIAA,
Reston, VA.
[83] Nelson, R. C., \Flight Stability and Automatic Control", Vol. 2, 1998, McGraw-Hill,
NY.
[84] Plaetschke, E., Mulder, J. A., and Breeman, J. H., \Results of Beaver Aircraft Pa-
rameter Identi cation," DFVLR-FB 83-10, March 1983.
[85] McRuer, D., Ashkenas, I., and Graham, D., \Aircraft Dynamics and Automatic Con-
trol", Princeton Univ. Press, Princeton, NJ, 1973, pp. 292-295.
[86] He ey, R. K., Jewell, W. F., Lehman, J. M., and Van Winkle, R. A., \A Compilation
and Analysis of Helicopter Handling Qualities Data, Volume One: Data Compilation,"
NASA CR-3144, Aug. 1979.
[87] Chowdhary, G. and Jategaonkar, R. V., \Aerodynamic Parameter Estimation from
Flight Data Applying Extended and Unscented Kalman Filter," AIAAAtmospheric
Flight Mechanics Conference and Exhibit; Keystone, Colorado, Aug. 2006, pp.
6146-6163.
[88] Jategaonkar, R. V., and Balakrishna, S., \E ects of Flap Position on Longitudi-
nal Parameters of HFB-320," TM SE-8602, NAL Bangalore, India, Feb. 1986 and
DFVLR-IB 111-85/15, June 1985.
[89] Mackie, D. B., \A comparison of Parameter Estimation Results from Flight Test Data
Using Linear and Nonlinear Maximum Likelihood Methods," DFVLR-FB 84-06, Dec.
1983.
[90] Ili , K. W., \Maximum Likelihood Estimation of Lift and Drag from Dynamic Aircraft
Maneuvers," Journal of Aircraft; Vol. 14, Dec. 1977, pp. 1175-1181.
[91] Hobby Outlet, 1700 S.E. Mile Hill Dr., Suite 101, Port Orchard, WA 98366.
[92] Jategaonkar, R. V., Fischenberg, D., and Gruenhagen, W., \Aerodynamic Modeling
and System Identi cation from Flight DataRecent Applications at DLR ," Journal
of Aircraft; Vol. 41, No. 4, 2004, pp. 681-691.
[93] Wang, K. C., and Ili , K. W., \Retrospective and Recent Examples of Aircraft Param-
eter Identi cation at NASA Dryden Flight Research Center ," Journal of Aircraft;
Vol. 41, No. 4, 2004, pp. 752-764.
97
[94] Morelli, E. A., and Klein, V., \Application of System Identi cation to Aircraft at
NASA Langley Research Center," Journal of Aircraft; Vol. 42, No. 1, 2005, pp.
12-25.
[95] Kuo, B. C., \Automatic Control Systems", Prentice Hall Inc., Englewood Cli s, NJ,
1962.
[96] Hildebrand, F. B., \Introduction to Numerical Analysis", McGraw-Hill, New York,
1956, pp. 319-323.
98
Appendix A
Discretization of continuous-time state equation
A.1 Need of Discretization
The state equation utilized in the thesis is given below as [1]:
_x(t) = Ax(t) +Bu(t) +bx (A.1)
This is a continuous time equation. The Kalman lter equations are discrete in nature.
Therefore, the state equation given above must be converted from a continuous-time state
equation to a discrete-time state equation. In the discrete-time case, it is assumed that the
input vector u(t) changes only at the equally spaced sampling instants. The discrete-time
state equation will be derived, which yields the exact values at t = kT, k = 0,1,2,3...N
The continuous-time state equation is Eq. A.1. Here, kT and (k+1)T are used instead
of k and (k+1)in order to clarify the analysis. The discrete-time representation of Eq. A.1
should be as follows [95]:
x(k+1)T = (t)xkT + (t)ukT (A.2)
The matrices and depend on the sampling period t. Therefore, if the period is xed,
and are constant matrices.
99
A.2 Determine and
The general integral solution of Eq. A.1 is [95]:
x(t) = eAtx(0) +eAt
Z t
0
e A Bu( )d (A.3)
Therefore
x(k+1)T = eA(k+1)Tx(0) +eA(k+1)T
Z (k+1)T
0
e A Bu( )d (A.4)
and
xkT = eAkTx(0) +eAkT
Z kT
0
e A Bu( )d (A.5)
If Eq. A.5 is multiplied by eAT and subtracted from Eq. A.4, the following result will be
obtained [95].
x(k+1)T = eATxkT +eA(k+1)T
Z (k+1)T
kT
e A Bu( )d (A.6)
That can be simpli ed to [95]:
x(k+1)T = eATxkT +
Z T
0
eA Bu( )d (A.7)
Comparing Eq. A.2 and Eq. A.7 [95].
= eAT and =
Z T
0
eA d (A.8)
100
is termed the state transition matrix, and is the integral of the state transition matrix.
= eAT
= A 1(eAT I) (A.9)
where I is an identity matrix.
The nal discrete-time state equation is given below:
ex(tk+1) = bx(tk) + Bu(tk) + bx (A.10)
101
Appendix B
Numerical Integration using Runge-Kutta method, Order 4
A Runge-Kutta algorithm was utilized in this thesis to integrate non-linear equations
numerically. Runge-Kutta algorithms use trial steps at the midpoint of each integration
interval to cancel out lower order error terms. The fourth-order formula is sometimes known
as RK44. This method is reasonably simple, robust and also a good general candidate for
numerical solution of di erential equations when combined with an intelligent adaptive
step-size routine [96]. Runge-Kutta Method of order 4 used here assumes the function to
be used is continuous-time. The Runge-Kutta method uses following formulae [96]
tk+1 = tk + t
xj+1 = xj + t6 (k1 + 2k2 + 2k3 +k4) (B.1)
Where k1, k2, k3, k4 are calculated as [96]:
k1 = k(x;u; )
k2 = k(x+k1( t=2);u; )
k3 = k(x+k1( t=2);u; )
k4 = k(x+k3 t;u; ) (B.2)
102