Hajiyev 2015
Hajiyev 2015
Sıtkı Yenal Vural
State Estimation
and Control
for Low-cost
Unmanned
Aerial Vehicles
State Estimation and Control for Low-cost
Unmanned Aerial Vehicles
Chingiz Hajiyev • Halil Ersin Soken
Sıtkı Yenal Vural
123
Chingiz Hajiyev Halil Ersin Soken
Faculty of Aeronautics and Astronautics Institute of Space and Astronautical Science
Istanbul Technical University (ISAS)
Istanbul, Turkey Japan Aerospace Exploration
Agency (JAXA)
Sıtkı Yenal Vural Sagamihara, Kanagawa, Japan
Faculty of Aeronautics and Astronautics
Istanbul Technical University
Istanbul, Turkey
v
vi Preface
In Chap. 5, a robust Kalman filter (RKF) with filter gain correction for the case
of measurement malfunctions is introduced. By the use of a defined variable named
the measurement noise scale factor (MNSF), the faulty measurements are taken into
the consideration with a small weight and the estimations are corrected without
affecting the characteristics of the accurate measurements. The RKF algorithms with
single and multiple MNSFs (R-adaptation procedure) are proposed and applied for
the state estimation process of the UAV platform in the presence of measurement
faults. The results of these algorithms are compared for different types of sensor
faults and recommendations about their utilization are given. Stability analysis for
proposed RKFs is also included.
In Chap. 6, the process noise covariance adaptation procedure (Q-adaptation) for
UAV state estimation is introduced and the robust adaptive Kalman filter (RAKF)
algorithm with the R- and Q-adaptations against sensor/actuator failures is proposed.
Thus, the filter stands robust against the faults and, even in case of sensor/actuator
failure, continues providing accurate estimation results. The performance of the
proposed RAKF is investigated by simulations for the state estimation procedure
of UAVs. The presented RAKF ensures that the parameter estimation system of
the UAV is not influenced by sensor/actuator failures and, therefore, autonomous
missions can be performed successfully.
In Chap. 7, Global Positioning System (GPS)/Inertial Navigation System (INS)
measurements, which have high accuracy, and air data system (ADS) measurements,
which have low accuracy but high frequency, are integrated using the Kalman
filtering technique in order to obtain high-accuracy measurement data at high
frequency. It is shown that the designed system, which is based on the indirect
Kalman filtering technique, is successful in finding the wind speed data using the
known error dynamics of the system and the determined statistical values. The
fault detection and isolation (FDI) algorithms are developed by using the proposed
system and diagnostic tests are performed in order to evaluate the performance of
the system under different sensor measurement fault scenarios. For fusing the data
coming from different measurement sets, a federated Kalman filter is used. The best
data are produced by using the FDI algorithms and leaving the faulty data out before
fusing with the federated filter.
The stability of UAVs is investigated in Chap. 8 for the purpose of controller
design. The transfer functions for the longitudinal and lateral stability analysis of
UAVs are introduced and discussed in this chapter.
The design procedure for the flight control system of UAVs based on the
classical controllers is presented in Chap. 9. Classical control methods, such as
the proportional-integral-derivative (PID) controllers, are investigated. It is observed
that the PID controller responses are good at controlling the longitudinal and lateral
flight dynamics of small UAVs.
In Chap. 10, an optimal controller is designed using the linear quadratic regulator
(LQR) method. In order to better evaluate the effect of disturbances on the obtained
measurements, a Kalman filter is used. Initially, the controller is tested for a case
where disturbances are absent. Then, a Kalman filter is designed and the system
Preface vii
under disturbances is tested with the designed controller and filter. The results reveal
the effectiveness of the Kalman filter and LQR controller.
In Chap. 11, a fuzzy controller for UAV flight dynamics is investigated. Fuzzy
logic-based longitudinal and lateral controllers are designed. The fuzzy logic
controller works well, even though no optimization technique is used to optimize
the system, and it was designed without dynamic model knowledge. In fuzzy logic
controllers, a careful arrangement of membership functions may lead to better
results. The stability of the fuzzy controllers is investigated. The effectiveness of
different control methods for UAV flight control is compared in this chapter.
This monograph is useful for both researchers in academia and professional
engineers in the aerospace industry. It is also an important reference book for
graduate research students. Moreover, Chaps. 2, 3, 4, 8, 9, and 10 can also be used
by undergraduate students studying aerospace and control and systems engineering.
We would like to thank Mr. Volkan Kargin for his contribution to the simulations
given in Chap. 7 and M.Sc. student Demet Cilden who helped us with the corrections
and formatting of the book.
We hope that this book will be helpful to researchers, students, our colleagues,
and all the readers.
ix
x Contents
Index . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . 225
Chapter 1
Introduction to Unmanned Aerial Vehicles
1.1 Introduction
Everything started with a passion and a desire for flying like birds. Was looking
to the world from above only a dream for humanity? Maybe that was the only
thought of Hezarfen Ahmet Celebi when he jumped down from Galata Tower four
centuries ago, for his first flight over Istanbul. Ahmet Celebi is a legendary character
in Turkish history and we will probably never be sure of his flight’s reality, but
there is one thing that we are certain of; much has changed since that date. Flight
technology has come to a point that he could not imagine, if he ever existed. Today,
aircraft are ideal vehicles of transportation, perfect weapons of warfare, and helpers
of human beings, which makes life easier in many areas. Now, think that they are
not very expensive or hard to produce and do not even need a pilot. UAVs, namely,
unmanned aerial vehicles, were developed for making this purpose real. Although
they do not have a genius director, a human aboard, for controlling extraordinary
situations without risking human life, UAVs are preferred to manned aircraft for
their properties of being more economical, smaller, and lighter.
From the beginning, as can be understood from its name, a UAV is a vehicle
without a pilot or crew aboard, and when a battle operation, risky mission, or a
special case like redirecting a guided missile is discussed, that specialty makes them
greatly superior. For years, wars have been a potential death machine for humanity
and, in that situation, having an assistant which can reduce your army’s loss of life
becomes very important. As Holder states, the availability of using them during
dangerous missions, when rapid information is needed, manned aircraft cannot be
operated, and in other situations that they are required, gives UAVs an important
role in warfare [1]. Supporting that, possibly the best example of the usage of UAVs
in battle is the Iraq war. The eyes and assistants of the ground-based coalition
forces in civil war areas of Iraq were UAVs [2]. If the Fallujah operation, a hard,
key operation of the war, is examined, that truth can be seen.West explains the
position before the operation; since the consultation between the USA and Iraqi
governments finished, they agreed on an operation to Fallujah for taking control of
that part of Baghdad again. The American and Iraqi military forces, comprising
12,000 persons, must fight with the 4,000 rebels without a serious loss of life,
which threatened both them and about 200,000 civilians in the city. At that point,
finding the location of the rebels became important and the only unit of the army
which could do that perfectly was Marine Unmanned Aerial Squadron (VMU-1),
nicknamed “Watchdogs”, and their Pioneer UAVs [2]. At the end of the operation,
after 11 days, the city was cleared and the result of the mission was nearly 30 times
lower life loss for multinational forces. This precisely demonstrates how useful
UAVs can be in a war scenario.
Another situation that suits a UAV is risky operations. Perhaps the attacking
capability of a UAV is limited, but their excellence in discovery missions, which
are described as the most risky operations, cannot be questioned [1]. For example,
Newcome declares that a life-threatening type of job for pilots is carrying out
discovery flights in nuclear test areas, which occurred with high possibility during
the Cold War, the result of which could bring about a nuclear war immediately. For
that reason, these missions were described as being “dirty”, so giving the job to
UAVs was the best approach [3]. They really are excellent in risky missions and,
from this point of view, Holder aids Newcome as he describes another advantage of
a UAV; there is no need for a risky rescue mission with the potential for more loss
of life when a UAV is lost [1]. In 2000, after the volcanic activity of Mount Usu in
Japan, instead of manned aircraft, UAVs flew over the area by the decision of the
Ministry of Construction. The reason for using UAVs for the discovery of dangerous
situations was the active volcano, deemed too risky for human pilots [3]. Also, there
have been some strange suggestions made as a result of the suitability of UAVs in
“dirty” duties. For instance, according to Holder, making UAVs a clear target for
guided missiles and not causing the loss of a pilot’s life is one of the other ideas
which is examined for the usage of UAVs [1]. When the most valuable thing of a
pilot, his or her life, is in danger, UAVs are their lifebuoys.
In today’s world, it is a truth that making an artificial bird fly is very expensive
and, because of that, more economical UAVs may be attractive. Firstly, the education
process of a pilot is very long and, related to that, not cheap, so operating a UAV for
the same jobs can be better for economy, simply because it does not need a pilot.
A combat pilot candidate can learn everything about fighting and flying in 6 weeks,
a period nearly 17 times longer than the duration of the building process of a new
Predator UAV [3]. There is no need to waste time and money. Also, they have the
possibility of being used as part of a pilots’ fighting skill development process. That
secures both a better as well as a cheaper education. Being a moveable target for
fighter aircraft weapons is a popular use for drone aircrafts. Holder defines as an
example and, he continues, for educating their fighter pilots and improving their
combat talent, the Air Force and Navy of the US Army have used various kinds of
drone aircrafts for years.
Another direct economical advantage of UAV usage is the price of a UAV. When
the verification between manned and unmanned aircraft is done in a security area,
1.1 Introduction 3
this fact can be seen more clearly. In 2005, Farivar noted that aircraft used by
police-like organizations for surveillance and crime-fighting missions demanded
a lot of money to buy, replace, and repair [4]. Hence, more economical UAVs
that can replace manned planes became more interesting. Again, according to
Farivar, the person responsible for UAV technology development for Los Angeles
County Sheriff’s Department, Cmdr. Sid Heal, explains “the plane’s low cost makes
it highly desirable, particularly when compared with the current alternative of a
helicopter” and, he continues, “that helicopter costs $450 an hour—for every 10 h
of use I’ve bought a new one of these that I own.” Also, he claims “they shoot one
of these down and I’ll have funds for 100 of them the next day” [4]. There is one
thing that is certain; people want aircraft which are cheap on account to both their
construction and flight costs, and that is valid even for armies with large budgets.
Thirdly, being unmanned reduces the size of an aircraft and, when compared
with large planes, a small plane has many more advantages. Being unmanned does
not only mean not having a pilot; it also means not carrying many more devices
necessary for the crew. For example, life support equipment is an essential part
of manned aircraft. Holder indicates that, when manned aerial craft and UAVs
are considered, not requiring life support equipment for a pilot puts UAVs a step
ahead [1]. That equipment is big and heavy, so if the same conditions are under
discussion, using a small UAV with less chance of being spotted is more reasonable;
nobody wants his or her spy plane to be seen. Some other advantages that small size
brings to a plane are long flight times and ease of transportation and taking off.
As a result, UAVs are considered a part of many missions, such as fire fighting.
DeGaspari tells his readers that, unlike the limited vision of ground-based fire
fighting personnel, aerial vehicles have the chance to explore wide fire areas.
For meeting that requirement, the laborers of the Idaho National Engineering and
Environmental Laboratory (INEEL) are working on a fire fighter UAV [5]. Scott
Bauer, program manager for the lab’s UAV program at the INEEL, says that carrying
a small UAV with a trailer, sending it airborne with the help of a catapult, and
controlling it with a suitcase-sized remote control is plausible, where he describes a
small UAV as a vehicle with 3–3.5-m wingspan, 12-kg payload capacity, and 10 h
of flight time [5]. At the same time, there are other examples of UAVs which are
easier to transport because of their portability and, surprisingly, they can be ready to
fly in a few minutes. Security UAVs, which are used by law-keeping organizations
like the police, are one of these kinds of vehicles. According to Farivar regarding
the size of a security UAV, Dr. Chang explains that the plane is portable enough to
carry in a bag and, when it is required for an operation, mounting it and making it
fly in a short time is possible [4].
Making a heavy object fly is always harder than flying a lighter one and, in that
case, the property of an UAV being lighter can be examined. First, for purposes
like security keeping, traffic controlling, and having vision for the press, operating
a UAV over cities is safer. Farivar explains the situation by using a quotation from
Prof. Kroo, the security UAVs director of the Aircraft Aerodynamics and Design
Group at Stanford University. Prof. Ilan Kroo states that, “the miniaturization of
these aircrafts is also driven by a concern for public safety” and he adds, “if
4 1 Introduction to Unmanned Aerial Vehicles
UAVs can be defined as aircraft without a human pilot aboard. According to the
International Civil Aviation Organization (ICAO), unmanned aircraft is classified
into two groups [7, 8]:
(A) Autonomous aircraft
Autonomous aircraft are considered not to be suitable for regulation due to
legal and liability issues
(B) Remotely piloted aircraft
Remotely piloted aircraft are subject to civil regulation under the ICAO and
under the relevant national aviation authority
There are also different names given to UAVs. They are called unmanned air
systems, unpiloted air systems, drones, unpiloted aerial vehicles, etc. The flight of
UAVs can be controlled using on-board computers autonomously (following the
online calculations or precoded control algorithms) or by remote control by a pilot
on the ground. UAV systems are complex systems that comprise:
• Control station used by system operators
• The UAV itself
• Communications system
UAVs can be used for many different military and civilian applications. These
include:
1.2 UAV Types and Applications 5
Civilian uses:
• Aerial photography
• Agricultural uses
• Meteorological purposes
• Traffic monitoring
• Logistics: UAVs for cargo operations
• Police authority
• Fire services
• Information services
• Research and development: as an experiment for further studies
Military uses:
• Reconnaissance
• Surveillance
• Target and decoy: functioning as an enemy missile or aircraft for aerial or
ground gunnery
• Combat: as fighter or bomber aircraft (in general) on high-risk missions
• Electronic intelligence
UAVs can be classified into groups depending on their size and capability for
specific missions. Endurance, range, weight, and maximum altitude characteristics
determine the capability of UAVs, and they are used to perform the classification.
However, these definitions are constantly changing as the technology advances and
new types of UAVs emerge. Different types of UAVs in use today are [7, 8]:
• Orbital UAVs – UAVs that may operate up to low Earth orbit. They can reach
speeds of Mach 25C and a primary mission example is launching spacecraft into
orbit.
• Hypersonic UAVs – UAVs that may reach supersonic (Mach 1–5) or hypersonic
(Mach 5C) speeds. They are suitable for flying at suborbital altitude.
• HALE (high altitude long endurance) – They operate over 15,000 m altitude and
for a period of 24 h or more. They are usually operated by the military.
• MALE (medium altitude long endurance) – They are operated at altitudes
between 5,000 and 15,000 m. They are generally used by the military.
• Tactical UAVs – Used by the military at smaller ranges. They are simpler systems
compared to MALE and HALE.
• Close-range UAVs – Used by mobile army groups. Their ranges are up to 100 km.
They are mostly used for surveillance, traffic monitoring, crop spraying, and
reconnaissance missions. Handheld UAVs may be regarded as a subcategory.
• Mini-UAVs – UAVs of up to 20 kg approximately. They may be hand-launched.
They are usually used by mobile army groups.
• Micro-UAVs – UAVs having a wingspan no greater than 150 mm. They are
required for operation in urban environments.
• Nano air vehicles – These are proposed to be on the same size scale as insects.
They can be used for ultrashort-range surveillance.
6 1 Introduction to Unmanned Aerial Vehicles
The listed types of aircraft can be either rotary wing aircraft (which are able to
land and take off vertically) or fixed-wing aircraft.
UAVs can also be classified according to their wing loading, engine type, and
power loading.
For calculating the wing loading, the total weight of the UAV is divided by the
wing area. Three classes have been created to classify UAVs according their wing
loading values. UAVs that have a wing loading value greater than 100 kg/m2 are
classified as having high loading. Those that have a wing loading value less than
100 kg/m2 but more than 50 kg/m2 are classified as having medium wing loading.
The other UAVs with wing loading values less than 50 kg/m2 are classified as having
low wing loading.
To complete different types of missions, different engines are used in UAVs.
These include turbofans, two-stroke engines, piston, rotary, turboprop, electric, and
propeller. Some examples of UAVs that use different engines are given below:
• Global Hawk – military UAV –turbofan
• Predator – military UAV – piston
• Dragon Eye – military UAV – electric
“Unmanned aerial vehicle” is a term which came into general use in the early 1990s
and replaced the previously preferred expression of “remotely piloted vehicle”
(RPV) [3]. According to Newcome [3], a UAV is:
A powered, aerial vehicle that does not carry a human operator, uses aerodynamic forces
to provide vehicle lift, can fly autonomously or be piloted remotely, can be expendable or
recoverable, and can carry a lethal or nonlethal payload. Ballistic or semiballistic vehicles,
cruise missiles, and artillery projectiles are not considered unmanned aerial vehicles.
Although UAVs is a recent concept which has become popular following uses
especially in NATO operations during the Kosovo War and American operations
like Operation Enduring Freedom, it has a history as old as manned planes from
World War I. In fact, misunderstanding their role in military missions and not seeing
their accomplishments clearly are key reasons why military planners did not think of
them as a part of operations earnestly until the mid-1980s [9], which brought about a
reduced public interest in those vehicles up to that date. However, nowadays, UAVs
are used in many public applications, such as traffic monitoring, fire fighting, etc., as
well as their crucial function in military missions, like reconnaissance, attack, etc.
Historically, the “Kettering Bug” built by Charles Kettering for the USA Navy
in 1918 may be accepted as the first UAV. This prototype was successful enough to
make the army order additional ones, but it has never flown operationally because of
the termination of World War I. On the other hand, in the light of common opinion,
it is possible to show the “Lightning Bug”, which was substantially used by the USA
during the Vietnam War between 1964 and 1972, as the ancestor of modern UAVs.
References 7
During that period, more than 3,400 combat UAVs were flown over North Vietnam,
Laos, and China by Strategic Air Command 100th Strategic Reconnaissance Wing
of the USAF [9].
In the near future, it is expected that UAV technology will make progress
dramatically. If this development phase is examined from the point of view of
this study, then certainly, it is required to be the focus of researches about the
autonomy level of UAVs. The performance of autonomous missions, where planes
have the ability to make decisions without human intervention, is directly affected
by wrong decisions and not so precise positioning, attitude determination, etc. In
this sense, first, an accurate method for the estimation of UAV states and, then, an
efficient algorithm for controlling these states are what we need for increasing the
performance of autonomous UAV missions. Several algorithms presented in this
book aim to solve these problems.
1.4 Conclusion
As stated, the properties of being unmanned, more economical, smaller, and lighter
make UAVs superior to manned planes and, at that point, not having a genius
director no longer becomes important. The thing which is fascinating is that the
visibility of UAVs’ supremacy is increasing day by day. Already, in many areas,
UAVs have taken the place of manned aircraft. Maybe one day, even passenger
planes will be UAVs and there will no longer be a need for pilots. However, a
crucial challenge awaiting us before we can use UAVs in all areas of human life
is to increase the autonomy level of these vehicles, leading to a minimum of human
interaction. Several algorithms introduced in this book aim to solve this issue.
References
1. Holder WG, Holder B (2001) Unmanned air vehicles: an illustrated study of UAVs. Schiffer
Publishing, Atglen
2. West B (2005) Nowhere to hide. Pop Mech 182(2):54–59
3. Newcome LR (2004) Unmanned aviation: a brief history of unmanned aerial vehicles. AIAA,
Reston
4. Farivar C (2005) A flying crime fighter (some assembly required). New York Times (Article
dated: January, 13, 2005), p. G7
5. DeGaspari J (2005) Unmanned fire spies [Electronic version]. Mech Eng 127(2):20–22
6. Long MW (ed) (1992) Airborne early warning system concepts. Artech House Publishing,
Norwood
7. Austin R (2010) Unmanned aircraft systems: UAVs design, development and deployment.
Wiley, Chichester
8. Fahlstrom PG, Gleason TJ (2012) Introduction to UAV systems, 4th edn. Wiley, Hoboken
9. Cook KLB (2007) The silent force multiplier: the history and role of UAVs in warfare. IEEE
aerospace conference, Big Sky, MT, USA, March 2007
Chapter 2
Equations of Motion for an Unmanned Aerial
Vehicle
Before going into deeper analyses of the equations of motion, it will be useful to
review the axis systems that will be dealt with in further steps.
The origin of the frame is located at the center of gravity of the aircraft. The x-axis
is directed toward the nose of the plane, while the y-axis is defined as coming out
of its right wing. The z-axis completes the system by the right-hand rule as it points
down through the bottom of the aircraft.
The z-axis of the system shows the center of the Earth. The system lies on the
local horizontal plane of the Earth, as the orthogonal x- and y-axes are defined in
the directions of north and east, respectively. As a general manner, the Earth axis
system is assumed to be the inertial axis for aircraft-related studies. Since aircraft
rotation rates are relatively high in comparison to the rotation rate of the Earth, this
assumption gives acceptably low error for aircraft problems and can be utilized [1,
2] (Fig. 2.1).
Under the assumptions that the aircraft is a rigid body and its mass remains constant
for a relatively short duration of time, Newton’s second law can be written as:
!
!
!
dVc
F Dm (2.1)
dt
i
where m is the mass of the aircraft, Vc is the velocity of the center of mass of
!
the airplane, F is the force acting on the plane, and subscript i indicates that the
derivation is taken with respect to the inertial frame.
If the equation is rewritten in order to fix the axis system to the aircraft body
frame, rotation of the plane with respect to the inertial frame must be taken into
!
consideration. The derivative of Vc , referred to a rotating body frame that has an
angular velocity of !
! b , can be shown by:
!
! !
! !
C!
dV c dV c
D ! b x Vc (2.2)
dt dt b
i b
!
Hence, F , the force that acts on the aircraft, becomes:
" !
! #
!
!
F Dm
dV c !
C ! b x Vc (2.3)
dt b
b
2.1 Rigid Body Equations of Motion 11
!
where Vc , the velocity of the aircraft in the body frame, is formed by u, v, and w,
b
velocities in the x, y, and z directions, respectively, as:
!
!
!
!
Vc D u i C v j C w k (2.4)
b
and !
! b , the angular velocity vector in the body frame, consists of p, r, and q, angular
rates about the x, y, and z axes (or roll, pitch, and yaw rates), respectively, as:
!
!
!
!
!b D p i Cr j Cq k (2.5)
As a result,the three force equations for the aircraft in the body axis may be given
as:
2 3 2 3 ˇ ˇ 2 3
uP ˇ! !
!
ˇ
uP C qw rv
Fx ˇ i j k ˇ
4 Fy 5 D m 4 vP 5 C m ˇˇ p q r ˇˇ D m 4 vP C ru pw 5 (2.6a)
ˇ ˇ
Fz wP ˇ u v wˇ wP C pv qu
Fx D m .Pu C qw rv/
Fy D m .vP C ru pw/ (2.6b)
Fz D m .wP C pv qu/
On the other hand, in order to derive the three moment equations of the aircraft,
a more complicated methodology may be followed based on examining an arbitrary
small elemental mass, ım, which is at a distance from the center of mass of the
aircraft and rotates about it [2] (Fig. 2.2).
r δm
δm
z
12 2 Equations of Motion for a UAV
First, the distance between the differential mass and the center of mass of the
aircraft is given by the following vector:
r D x!
! !
!
i Cy j Czk (2.7)
ım
Here, x, y, and z are the distances in the x, y, and z axes. Now, the velocity of that
differential mass with respect to the center of mass should be written as:
!
!
d!
r
C!
!b x!
r
ım
V ım D ım (2.8)
dt
b
!
d!
r
ım
D0 (2.9)
dt
b
Therefore:
ˇ ˇ
ˇ! !
!
ˇ
!
ˇ i j k ˇ
ˇ ˇ
V ım D !
!b x!
r
ım Dˇ p q rˇ (2.10a)
ˇ ˇ
ˇ x y zˇ
!
!
!
!
) V ım D .qz ry/ i C .rx pz/ j C .py qx/ k (2.10b)
Hence, the three terms of the angular momentum can be obtained as:
2 3 ˇˇ !
!
!
ˇ
ˇ
ıHx ˇ i j k ˇ
!
ˇ ˇ
ı H ım D 4 ıHy 5 D ˇ x y z ˇ (2.12a)
ˇ ˇ
ıHz ˇ ım .qz ry/ ım .rx pz/ ım .py qx/ ˇ
ıHx D p y2 C z2 ım qxy ım rxz ım
) ıHy D q x2 C z2 ım ryz ım pxy ım (2.12b)
ıHz D r x2 C y2 ım pxz ım qyz ım
After that, the expressions should be integrated over the whole aircraft to
determine the angular momentum equations of plane itself. As long as p, r, and
q are not dependent on mass, then:
2.1 Rigid Body Equations of Motion 13
• • • •
2 2
Hx D ıHx D p y C z ım q xy ım r xz ım
• • • •
Hy D ıHy D q x2 C z2 ım r yz ım p xy ım (2.13)
• • • •
2 2
Hz D ıHz D r x C y ım p xz ım q yz ım
•
Iyz D yzım (2.15b)
•
Ixz D xzım (2.15c)
or in matrix notation:
2 3 2 32 3
Hx Ixx Ixy Ixz p
4 Hy 5 D 4 Ixy Iyy Iyz 5 4 q 5 (2.16b)
Hz Ixz Iyz Izz r
14 2 Equations of Motion for a UAV
In general, the x–z plane of the aircraft is assumed to be a plane of symmetry and,
so:
Hx D pIxx rIxz
Hy D qIyy (2.18)
Hz D rIzz pIxz
As the final stage, the time derivative of the angular momentum vector should
!
be determined. In a similar way to Eq. (2.2), the time derivative of H , referred to a
rotating body frame that has an angular velocity of !
! b , can be expressed as:
!
! !
!
!
C!
dH dH
D !b x Hb (2.19)
dt dt
i b
Since it is assumed that the mass and the mass distribution of the plane are
constant, the time moments of inertias and the products of inertias do not vary with
time:
These yields:
2 3
!
! !
! pP Ixx C qr Izz Iyy .Pr C pq/ Ixz
dH dH !
D C!
! b x H b D 4 qP Iyy pr .Izz Ixx / C p2 r2 Ixz 5
dt dt
i b rP Izz C pq Iyy Ixx C .qr C pP / Ixz
(2.23)
2.2 Orientation and Position of an Aircraft 15
As a result, the three moment equations of the motion of the aircraft in the body
axis system can be written as:
L D pP Ixx C qr Izz Iyy .Pr C pq/ Ixz
2 2
M D qP Iyy pr .I
zz Ixx /C p r Ixz
(2.24)
N D rP Izz C pq Iyy Ixx C .qr C pP / Ixz
Here, L is the moment around the x-axis, M is the moment around the y–axis, and
N is the moment around the z-axis. Also, it is good to emphasize that the first terms
pIxx ; qIyy ; rIzz , are related to the angular
of all three equation, acceleration, while the
second terms, qr Izz Iyy ; pr.Izz Ixx / ; pq Iyy Ixx , denote
gyro precision, and
2 2
the third expressions, r C pq Ixz ; p r Ixz ; qr C p Ixz , are coupling terms.
Therefore, in summary, the three force and three moment equations for an aircraft
are:
Fx D m .Pu C qw rv/
Force Eq: )
Fy D m .vP C ru pw/
Fz D m .wP C pv qu/
L D pP Ixx C qr Izz Iyy .Pr C pq/ Ixz
Moment Eq: ) M D qP Iyy pr .Izz Ixx / C p2 r2 Ixz
N D rP Izz C pq Iyy Ixx C .qr C pP / Ixz
The equations of motion of an aircraft examined thus far are derived for the body
frame of the plane. However, the orientation and position of the aircraft cannot be
defined in a relatively moving frame and, instead, a fixed frame such as the Earth
axis system must be used.
Hence, it is necessary to introduce some methodology to relate these two axis
systems: the fixed inertial frame and the aircraft body frame. To achieve this, the
orientation of the aircraft can be described by three consecutive rotations, and these
angular rotations are called Euler angles [1], a concept developed by Leonhard Euler
to define the orientation of a rigid body in 3D Euclidean space.
For such a problem, the order of rotations is important and the 3-2-1 Euler angle
representation is used. That means:
• First, rotate about the z-axis through the yaw angle
• Then, as the second step, rotate about the y-axis through the pitch angle
• And, finally, rotate about the z-axis through the roll angle of
As a consequence, it is possible to associate flight velocity components relative
to the fixed reference frame and the velocity components in the body frame:
16 2 Equations of Motion for a UAV
Fig. 2.3 Euler angles, velocity components, and angular velocities for an unmanned aerial vehicle
(UAV)
2 dx 3 2 32 3
dt cc ssc cs csc ss u
4 dy 5 D 4 cs sss C cc css sc 5 4 v 5 (2.25)
dt
dz
dt
s sc cc w
or, inversely:
2 3 2 32 3
P 1 s tan c tan p
4 P 5 D 4 0 c s 5 4 q 5 (2.26b)
P 0 s sec c sec r
As long as the six equations of motion of the plane, as derived in the previous
section, are nonlinear and must be linearized for solving behavior problems of
aircraft, small perturbation theory must be used to progress. First of all, in order
to apply the method, it is assumed that the motion of the aircraft consists of small
perturbations about the steady flight condition [1]. Then, as a base for the theory,
small perturbation terms are added to the steady states:
u D u0 C u p D p0 C p ™ D ™0 C ™
v D v0 C v q D q0 C q D 0 C (2.27)
w D w0 C w r D r0 C r D 0 C
u D u0 C u p D p ™ D ™0 C ™
v D v q D q D (2.28)
w D w r D r D
X D m .Pu/ (2.31)
If the same method is followed for all six equations of motion, linearized
equations may be obtained as follows:
2 3 2 3
X Pu
14
Y 5 D 4 vP ru0 5 (2.32)
m
Z wP qu0
18 2 Equations of Motion for a UAV
2 3 2 3
L Ixx pP C Izz rP
4 M 5 D 4 Iyy qP 5 (2.33)
N Izz rP C Ixz pP
@X @X @X @X
X D u C w C ıe C ıT (2.34)
@u @w @ıe @ıT
where ı e , ı T are control-related terms and stand for the elevator deflection and
change in thrust, respectively. Then, these two equations can be related as follows:
@X @X @X @X
m .Pu/ D u C w C ıe C ıT (2.35)
@u @w @ıe @ıT
If the same kind of association is built for all six linearized equations and, after
that, necessary simplifications are realized under the light of assumptions, the result
will be the longitudinal and lateral equations of motions for the aircraft, which is
also a consequence of small perturbation theory. To learn more about the technique,
[1] and [2] are good references.
In this final part of the chapter, first, linearized equations of motion of an airplane
are given in state space form and, then, by taking the specifications of the Zagi
unmanned aerial vehicle (UAV) [3] as a guideline, the equations are rewritten. Zagi
is a flying wing type radio-controlled aircraft that is popular amongst hobbyists.
At this point, it might be useful to indicate that we mainly take references [3–
5] as an essence to obtain scalar quantities of stability derivatives and, therefore,
stability coefficients. Hence, the values are accepted as true and used as they are in
general. Some small modifications are applied, like changing the initial conditions,
if necessary.
2 3 2 32 3
Pu Xu Xw 0 g u
6 wP 7 6 07 6 7
6 7 6 Zu Zw u0 7 6 w 7
4 Pq 5 D 4 Mu C MwP Zu Mw C MwP Zw Mq C Mw u0 0 5 4 q 5
P 0 0 1 0
2 3
Xıe XıT
6 Z ZıT 7 ıe
C4 6 ıe 7 (2.36a)
Mıe C MwP Zıe MıT C MwP ZıT 5 ıT
0 0
Fig. 2.4 The Zagi UAV used by the students of Brigham Young University [4]
20 2 Equations of Motion for a UAV
In Table 2.1, m is the mass, S is the wing area, b is the span length, c is the mean
chord length, and J is the inertia matrix of the plane. Notwithstanding, the stability
coefficients of the Zagi UAV, which are used in order to obtain stability derivatives,
may be tabulated as in Table 2.1.
In this last step, stability derivatives calculated for the Zagi UAV, in the light of
Tables 2.1 and 2.2 and related equations from [5], are substituted in the linearized
lateral and longitudinal equations of motion. However, before that process, the
2.4 Linearized Equations of Motion 21
equations are modified in such a way that observing the height of the UAV is also
possible and the sideslip angle is used instead of the starboard velocity v:
2 3 2 32 3
Pu Xu Xw 0 g 0 u
6 w 7 6 0 07 6 7
6 7 6 Z u Z w u 0 7 6 w 7
6 7 6 76 7
6 Pq 7 D 6 Mu C MwP Zu Mw C MwP Zw Mq C Mw u0 0 0 7 6 q 7
6 P7 6 76 7
4 5 4 0 0 1 0 0 5 4 5
hP 0 1 0 u0 0 h
2 3
Xıe XıT
6 ZıT 7
6 Zıe 7
6 7 ıe
C 6 Mıe C MwP Zıe MıT C MwP ZıT 7
6 7 ıT
4 0 0 5
0 0
(2.37a)
2 3 2 3
2 3 Yˇ Yp u0 Yr g cos .0 / 2 3 Yır
ˇP 6 u0 7 ˇ 6
0
6 Pp 7 6 u0 u0 u0 7 6 w 7 6 u0 77
6 7 6L 0 7 6 7 C 6 Lıa Lır 7 ıa
4 Pr 5 D 6 ˇ
Lp Lr 74 5 6 7
4 Nˇ Np Nr 0 5 q 4 Nıa Nır 5 ır
P 0 1 0 0 0 0
(2.37b)
2 3
0:7436 6:8728
6 3:7855 0 7
6 7
6 7 ıe
C 6 47:917 0 7
6 7 ıT
4 0 0 5
0 0
(2.38a)
22 2 Equations of Motion for a UAV
2 3 2 32 3
ˇP 0:1069 0:1962 1 0:9984 ˇ
6 Pp 7 6 1:2213 1:9155 1:0096 0 7 6 7
6 7 6 7 6 w 7
4 Pr 5 D 4 1:7255 0:0919 1:7198 0 5 4 q 5
P 0 1 0 0
2 3
0 0:1855
6 8:348 0 7
C6 7 ıa (2.38b)
4 4:24 2:1272 5 ır
0 0
Note that it is also possible to show the equations in a joined way, as follows:
2 3 2 3
Pu 0:3356 1:3181 0 9:80665 0 0 0 0 0
6 wP7 6 0 7
6 7 6 1:7916 3:9003 9:8215 0 0 0 0 0 7
6 Pq 7 6 0:702 3:5375 11:392 0 0 0 0 0 0 7
6 7 6 7
6 P 7 6 7
6 7 6 0 0 1 0 0 0 0 0 0 7
6 7 6 7
6 hP 7 D 6 0 1 0 9:8215 0 0 0 0 0 7
6 7 6 7
6
6 ˇP 7
7 6
6 0 0 0 0 0 0:1069 0:1962 1 0:9984 7
7
6 Pp 7 6 0 0 0 0 0 1:2213 1:9155 1:0096 0 7
6 7 6 7
4 Pr 5 4 0 0 0 0 0 1:7255 0:0919 1:7198 0 5
P 0 0 0 0 0 0 1 0 0
2 3
0:7436 6:8728 0 0
6 3:7855 0 0 0 7
6 7
6 47:917 0 7
6 0 0 7 2 ı 3
6 7 e
6 0 0 0 0 76
6 7 6 ıT 7
7
C6 0 0 0 0 74
6 7 ıa 5
6 0 0 0 0:1855 7
6 7 ır
6 0 0 8:348 0 7
6 7
4 0 0 4:24 2:1272 5
0 0 0 0
(2.39)
where A is the system matrix, B is the control distribution matrix, U is the control
vector, and X is the state space vector, one can write the following:
References 23
2 3
0:3356 1:3181 0 9:80665 0 0 0 0 0
6 1:7916 3:9003 9:8215 0 0 0 0 0 0 7
6 0:702 3:5375 11:392 0 0 0 0 0 0 7
6 7
6 0 0 1 0 0 0 0 0 0 7
6 7
AD6
6 0 1 0 9:8215 0 0 0 0 0 77
6 0 0 0 0 0 0:1069 0:1962 1 0:9984 7
6 7
6 0 0 0 0 0 1:2213 1:9155 1:0096 0 7
4 5
0 0 0 0 0 1:7255 0:0919 1:7198 0
0 0 0 0 0 0 1 0 0
(2.41a)
2 3
0:7436 6:8728 0 0
6 3:7855 0 0 0 7
6 7
6 47:917 0 0 0 7
6 7
6 0 0 0 0 7
6 7
6 7
BD6 0 0 0 0 7 (2.41b)
6 7
6 0 0 0 0:1855 7
6 7
6 0 0 8:348 0 7
6 7
4 0 0 4:24 2:1272 5
0 0 0 0
T
U D ıe ıT ıa ır (2.41c)
T
X D u w q h ˇ p r (2.41d)
References
1. Nelson RC (1998) Flight stability and automatic control, 2nd edn. WCB/McGraw-Hill, Boston
2. Yechout TR (2003) Introduction to aircraft flight mechanics. AIAA Education Series, Reston
3. Matthews JS (2006) Adaptive control of micro air vehicles. M.Sc. thesis, Department of
Electrical and Computer Engineering, Brigham Young University, Provo, UT, USA
4. Christiansen RS (2004) Design of an autopilot for small unmanned aerial vehicles. M.Sc. thesis,
Brigham Young University, Provo, UT, USA
5. Vural Y (2007) Autopilot system design for a small unmanned aerial vehicle (in Turkish). M.Sc.
thesis, Department of Aeronautical Engineering, Istanbul Technical University, Istanbul, Turkey
Chapter 3
Navigation Systems for Unmanned Aerial
Vehicles
Rate gyros measure angular velocities, in contrast to free gyros measuring attitude
angles. Rate gyros sense the vehicle’s angular rate relative to the inertial space [6].
These rate components are the craft angular rate relative to the Earth ! nb , an angular
rate as it moves about the spherical Earth ! en , and the angular rate of the Earth itself
! ie . The vector sum of these angular rates ! ib is given by:
There are three types of gyro technology used in today’s IMU systems:
• Ring laser gyro (RLG)
• Fiber optic gyro (FOG)
• MEMS
Most current RLG sensors are single DOF sensors requiring three mechaniza-
tions for an INS implementation. A single DOF RLG is shown schematically in Fig.
3.1. This figure illustrates a triangular version of the RLG. The gyro includes a laser
as a source, a closed-path cavity and mirrors at each intermediate corner in the path,
and an interferometer/photodetector. The operation of the gyro is based on optical
and electronic phenomena rather than mechanical phenomena [6, 7].
The Sagnac effect can cause a certain rate of rotation inducement with a small
difference in the time it takes light to traverse the ring in two directions. Small
separation can occur between the frequencies of the counterpropagating beams. Two
main patterns can be introduced: standing wave pattern within the ring and beat
pattern when those two beams are interfered outside the ring. The net shift of that
interference pattern follows the rotation of the unit in the plane of the ring. At very
slow rotation rates, RLGs can have the “lock-in” condition. However, they are more
accurate than mechanical gyroscopes in most cases. Also, the frequencies of the
counterpropagating laser modes get closer and become almost identical if the ring
laser is hardly rotating.
3.3 Inertial Measurement Unit 29
Laser
Interferometer
Mirror
3.3.2 Accelerometers
Accelerometers are devices that measure acceleration along the sensing axis. The
measurement of acceleration can be integrated into an onboard computer used to
derive the UAV’s velocity and position. All accelerometers use the principle of
sensing the force on a loosely suspended mass, from which the acceleration may
be calculated. In comparison to gyros, they have more developed technology. The
accelerometer converts acceleration into an electrical signal.
Dynamic acceleration can be measured using an accelerometer. Dynamic accel-
eration is due to any force except for the gravitational force applied on a rigid body
30 3 Navigation Systems for UAV
and the static acceleration (or gravitational acceleration) is due to the gravitational
force. The output of an accelerometer can be analog or digital. In the analog case,
the output voltage is directly proportional to the acceleration. On the other hand, the
output of a digital accelerometer can be directly accessed using protocols such as
serial interfaces [9]. The principle of an accelerometer is illustrated in Fig. 3.2a, b.
If the body accelerates in the direction indicated in Fig. 3.2, the mass displace-
ment with respect to the body of the device supported on a bar by springs can
be recorded. This displacement is proportional to the measured acceleration on
the sensing axis. Displacement of the mass generates an electrical signal. For this
purpose, potentiometric, self-capacitance converters can be used.
The accelerometer senses the UAV’s inertial dynamic acceleration ai , but is
insensitive to the gravitational acceleration. The gravitational field vector gm , which
is the acceleration of mass attraction to the Earth, should be taken into account
to determine the total acceleration. Therefore, the acceleration of the UAV can be
determined as:
a D ai gm (3.2)
All of the air data parameters that are relevant to flight performance are derived
by sensing the pressures, temperatures, and flow directions surrounding the vehicle.
Free-stream pressures and temperatures are required for the computation of static
air temperature, altitude, airspeed, and Mach number [3].
An air data system provides information on quantities such as pressure, altitude,
vertical speed, calibrated airspeed, true airspeed, Mach number, static air temper-
3.4 Air Data System 31
ature, and air density ratio [1]. It is, thus, one of the key avionic systems in its
own right and forms part of the essential core of avionic subsystems required of all
modern aircraft, civil or military. In this section, a short description of the air data
system is given.
Air data quantities, such as pressure, altitude, vertical speed, calibrated airspeed,
true airspeed, true Mach number, etc., are derived from three basic measurements
by sensors connected to probes which measure [1]:
• Total (or pitot) pressure
• Static pressure
• Total (or indicated) air temperature
A basic air data system scheme is presented in Fig. 3.3.
The total pressure PT is measured by means of an absolute pressure transducer
connected to a pitot tube facing the moving airstream. This measures the impact
pressure QC exerted to bring the moving airstream to rest relative to the pitot tube
plus the static pressure PS of the free airstream, i.e.:
P T D QC C P S : (3.3)
Pressure Altitude HP
Tm Vertical Speed
HP
Pitot static probe
Total Pressure Sensor/
PT Transducer Calibrated Airspeed
VC
PT
AIR DATA
COMPUTER
Mach number M
PS PS
True Airspeed VT
H
ΔPT
Compute Pressure
Error Corrections
M
Correct PT QC = PT − PS Airspeed VC
PT PT
PS PT / PS Mach Number M
Compute
Compute Mach No.
Total Pressure PT / PS
PS
Ti True Airspeed
Compute Static Air Compute True
Ti Temp Airspeed
VT
The information obtained from the air data system is as follows: true airspeed VT ,
angle of attack ˛, and sideslip angle ˇ. Using these data, the X, Y, and Z components
of the true airspeed can be obtained [10]:
q
VT D Vx2 C Vy2 C Vz2 ; (3.4)
Vz
˛ D arctan ; (3.5)
Vx
Vy
ˇ D arcsin ; (3.6)
VT
where:
Vx D VT cos ˛ cos ˇ
Vy D VT sin ˇ (3.7)
Vz D VT sin ˛ cos ˇ
These speeds are derived according to the aircraft reference system. Later
on, these coordinates can be converted to a ground-based reference system for
navigational usage.
Measurement of angle and measurement of range methods are used at the same
time in the surface radar [11]. An integrated angle/range measurement method
is usually used in radar systems, and it determines the range D of the UAV,
azimuth angle ˛, and elevation angle ˇ. When this method is used, the coordinates
of the UAV is determined as the intersection point of the sphere state surface
(D D constant), cone state surface (ˇDconstant), and vertical plane suitable for
˛Dconstant state surface. UAV coordinates are determined with a single point
(ground station) by the help of this method and does not require difficult calcula-
tions. The following formulas are used to calculate the UAV’s Cartesian coordinates:
x D D cos ˇ cos ˛
y D D cos ˇ sin ˛ (3.11)
z D D sin ˇ
An altimeter is a device that measures the flight altitude of aircraft on the Earth.
Three types of altitudes are considered during flight:
• Absolute altitude – flight altitude according to sea level
• Relative altitude – flight altitude with respect to the takeoff or landing site
• Real altitude – altitude over the flight area
Descriptions of the altitudes are given in Fig. 3.6.
The value of absolute altitude is required for setting corridors to flight routes and
flight tests of aircrafts and engines.
The relative altitude must be known for taking off and landing.
The real altitude is for arbitrary flight status.
Real altitude
Relative
altitude Absolute
altitude
Sea level
Take off
position
2H
t1 D (3.12)
c
where H is the flight altitude and c D 3 108 m/s is the propagation velocity of radio
signals.
The registration of two signals’ receiving times is used to measure the propaga-
tion time of radio signals: the signal from the transmitter and the reflected signal.
The first signal reaches the receiver antenna after time t2 D l=c and the second after
time t1 D 2H=c. At that time, the real altitude of the flight will be equal to:
.c C l/
HD ; (3.13)
2
where D t1 – t2 and l is the range between the transmitter and receiver antennas.
Expression (3.13) is used to steady radio altimeters or radio altimeters with
impulse teleportation.
When high altitudes are measured .H 1; 500 m/, the power of the transmitted
signal decreases significantly, so, for such altitudes, impulse teleportation must be
run after steady teleportation.
An altimeter’s receiver takes the reflected impulses from the aircraft and the
Earth and, after strengthening the impulses, they are passed to an electron gun (see
Fig. 3.7b). A circular movement with angular velocity ! is acted upon the beam to
find the first impulse towards the screen and the reflected second impulse. The angle
between the impulses ® (Fig. 3.7b) is calculated via the following formula:
2H!
' D ! D (3.14)
c
For increasing the measurement accuracy, the angular velocity ! of the beam
must be increased. In this case, the height of the device’s whole pointer must be
appropriate for the measurement interval.
3.6 Altitude Measurements 37
where hB is the barometric altitude, T0 is the temperature at sea level, Tgrad is the
lapse rate, p0 is the pressure at sea level, p is the pressure measured at altitude hB , R
is the universal gas constant, and g is the local gravity.
The variation of atmospheric pressure with height is given in Fig. 3.8.
1
Atmospheric
Pressure 0.8
(x105 Pa)
0.6
0.4
0.2
0
0 2 4 6 8 10
Altitude (km)
The conversion of the measured air pressure to altitude given in Eq. (3.15) is
based on a theoretical standard atmosphere and the assumption that air is an ideal
gas. More precisely, the standard atmosphere is defined as follows [15]:
• Air is an ideal gas with gas constant R D 287 J= .Kg K/
• The pressure at sea level is p0 D 29:92 in: Hg
• The air temperature at sea level is T0 D 15 ı C
• The temperature gradient (lapse rate) is Tgrad D 0:0065 ı C=m
The pressure at sea level p0 is not constant but varies from day to day and location
to location. In practice, therefore, a value of p0 which makes the altitude estimated
using (1) correct is input into the altimeter. This value of p0 is known as the altimeter
setting.
The main source of error in barometric altitude measurement arises from dif-
ferences between the true and modeled atmospheric temperature and pressure. For
standalone barometry, altitude errors can be several hundred meters. For differential
barometry, the error increases with distance from the reference station and age of
the calibration data [16].
By calculating the sum and subtracting the divided frequencies, we can write:
fd.1C3/ C fd.2C4/ D k1 W
(3.18)
fd.1C3/ fd.2C4/ D 0
are valid.
40 3 Navigation Systems for UAV
The second equation is the difference signal. This signal is used in schemes which
determine W and ˇ automatically.
The standard deviations of Doppler measuring devices is 0.2 % in respect of the
velocity and 0.1 % in respect of the drift angle.
Typically, a Doppler sensor operating over land is able to provide measurements
of velocity to an accuracy of about 0.25 % or better. Performance is degraded during
flights over water owing to poor reflectivity, scattering of the reflected signal, giving
rise to a bias in the measurement, wave motion, tidal motion, and water currents.
However, this navigation aid offers good long-term stability and a chance to bind
the position and velocity estimates provided by an inertial navigation system [16].
The variations in the Earth’s geomagnetic field can be internal, which are related to
the dynamo within the Earth’s core and crust. Inner dynamo variations, which are
long-term variations within the core, are called secular variations. On the other hand,
solar wind–geomagnetic dipole interaction variations are relatively short-duration
variations on the order of months, days, and hours. The Earth has a magnetic field
similar to that of a bar magnet, with poles located close to its axis of rotation. This
means that the direction of the horizontal component of the Earth’s magnetic field
lies close to true north, and the magnetic north, as determined by a magnetic field
sensor, or compass, can be used as a working reference [16]. The angle between
true north and the magnetic north is not constant. It varies with the position of the
observation on the Earth and slowly with time, although it is possible to compensate
for this. The direction of the Earth’s magnetic field at any point on its surface is
defined in terms of its orientation with respect to true north, known as the angle
of “magnetic declination”, and its angle with respect to the horizontal is called
the angle of “dip”. These angles are illustrated in Fig. 3.5. Magnetometers are
widely used in aerospace as attitude determination devices. These sensors provide
the direction of the Earth’s magnetic field and its magnitude; therefore, they are
useful for determining the attitude of flying vehicles. Moreover, their low mass and
low power consumption make them attractive for UAVs. For these reasons, most
UAVs have magnetometers as part of their basic sensor package.
In the last decade, the popularity of UAVs, due to their low cost and weight, has
increased significantly. That has brought about a search for lighter but more accu-
rate sensors. Under these circumstances, the three-axis magnetometer (TAM) has
become very attractive because of its advantages, such as providing continuously
available two-axis attitude measurements, relative low cost, and almost insignificant
power demand (Fig. 3.10).
In the absence of local magnetic disturbances, magnetometers sense the compo-
nents of the Earth’s magnetic field (H0 ) acting along its sensitive axis. A TAM may
be mounted in a UAV to sense the components of the Earth’s magnetic field along its
principal axes, (Hx , Hy , Hz ). Expressing a UAV’s attitude, with respect to the local
3.9 Satellite Radio Navigation 41
coordinate frame, can be done as a direction cosine matrix, Cbn . The relationship
between the magnetic measurements and the attitude of the UAV may be written as:
2 3 2 3
Hx H0 cos ı cos
4 Hy 5 D Cnb 4 H0 cos ı sin 5 ; (3.19)
Hz H0 sin ı
where ı is the angle of dip and ” is the angle of declination. If the angle of dip is
known, as well as the variation between true and magnetic norths, estimates of a
UAV’s attitude can be made from the magnetometer measurements.
The aircraft in which the magnetic sensor is mounted will almost always have a
magnetic field that cannot be distinguished from that of the Earth’s magnetic field
and, consequently, it is necessary to compensate for this effect as well. Therefore,
in order to estimate the components of the geomagnetic field vector for determining
the attitude of a UAV accurately, the magnetometer should be calibrated precisely.
total of 24 satellites. Four satellites are necessary to calculate the position of the user.
There are also other similar systems used or planned to be used, such as GLONASS
(Russia) and Galileo (European Union).
In a GPS system, the GPS receiver chooses the best available satellites to obtain
the measurements. The system is available as two services, a standard positioning
system (SPS) for civilians and a precise positioning service (PPS) for military users.
The order of accuracy of SPS is 10 m and that of PPS is around 3 m.
The GPS consists of three parts: the space segment, the control segment, and the
user segment [3]. The structure of the GPS is shown in Fig. 3.11.
The U.S. Air Force develops, maintains, and operates the space and control
segments. GPS satellites broadcast signals from space, which each GPS receiver
uses to calculate its three-dimensional location (latitude, longitude, and altitude)
plus the current time.
The space segment is composed of 24–32 satellites in medium Earth orbit and
also includes the boosters required to launch them into orbit. The control segment
is composed of a master control station, ground antennas, and monitor stations. The
user segment is, in our case, composed of UAVs.
A GPS receiver calculates its position by precisely timing the signals sent by the
GPS satellites. Each GPS satellite continually transmits messages, which include:
• The time the message was transmitted
• Precise orbital information (the ephemeris) of satellites
• The general system health
The receiver utilizes the messages it receives to determine the transit time of
each message and computes the distances to each satellite. These distances, along
with the satellites’ locations, are used to compute the position of the receiver. This
position is determined as the geographical coordinate system latitude, longitude,
and elevation. Many GPS units give derived information, such as the direction and
speed of the user, calculated from position changes.
A very small clock error, which is the discrete time scale of the user and the
satellite system multiplied by the very large speed of light—the speed at which satel-
lite signals propagate—results in a large positional error (clock bias). Therefore,
receivers use four or more satellites to solve for the four-parameter determination:
clock bias and receiver’s location as latitude, longitude, and elevation.
Today, satellite navigation systems are used for solving the navigation problems
of most of UAVs. In this technique, the position of the object is generally detected
with the help of distance measuring methods (distances from object to n 4
satellites are measured) (see Fig. 3.12).
Then, the coordinates of the object will be on the intersection of the sphere
surfaces with radiuses equal to the measured distances; see [3, 11, 17]. The
equation system that consists of n 4 nonlinear equations, which are obtained by
this approach, can be solved using a computer with the aid of iteration methods.
However, the distance equations used for formerly stated purposes have a term
other than the unknown coordinates of the user, which results from the discrete
time scale of the user and the satellite system [3, 18]. In this case, this term and the
user coordinates are determined as the intersection point of the n different sphere
position surfaces (R D constant). In a simple case, the user coordinates and clock
bias can be obtained by solving the four nonlinear algebraic equations below:
where xi , yi , and zi are the Descartes coordinates of the ith satellite, in the geocentric
axis system, 4 is the minimum number of satellites, x, y, and z are the Descartes
coordinates of the object (users) in the geocentric coordinate system, and ıt is the
distance measuring errors, resulting from the time differences between the user and
the satellite (clock bias).
UAVs are now able to work in different environments easily using a GPS receiver.
This allows the use of MALE UAVs (a medium-altitude, long-endurance UAV) and
HALE UAVs (a high-altitude, long-endurance UAV) in every part of the world.
However, in some places, by jamming or with local conditions, GPS receivers may
receive weak signals. Thus, depending solely on GPS may cause problems in the
navigation of UAVs. To overcome this problem, dead reckoning systems can be
used. A filter, such as a Kalman filter, may be used to integrate the measurements of
GPS with dead reckoning devices to obtain better measurements.
Visual sensing can also provide data for position estimation. These kinds of systems
have the advantage of using visual sensors that actually sense the real physical data
to comprehend the results in an easier way. The data from the visual sensors can also
be combined with the data from INS, GPS, and magnetometers in order to supply a
better understanding to the operator. Improved digital imaging techniques make the
use of visual sensors in small UAVs navigation more practical and simpler. Also,
using methods for emulating flying insects is a new approach of developed systems
in visual sensing. One method is the use of a newly developed artificial compound
insect eye, which is a system that replicates the eye structure of an insect. This
system uses many small-sized lenses with hemispherical shapes. The camera chips
in this system are made of hard and soft composite materials, which allow stretching
and bending, therefore acting like an insect’s eye, as desired [19–21].
SLAM algorithms are used to develop systems that are terrain aided. These are used
for map building. Using generated maps, these systems help other sensors. A few
of these systems are used on robotic systems that depend on visual sensors. Visual
motion estimation techniques that use stereovision and visual features tracking have
3.12 Measurement Fault Classification and Fault Modeling 45
been proposed in the context of ground rovers [22]. Stereovision and monocular
systems are also used in different studies [23]. These types of systems allow UAVs
to operate in all environments, including those in which they face GPS issues.
y D Km x C B; (3.21)
yi .t/ D xi .t/ C
i .t/vi .t/ (3.22)
a b
Actual Measured Actual Measured
y (t) y (t)
time time
c d
Actual Measured Actual Measured
y (t) y (t)
tF
tF
time time
e f
Measurement Noise Increment
Actual Measured
Measurement Noise
y (t)
v (t)
time
x(t)
Fig. 3.13 Several types of sensor faults: (a) sensor bias, (b) sensor drift, (c) loss of accuracy, (d)
frozen sensor (occurring after tF ), (e) calibration error, and (f) noise increment
3.12 Measurement Fault Classification and Fault Modeling 47
For example, all types of accelerometers and gyros exhibit biases, scale factor
and misalignment errors, and random noise to a certain extent [27]. The bias is a
constant error which manifests in the measurement results. Mostly, accelerometers
and gyros have bias errors. It is independent of the underlying specific force and
angular rate.
The scale factor error is the departure of the input–output gradient of the
instrument from unity following unit conversion by the IMU. The accelerometer
output error due to the scale factor error is proportional to the true specific force
along the sensitive axes, while the gyro output error due to the scale factor error is
proportional to the true angular rate about the sensitive axis [27].
Misalignment errors in all types of IMUs arise from the misalignment of the
sensitive axes of the inertial sensors with respect to the body frame orthogonal axes
because of manufacturing limitations.
To give an example regarding this subject, the reasons for INS sensor faults may
be listed as follows:
• Gyroscope drift faults caused by temperature changes, acceleration, magnetic
field, and vibration
• Gyroscope measurement faults with nonlinearity and asymmetry
• Accelerometer bias because of temperature changes and vibrations
• Accelerometer measurement faults with nonlinearity and asymmetry
• The input axes of gyroscopes and accelerometers are not perpendicular to each
other
• Initial state faults (must be considered that position and velocity errors are
proportional to time) etc.
Sensor faults may cause substantial performance degradation of a decision-
making system. For example, in a UAV flight control system, sensors are used
to measure directly the UAV states. Thus, the presence of faults in sensors may
deteriorate state estimates and, consequently, cause inaccurate flight control.
A fault monitoring system which is used to detect faults and diagnose their
location and significance within a system is called a “fault diagnosis system”. Such
a system normally consists of the following tasks [28]:
• Fault detection: to make a binary decision – either that something has gone wrong
or that everything is fine
• Fault isolation: to determine the location of the fault, e.g., which sensor or
actuator has become faulty
• Fault identification: to estimate the size and type or nature of the fault
A fault-tolerant flight control system is capable of controlling UAV dynamics
with satisfactory performance, even if there is one or more faults.
The fault-tolerant state estimation system of UAVs in this book is capable
of detecting and diagnosing faults in the UAV control system and adequately
compensate for failures.
48 3 Navigation Systems for UAV
References
4.1 Introduction
The equations of motion for an unmanned aerial vehicle (UAV) are given in
Chap. 2. In this chapter, the equations of motion are further examined and their
linearization is discussed in detail for the stability analysis. In this procedure, the
first step is investigating the trim conditions.
In this chapter, we introduce and briefly discuss the mathematical model and
the relations of the linear discrete-time Kalman filter and its application to UAV
dynamics.
Let us take a linear discrete dynamic system. The system’s dynamic state
equation defines the dynamics of the system, whereas the measurement equation
defines the generation mechanism of the measurement. The equations for a linear
system are written as follows:
State equation:
Measurement equation:
zk D Hk xk C vk (4.2)
1; k D j;
ıkj D
0; k ¤ j:
Depending on the application, one might want to obtain an estimate of the state
at a certain time. If the state is estimated for some time in the future, the process is
called prediction. If the estimate is made using all measurements up to and including
the current moment, it is called filtering. If an estimate is made for some time in the
past using measurements up to the current moment, the process is called smoothing.
In this chapter, we limit ourselves to prediction and filtering.
Generally, the filtering process of Kalman filters can be examined in two distinct
phases and if the equations are presented in the following way [3]:
Time Update (Prediction) The phase where the estimations of the preceding step
is used for producing the estimations of the present step. It can be thought of as
preparation to the estimation.
Predicted state:
xQ k=k1 D Fkb
xk1=k1 C Bk uk1 (4.3)
Here, bxk1=k1 is the estimated state from the previous step, exk=k1 is the predicted
state of the current step, Pk=k1 is the predicted covariance matrix of the current step,
Qk1 is the covariance of process noise w in discrete form, Fk and Bk are system
dynamics and control distribution matrices, respectively, in discretized form, and
uk1 is the control vector . Subscript k=k 1 denotes that the computation is done
in the current step k by using measurements of step k 1: Similarly, the subscript
k 1=k 1 means that the variable is computed at step k 1 by taking measurements
of step k 1 into consideration.
At the same time, Qk1 should be introduced at that point. First, the continuous
process noise covariance Q is:
Q D E wwT (4.5)
where w is the process noise derived from a zero-mean Gaussian distribution and
E(.) is the operator of expected value.1
1
In mathematics, the expected value of a random variable is the sum of the probabilities of each
possible outcome multiplied by the numerical value of the outcome.
54 4 Estimation of the UAV Dynamics
Zt
Qk1 D F ./ QF T ./ d (4.6)
0
Measurement Update The phase that the estimation of the present step is realized
by using prediction phase outputs and, also, measurements information operation,
where the predictions are improved.
Innovation or measurement residual:
eQ k D zk Hk xQ k=k1 (4.7)
xk=k D xQ k=k1 C Kk eQ k
b (4.10)
Here, eek is the innovation (or measurement residual), zk is the measurement vector,
Hk is the measurement matrix, Kk is the optimal Kalman gain, xO k=k is the estimated
state, Pk/k is the covariance matrix of the present step, and Rk is the covariance of
measurement noise v in discrete form.
If discrete measurement noise covariance Rk is introduced, then similarly to
Qk1 :
Rk D E vk vkT (4.12)
The OKF uses a filter gain where the expected value of the square of the magnitude
of error in posterior state estimation is minimized. In other words, the filter runs
under some certain optimization law defined by the minimization rule of the
indicated vector and, so, it has the optimal gain. If this optimal gain is modified
in order to adapt the filter to the changing conditions, this means that the filter is no
longer optimal and can be called an adaptive filter.
To make the meaning of optimal gain for an OKF more understandable, the gain
derivation process can be described as the following [4].
First, it is known that the matrix Pk/k is the covariance of estimation and, so:
Pk=k D cov xk b
xk=k (4.13)
where cov() is the covariance operator, xk is the real value of the states in step k,
and b
xk=k is the estimation of the states in step k.
Then, by using the definition of state estimation xO k=k (Eq. 4.10):
Pk=k D cov xk xQ k=k1 C Kk eQ k (4.14)
and zk as:
zk D Hk xk C vk (4.16)
Pk=k D cov xk xQ k=k1 C Kk Hk xk C vk Hk xQ k=k1 (4.17)
56 4 Estimation of the UAV Dynamics
By collecting the necessary terms, the equation can also be written as:
Pk=k D cov .I Kk Hk / xk xQ k=k1 Kk vk (4.18)
This expression of updated covariance estimation is valid for every Kk and can
be used for obtaining the optimal Kalman gain.
The Kalman filter is a minimum mean square error estimator and, as previously
mentioned, to run the filter optimally, the minimum expected value of the square of
the magnitude of state estimation
hˇ is sought [4]. The estimation error is xk b
xk=k and
ˇ2 i
the minimized quantity is E ˇxk b xk=k ˇ . However, since Pk=k D cov xk b xk=k ,
that is the same as minimizing the trace of covariance matrix Pk/k . It has been already
found as:
Initial Conditions
P(0)
P( k − 1/ k− 1)
calculation of P( k / k− 1)
Q( k)
P ( k / k− 1) Aprior data
calculation of P( k / k)
R ( k)
Initial Conditions
P (k / k)
z( k) x (0)
x ( k / k)
x (k )
System
Model
Measurement
Model ∑ K ( k) +
z( k / k− 1)
xˆ ( k − 1/ k− 1)
w(k) v ( k)
H ( k) φ (k, k − 1)
xˆ( k / k− 1)
As can be seen, this equation is the same as Eq. (4.9) and the gain derived so far
is called the optimal Kalman gain. Besides, if this expression of Kk is substituted
into Eq. (4.22), then the equation:
Pk=k D .I Kk Hk / Pk=k1
The structural diagram of the Kalman filter [5] is shown in Fig. 4.1.
The initial values x0 and P0 , correlation matrix of the system noise Q, and
correlation matrix of the measurement noise R must be known previously in order
for the filter to work.
According to formula (4.10), the estimation value is equal to the sum of
the extrapolation value xQ k=k1 and the correction term Kk ẽk . The extrapolation
(prediction) value is calculated by multiplying the estimation value found in the
previous step by the transfer matrix of the system. After that, correction is made
to the extrapolation value. Hence, the Kalman filter works on the principle of
correction of the prediction.
58 4 Estimation of the UAV Dynamics
y( k)
x(k 1)
xˆ ( k / k)
z( k / k 1)
x (k)
x ( k 1) ion
xˆ ( k 1/ k 1) polat y( k 1)
Extra
y( k 2)
K ( k) z ( k/ k 1)
x ( k 2) xˆ ( k / k 1)
y ( k 1)
Discrete Time
k 2 k 1 k k 1
Fig. 4.2 Time-dependent diagram of the value estimation mechanism in the Kalman filter
This time evaluation operation is combined with the identification operation of the
parameters and/or the model structure.
Nonlinear equations must be linearized first.
As seen above, the Kalman filter is a time-domain filter, because the design is done
in the time domain rather than the frequency domain. One of the advantages of the
Kalman filter is its ability to estimate time-variable parameters. When the Kalman
filter is used, the convergence of estimated values to the actual values of parameters
depends on the stability condition of the filter. The stability of conventional digital
filters is easily analyzed with z-transform methods [7]. We shall do the same for the
optimum Kalman filter.
It is assumed that the optimum discrete Kalman filter [Eqs. (4.3), (4.4), (4.7),
(4.8), (4.9), (4.10), and (4.11)] has reached a constant-gain condition. The estima-
tion equation (4.10) can also be expressed for convenience:
˚
b
xk=k D xQ k=k1 C Kk zk Hk xQ k=k1 (4.25)
xk=k D fFk Kk Hk Fk gb
b xk1=k1 C Kk zk (4.26)
bz
Xk=k D fFk Kk Hk Fk g z1bz
Xk=k C Kk Zkz (4.27)
ŒzI fFk Kk Hk Fk g bz
Xk=k D zKk Zkz (4.28)
where z denotes the usual z-transform variable and Z zk refers to the z-transform of
the measurement vector.
It is well known from linear system theory that the bracketed quantity on the left
side of Eq. (4.28) describes the natural modes of the system. The determinant of the
bracketed n n matrix gives the characteristic polynomial for the system; that is:
The roots of this polynomial provide information about the filter’s stability. If
all the roots lie inside the unit circle in the z-plane, the filter is stable; conversely,
if any root lies on or outside the unit circle, the filter is unstable. As a matter
of terminology, the roots of the characteristic polynomial are the same as the
eigenvalues of fFk Kk Hk Fk g.
When the Kalman filter algorithm is built for UAV state estimation, it is used for the
combined longitudinal and lateral dynamics of aircraft. Hence, the state vector to be
estimated is formed of n D 9 states as:
T
x D u w q h ˇ p r ; (4.30)
zk D Hk xk C vk ; (4.33)
E vk vjT D Rk ıkj ; (4.34)
E wk vjT D 0: (4.35)
4.5 Simulations 61
4.5 Simulations
In this section, simulation results for the state estimation of the Zagi UAV are given.
Simulation is realized in 500 steps for a period of 50 s with 0.1 s of sampling
time, t. For a joined state vector of longitudinal and lateral equations (4.30), the
measurement noise v is taken as:
T
v D Nr 0:92 0:92 0:083 0:17 0:2 0:05 0:083 0:083 0:17 (4.36)
where Nr is a random number with zero mean and unit variance Gaussian distri-
bution. Afterwards, R, the covariance of measurement noise, a 9 9 matrix for all
steps, is given as:
2 3
0:922 0
6 :: : : :: 7
RD4 : : : 5 (4.37)
2
0 0:17
For simulating measurements and also introducing system dynamics into the
filter, the equations of motion of the UAV are put into discretized form. If the state
space form of equations is shown as:
XP D AX C BU (4.38)
Here, A is the system matrix, B is the control distribution matrix, X is the state vector
T
as given above, and U is the control vector, given as u D ıe ıT ıa ır .
During simulations, the initial values of all states are assumed to be 0. On the
T
other hand, the control vector is chosen as u D 0:2 0 0 0:2 rad.
Measurements are simulated by adding noise values to the real values of dynamic
variations as:
Besides for showing the drawbacks of the OKF and the necessity for filter
adaptation in certain cases, errors are implemented in the measurements every 100th
step. These abrupt errors, which represent a failure condition at the measurement
channel, are formed by adding a constant term to the measurements of two states,
and ˇ.
Here, the simulation results are given for two states, u and in Figs. 4.3, 4.4,
4.5, 4.6, and 4.7. The results given in Figs. 4.5, 4.6, and 4.7 are for OKF in case
of faulty measurements. Also Table 4.2 presents the estimation errors for all states
at several time steps in case of faulty measurements. For any state, three graphs
are introduced in a combined fashion: a plot that shows the Kalman estimation
and simulated measurement values, a plot that gives the error of measurement
and Kalman estimation according to the real value found by discretized dynamic
equations (Eq. 4.39), and a plot that shows the variation of variance of Kalman
estimations.
The simulation results show that the OKF gives sufficiently good estimation
results in the case of normal operation of the estimation system (without faults).
However, the optimal filter fails in providing accurate estimations in case of
measurement faults, as expected. That signifies the importance of Kalman filter
adaptation, which will be discussed in the next section.
-5
-10
0 5 10 15 20 25 30 35 40 45 50
-5
0 5 10 15 20 25 30 35 40 45 50
1
Kalman Estimation Variance
variance
0.5
0
0 5 10 15 20 25 30 35 40 45 50
time(sec)
1 Measurement Simulation
-1
0 5 10 15 20 25 30 35 40 45 50
1
Kalman Estimation Error
error(rad)
-0.5
0 5 10 15 20 25 30 35 40 45 50
0.03
Kalman Estimation Variance
variance
0.02
0.01
0
0 5 10 15 20 25 30 35 40 45 50
time(sec)
The Kalman filter examined in the previous section is valid in the presence of
whole prior statistical data. In practice, situations where a priori data are known
with estimations or which exist completely are rarely encountered. For this case,
the developed algorithms are not optimal and values estimated via these algorithms
may not converge. Therefore, another approach is required to synthesize estimation
algorithms.
A priori uncertainty degree can vary:
1. Exact a priori statistical uncertainty – In this case, the shape and parameters of
the probability distribution rules of the components of measured and estimated
random processes are not understood. However, allowable limited change areas
of appropriate components of random processes are given. In such a case, the
synthesis of estimation algorithms can be performed only on the basis of a
guaranteed approach [8].
64 4 Estimation of the UAV Dynamics
Measurement Simulation
-10
0 5 10 15 20 25 30 35 40 45 50
10
Kalman Estimation Error
error(m/s)
-10
-20
0 5 10 15 20 25 30 35 40 45 50
1
Kalman Estimation Variance
variance
0.5
0
0 5 10 15 20 25 30 35 40 45 50
time(sec)
2. Partial prior statistical uncertainty – In this case, the probability distribution rules
of some parts of measured and estimated random processes’ components are
known. The number of parameters which have unknown probability characteris-
tics should not be too many. Increment in the number of these parameters causes
the quality of the problem solution to decrease. On account of this, when there
is a priori uncertainty in parameters, a distributions group is given instead of the
probability distribution rule of random processes. The estimation algorithm from
the distributions group must be chosen to provide optimal criterion. It means that
the estimation algorithm should be adaptive.
There are three approaches for the adaptation problem: parametric approach,
invariant principle-based approach, and structural approach. The parametric
approach is the most commonly used. The parametric adaptive estimation algorithm
is an algorithm that estimates the required components of random processes
at the end of the measurement data transaction and restores a priori statistical
characteristics of the dynamical system and measurements.
A self-tuning circuit is added to an ordinary Kalman filter in most parametric
adaptive estimation algorithms. When the Kalman filter designed, it is supposed
that the statistical characteristics of the system model, measurements, and noise
4.6 Necessity for Kalman Filter Adaptation 65
20
Measurement Simulation
0
-20
0 5 10 15 20 25 30 35 40 45 50
20
Kalman Estimation Error
error(rad)
0
Simulated Measuremnt Error
-20
-40
0 5 10 15 20 25 30 35 40 45 50
0.03
Kalman Estimation Variance
variance
0.02
0.01
0
0 5 10 15 20 25 30 35 40 45 50
time(sec)
20
15
teta(rad)
10
0
8 8.5 9 9.5 10 10.5 11 11.5 12
time(sec)
Fig. 4.7 Focused graph for estimation of the OKF in case of faulty measurements
66 4 Estimation of the UAV Dynamics
are known as being true. However, this supposition is not often accomplished, so
the applied filter is suboptimal. Another reason for the suboptimality of the filter
is simplifying the process calculations in the filter algorithm. In some cases, the
problem of concretizing the fault caused by suboptimality appears. For this reason,
noise covariances are estimated via the Kalman filter. This type of filter is called an
adaptive filter.
Some approaches to designing adaptive filters are presented below:
• Multiple model-based adaptive estimation (MMAE). In the MMAE approach,
adaptive filtering of a bank of Kalman filters runs in parallel under different
models for the statistical filter information matrices, i.e., the process noise matrix
Q and/or the update measurement noise matrix R.
• Unknown noise covariances of the Kalman filter are determined by the statistical
analysis of the innovation or residual series. This might be either as a direct
estimation of the covariance matrices or as an adaptation performed on the basis
of covariance scaling [9] (Fig. 4.8).
• Noise estimation is performed via probability methods. Estimation values of
unknown noise covariances are periodically used to renew the noise structure
in filtering algorithms.
• By processing the same values via an iteration procedure, the determination of
values of unknown covariances is based on the analysis of residual faults from
the previous data processing circuit. This method is generally used when huge
computational burden is required for a test.
Next, we review some of the existing common direct estimation techniques
for Kalman filter adaptation. Newly proposed scaling-based techniques will be
discussed in the next two chapters in detail. Readers are also referred to [9] for a
comparison between these different methods for Kalman filter adaptation.
4.6 Necessity for Kalman Filter Adaptation 67
Fig. 4.8 The Kalman filter (KF) adaptation by the statistical analysis of the innovation and
residual series: (a) innovation-based KF adaptation, (b) residual-based KF adaptation. For both
adaptation techniques, updating of the process and measurement noise covariances might be either
by estimation or by scaling
1 X
k
^
C ek D eQ j eQ T (4.41)
N jDkNC1 j
^
b
Rk D Cek Hk Pk=k1 HkT (4.42)
^
b
Qk D Kk Cek KkT (4.43)
68 4 Estimation of the UAV Dynamics
In the residual-based approach [16, 17], the measurement covariance matrix and/or
the system noise that are affected by the change in residual series are directly
updated.
The residual series can be described by:
k D zk Hb
xk=k ; (4.44)
1 X
k
b
Ck D j T (4.45)
N jDkNC1 j
b
Rk D b
Ck C Hk Pk=k HkT (4.46)
where Pk1=k1 is the covariance matrix of the estimation error in the previous step,
F is the system transformation matrix, and xk is the state correction series (the
difference between the filter estimation and the extrapolation estimation):
xk D b
xk=k b
xk=k1 (4.48)
In a stable state, only the first part of Eq. (4.47) is considered, in this case, taking
into account that xk D Kk k in Eq. (4.47) can be approximated by Eq. (4.43).
It is understood that, for some practical applications [(like an integrated inertial
navigation system (INS)/Global Positioning System (GPS)] [15], residual-based
adaptive estimation methods are more effective than innovation-based approaches.
Some of the disadvantages of the described adaptive estimation methods are
given below:
(a) In the estimation of unknown system and measurement noise covariance
matrices using the residual- or innovation-based approaches, the innovation
References 69
and residual vectors should be used for N cycles. This consumes a vast
amount of computing power and makes determining appropriate sliding window
dimension values a must.
(b) The measurement quantities, types, and distributions must be consistent in
order to be able to use the residual- or innovation-based estimation methods.
Otherwise, the covariance matrices of the system and measurement noise could
not be estimated using the residual or innovation matrices.
(c) Innovation-based estimators may cause a negative covariance matrix if the value
^
on the right-hand side of the formula b
Rk D Cek Hk Pk=k1 HkT is bigger than that
on the left-hand side.
4.7 Conclusion
In this chapter, an OKF algorithm for the state estimation of a specific UAV’s
dynamics is presented and examined. Simulations realized by the use of a MATLAB
program shows that the OKF algorithm gives accurate estimations under normal
operation of the system (without faults). The OKF algorithm fails when faults occur
in the estimation system.
Mathematically, the Kalman filter is a system of first-order ordinary differential
equations with quadratic nonlinearities, which are solved on digital computers.
When implementing navigational and other algorithms on a maintenance digital
computer, incorrect results may occur due to equipment failure, malfunctioning
of the computer, and noise during information transmission. For these cases, the
adaptive Kalman filter can be used.
In accordance to the main motivation behind the study, supplying reliable
parameter estimations to the control system of an autonomous UAV and ensuring
that it completes its mission successfully is the objectives of this monograph.
References
1. Kalman RE (1960) A new approach to linear filtering and prediction problems. Trans ASME J
Basic Eng 82:35–45
2. Hajiyev C, Caliskan F (2003) Fault diagnosis and reconfiguration in flight control systems.
Kluwer Academic Publishers, Boston
3. Zarchan P, Musoff H (2000) Fundamentals of Kalman filtering: a practical approach. Progress
in aeronautics and astronautics. AIAA, Reston
4. Sage AP, Melsa JL (1971) Estimation theory with applications in communication and control.
McGraw-Hill, New York
5. Hajiyev CM (1999) Radio navigation. Istanbul Technical University, Istanbul (in Turkish)
6. Anderson BDO, Moore JB (1979) Optimal filtering. Prentice Hall, Englewood Cliffs
7. Brown RG, Hwang PYC (1997) Introduction to random signals and applied Kalman filtering
with MATLAB exercises and solutions, 3rd edn. Wiley, New York
70 4 Estimation of the UAV Dynamics
5.1 Introduction
When there are sufficient measurements, the necessary states of an unmanned aerial
vehicle (UAV) can be estimated via a Kalman Filter. That is a desired procedure,
since it is important to know precisely parameters like velocity, altitude, attitude,
etc. When these states of the UAV are obtained without any problems, control of the
aircraft can be achieved successfully. On the other hand, it is a process dependent
upon the accuracy of measurements. If the measurements are not reliable due to any
kind of malfunction in the estimation system, the filter gives inaccurate results and
diverges over time. Since achieving fault tolerance in the design of a UAV flight
control system is important, the filter should be built robustly in order to avoid such
problems.
The Kalman filter approach to state estimation is quite sensitive to any mea-
surement malfunctions (abnormal measurements, sudden shifts in the measurement
channel, and other difficulties, such as decrease of instrument accuracy, an increase
in background noise, etc.). If the conditions of operation of the measurement
system do not correspond to the models used in the synthesis of the filter, then
these changes resulting from some possible failures at the measurement channels
significantly decrease the effectiveness of the estimation systems. When dealing
with the measurement faults in previous estimation steps, rather than the current
one, adaptive Kalman filters can be used to recover from the possible malfunctions.
A Kalman filter can be made adaptive and, hence, insensitive to the a priori
measurements or system uncertainties by using various different techniques. The
basic approaches to the adaptive Kalman filtering problem are multiple model-
based adaptive estimation (MMAE) [1–3], innovation-based adaptive estimation
(IAE) [4–6], and residual-based adaptive estimation (RAE) [7, 8]. While in the
first approach a bank of Kalman filters run in parallel under different models for
the filter’s statistical information, in the other approaches, the adaptation is done
directly to the covariance matrices of the measurement and/or system noises based
on the changes in the innovation or residual sequences.
In the methods described in [1–3], the faults are assumed to be known and the
Kalman filters are designed for the known sensor/actuator faults. As the MMAE
approach requires several parallel Kalman filters, and the faults should be known, it
can only be used in a limited number of applications.
Estimation of the covariance matrices by IAE and RAE requires the usage of
innovation vectors or residual vectors of m epoch. This increases the storage burden
and presents the determination of the width of the moving window, m, as another
problem. Furthermore, IAE and RAE estimators require that the number, type, and
distribution of measurements for all epochs within a window should be consistent. If
they are not, the covariance matrices of the measurement noises cannot be estimated
based on the innovation or residual vectors.
The adaptive Kalman filter presented in [9] has been applied to fuse position
signals from the Global Positioning System (GPS) and an inertial navigation system
(INS) for autonomous mobile vehicles. The extended Kalman filter (EKF) and the
noise characteristics have been modified using the fuzzy logic adaptive system. In
[10], a method for multisensor data fusion based on the adaptive fuzzy Kalman
filter is presented. This method is applied for fusing the position and orientation
signals from the dead reckoning (DR) system and the GPS for landing vehicle
navigation. The EKF and the characteristics of the measurement noise are modified
using the fuzzy adaptive system, which is based on a covariance matching technique.
It has been demonstrated that the fuzzy adaptive Kalman filter gives better results
(more accurate) than the EKF [9, 10]. In [11], a fuzzy logic-based adaptive
Kalman filter is used to build adaptive centralized, decentralized, and federated
Kalman filters for adaptive multisensor data fusion. The adaptation carried out
is in the sense of adaptively adjusting the measurement noise covariance matrix
of each local filter to fit the actual statistics of the noise profiles present in the
incoming measured data. A fuzzy inference system based on a covariance matching
technique is used as the adaptation mechanism in the paper. The simulation results
show that the architectures proposed by the authors are effective in situations
where there are several sensors measuring the same parameters, but each one has
different measurement dynamics and noise statistics. Although fuzzy logic-based
adaptive Kalman filter algorithms perform well under specific circumstances, they
are knowledge-based systems operating on linguistic variables, and these methods,
which are based on human experiences, are not widely applicable to vital systems
such as aircraft flight control systems.
Another concept is to scale the noise covariance matrix by multiplying it with a
time-dependent variable. One of the methods for constructing such an algorithm
is to use a single adaptive factor as a multiplier to the process or measurement
noise covariance matrices [12–15]. This algorithm, which may be called an adaptive
fading Kalman filter (AFKF), can be used when information about the dynamic
process or the a priori measurements are absent [16]. However, when the point
at issue is the current measurements, another technique to scale the measurement
noise covariance matrix and make the filter robust (insensitive to measurement
5.2 RKF with a Single Measurement Noise Scale Factor 73
the filter must be run robustly. Here, tr fg denotes the trace of the related matrix.
In order to determine the scale factor, Sk , let us substitute Eq. (5.1) into Eq. (5.3)
and use as input the scale factor calculated at the point where condition (5.3) is
satisfied:
˚ ˚
tr eQ k eQ Tk D tr Hk Pk=k1 HkT C Sk tr fRk g : (5.4)
˚
Then, in light of the equality tr eQ k eQ Tk D eQ Tk eQ k , Sk can be written as:
˚
eQ Tk eQ k tr Hk Pk=k1 HkT
Sk D : (5.5)
tr fRk g
the threshold value,
2˛,M , can be found. Hence, when the hypothesis 1 is correct,
the statistical value of ˇ k will be greater than the threshold value
2˛,M , i.e.:
0 W ˇk
2˛;M 8k
(5.8)
1 W ˇk >
2˛;M 9k:
1 X
k
eQ j eQ T D Hk Pk=k1 HkT C Sk Rk ; (5.9)
jDkC1 j
the use of Eq. (5.10), may not be diagonal and may have diagonal elements which
are “negative” or less than “one” (actually, that is physically impossible).
Therefore, in order to avoid such situations, it is suggested that the scale matrix
is composed by the following rule:
S D diag s1 ; s2 ; : : : ; sn ; (5.11)
where:
Here, Sii represents the ith diagonal element of the matrix S. Apart from that point, if
the measurements are faulty, Sk will change and, so, affect the Kalman gain matrix:
1
Kk D Pk=k1 HkT Hk Pk=k1 HkT C Sk Rk : (5.13)
The R-adaptation procedures with single and multiple measurement noise scale
factors are compared through simulations. Simulations are realized in 1,000 steps
for a period of 100 s with 0.1 s of sampling time, t. As an experimental platform,
the Zagi UAV is chosen and Kalman filter applications are performed with respect
to the dynamics and characteristics of this UAV.
During simulations, for testing the RKF algorithm, four kinds of measurement
malfunction scenarios are considered: instantaneous abnormal measurements, con-
tinuous bias, measurement noise increment, and fault of zero output. Also, in case
of measurement faults, simulations are also performed for the OKF to compare the
results with the RKF algorithms and understand the efficiency of the RKF in more
detail. Nonetheless,
2˛,s is taken as 16.919, and this value comes from a Chi-square
distribution when the number of degrees of freedom is 9 and the reliability level is
95 %.
5.4 Comparison of the R-adaptation Techniques 77
In this chapter, the first set of figures gives the OKF or RKF state estimation
results and the actual values are compared. The second set of figures shows the
error of the estimation process based on the actual state values of the UAV. The last
set of figures indicates the variance of the estimation that corresponds to the related
diagonal element of the covariance matrix of estimation errors.
Fig. 5.1 Estimation result for the pitch angle by the regular optimal Kalman filter (OKF) in case
of instantaneous abnormal measurements
78 5 Estimation of the UAV Dynamics in the Presence of Sensor Faults
Fig. 5.2 Estimation result for the pitch angle by the robust Kalman filter (RKF) with a single
measurement noise scale factor (SMNSF) in case of instantaneous abnormal measurements
The superiority of the RKF algorithms may also be seen in Table 5.1, which gives
the absolute estimation errors of the filters for the indicated time step. Note that, for
all the tables given hereafter, the highlighted results are obtained at time steps where
abnormal measurements are implemented.
When the measurements are faulty, the RKF with SMNSF compensates by
increasing its single scale factor and all of the measurements are taken into
consideration with a small weight for these time steps. Besides, the RKF with
MMNSF ensures the robustness of the filter by increasing related scale factors of the
matrix individually. Increment of the related scale factors brings about a decrement
in the related components of the Kalman gain, so as to reduce the corrective effect
of the innovation sequences of the faulty measurements on the state estimation
process. Progress may be understood better by examining the scale factor of the
RKF with SMNSF and one of the scalar measures (trace, determinant, and the
maximum eigenvalue) of the scale matrix of the RKF with MMNSF. The trace of
the scale matrix is examined in Table 5.2 as a scalar measure. Similar results have
been obtained when the fault is implemented in the other measurement channels.
5.4 Comparison of the R-adaptation Techniques 79
Fig. 5.3 Estimation result for the pitch angle by the RKF with multiple measurement noise scale
factors (MMNSF) in case of instantaneous abnormal measurements
Even though both of the RKFs give similar estimation results in case of
instantaneous abnormal measurements, the RKF with MMNSF can be thought
of as a more advantageous algorithm, as it takes the faulty measurements into
account individually. Disregarding all of the measurements, as the RKF with
SMNSF does, affects the estimation procedure of all the states. However, since
the abnormal measurement is implemented in only one measurement channel,
leaving out the related measurement would be more significant, especially when
the filter mathematical model of the UAV, which is constructed with the decoupled
longitudinal and lateral dynamics, is taken into consideration.
The continuous bias term is formed by adding a constant term to the measurement
of pitch angle, , in between the 30th and 60th seconds. As Fig. 5.4 shows,
again, the OKF fails in estimating states accurately. In contrast, the RKF with
SMNSF derogates estimation errors to a perceptible degree and gives relatively
better estimation outputs (Fig. 5.5). In this case, the RKF with MMNSF gives
80 5 Estimation of the UAV Dynamics in the Presence of Sensor Faults
Table 5.1 Comparison of absolute estimation errors of regular the optimal Kalman filter (OKF),
robust Kalman filter (RKF) with a single measurement noise scale factor (SMNSF), and RKF
with multiple measurement noise scale factors (MMNSF) in case of instantaneous abnormal
measurements
Parameter Absolute values of Absolute values of Absolute values of
error for regular error for RKF with error for RKF with
OKF SMNSF MMNSF
40 s 70 s 40 s 70 s 40 s 70 s
Δu (m / s) 3.4509 0.0609 0.0738 0.1009 0.0746 0.0093
Δw(m / s) 1.0658 0.0715 0.0419 0.177 0.0148 0.009
Δq (rad /s ) 0.1541 0.0402 0.0011 0.1155 0.0905 0.0407
Δθ (rad ) 2.6624 0.0139 0.003 0.0504 0.0165 0.0073
Δh(m) 2.9367 0.0266 0.0203 0.1519 0.0492 0.0075
Δβ (rad ) 0.0318 0.0283 0.0448 0.0015 0.0137 0.0145
Δp (rad / s ) 0.006 0.0009 0.0047 0.0193 0.0894 0.0833
Δr (rad / s ) 0.0085 0.0587 0.0424 0.0241 0.0195 0.0049
Δφ (rad ) 0.04 0.0931 0.1455 0.0156 0.0037 0.0606
Parameter 20 s 40 s 45 s 50 s 60 s 90 s
Scale factor of 346 345 1 1 343 1
RKF with
SMNSF
Sk
superior estimation results to the RKF with SMNSF by reducing the effect of the
innovation term of measurement and totally eliminating the estimation error
caused by the biased measurements of this channel (Fig. 5.6). As it disregards
only the related measurements, by the use of the innovative terms taken from all
other sensors operating correctly, it ensures accurate estimation outputs throughout
this period. It is possible to obtain the same results when the continuous bias is
implemented in the other measurement channels.
Comparison of the performance of the regular OKF, RKF with SMNSF, and
RKF with MMNSF in case of continuous bias of measurement can be seen
in Table 5.3.
5.4 Comparison of the R-adaptation Techniques 81
Fig. 5.4 Estimation result for the pitch angle by the regular OKF in case of continuous bias
In this table, the time means of absolute error values are given and they are
calculated for 30 s (between the 30th and 60th seconds for the fault case and between
the 65th and 95th seconds for the no fault case).
Again, it is possible to investigate the behavior of scale factors of both RKFs
(Table 5.4).
Fig. 5.5 Estimation result for the pitch angle by the RKF with SMNSF in case of continuous
bias
all measurements at the same time. Simulations give similar results when the noise
increment is implemented in the other measurement channels.
Table 5.5 shows a comparison of the mean of absolute estimation errors of the
regular OKF, RKF with SMNSF, and RKF with MMNSF in case of measurement
noise increment.
Table 5.6 presents the variation of the scale factor and the scalar measures of the
scale matrix.
In the last measurement malfunction scenario, it is assumed that height, h, cannot
be measured and the related sensor gives “0” as the output. This sort of fault is
easily simulated by taking the h measurement as “0” for the filter algorithm
in between the 30th and 60th seconds. As shown in Fig. 5.10, the OKF fails in
estimating the height and gives a result of “0” in light of the measurement. Similarly,
despite working robustly at the initial steps of the measurement fault implementation
5.4 Comparison of the R-adaptation Techniques 83
Fig. 5.6 Estimation result for the pitch angle by the RKF with MMNSF in case of continuous
bias
Table 5.3 Comparison of the mean absolute estimation errors of the regular OKF, RKF with
SMNSF, and RKF with MMNSF in case of continuous bias of measurements (for 30 s)
Parameter Absolute values of Absolute values of Absolute values of
error for regular error for RKF with error for RKF with
OKF SMNSF MMNSF
Fault No fault Fault No fault Fault No fault
Δu (m / s) 23.091 0.0906 1.9103 0.0982 0.0846 0.1066
Δw(m /s) 6.9494 0.0605 0.3337 0.0603 0.0637 0.06
Δq (rad /s ) 0.7806 0.0535 0.2153 0.0518 0.0512 0.0515
Δθ (rad ) 6.0079 0.027 0.1473 0.0261 0.0225 0.0276
Δh(m) 12.9792 0.0792 2.6621 0.0855 0.082 0.086
Δβ (rad ) 0.0248 0.0227 0.0083 0.0212 0.023 0.0215
Δp (rad /s ) 0.0314 0.0327 0.0016 0.0311 0.0267 0.0332
Δr (rad / s ) 0.0238 0.0255 0.0078 0.0233 0.0219 0.0249
Δφ (rad ) 0.0424 0.0441 0.009 0.0375 0.0327 0.0403
84 5 Estimation of the UAV Dynamics in the Presence of Sensor Faults
Fig. 5.7 Estimation result for the velocity component u by the regular OKF in case of measure-
ment noise increment
(between the 30th and 40th seconds approximately), the RKF with SMNSF accepts
the h measurement as a true value and estimates it as “0” (Fig. 5.11).
Clearly, not taking healthy measurements into consideration throughout this
period affects the estimation procedure of the RKF with SMNSF. However, the RKF
5.5 Remark on Stability 85
Fig. 5.8 Estimation result for the velocity component u by the RKF with SMNSF in case of
measurement noise increment
with MMNSF performs well under these circumstances and, by disregarding only
the faulty h measurement and using innovative terms taken from all other sensors
operating correctly, accurately estimates all of the states (Fig. 5.12). This behavior
can also be seen from Table 5.7.
Nonetheless, the behavior of the RKFs can be observed by investigating the
variation of scale factor and the scalar measures of the scale matrix (Table 5.8).
As can be seen, the scale factor takes a higher value than one only at the initial
steps of the malfunction scenario, while the matrix scalar measures conserve that
increment for the whole process.
Indeed, the proposed RKFs are not so different from the linear Kalman filter from
the point of view of structure; they may be assumed to be a modification. In
this sense, a similar approach to that of a linear Kalman filter may be followed
86 5 Estimation of the UAV Dynamics in the Presence of Sensor Faults
Fig. 5.9 Estimation result for the velocity component u by the RKF with MMNSF in case of
measurement noise increment
Table 5.5 Comparison of the mean absolute estimation errors of the regular OKF, RKF with
SMNSF, and RKF with MMNSF in case of measurement noise increment (for 30 s)
Parameter Absolute values of Absolute values of Absolute values of
error for regular error for RKF with error for RKF with
OKF SMNSF MMNSF
Fault No fault Fault No fault Fault No fault
Δu (m / s) 4.0183 0.0901 0.1571 0.0959 0.0836 0.112
Δw(m / s) 1.2699 0.0656 0.0302 0.0644 0.052 0.0656
Δq (rad / s ) 0.1462 0.056 0.0197 0.0562 0.0524 0.0557
Δθ (rad ) 0.4319 0.0287 0.0177 0.0266 0.0243 0.0275
Δh(m) 1.336 0.0869 0.1411 0.0913 0.0818 0.0905
Δβ (rad ) 0.0227 0.0226 0.0328 0.0222 0.0225 0.0224
Δp (rad / s ) 0.0295 0.0314 0.0125 0.0308 0.0283 0.0316
Δr (rad / s ) 0.0233 0.0231 0.03 0.0214 0.0198 0.0221
Δφ (rad ) 0.0455 0.0344 0.0443 0.0355 0.0388 0.0395
5.5 Remark on Stability 87
Parameter 20 s 40 s 45 s 50 s 60 s 90 s
Fig. 5.10 Estimation result for height h by the regular OKF in case of fault of zero output
for the stability analysis of the RKFs. From [21], it is known that the following
characteristic polynomial of the system can be used for stability analyses of a linear
Kalman filter:
Fig. 5.11 Estimation result for height h by the RKF with SMNSF in case of fault of zero output
The roots of this polynomial provide information about the filter stability. If all
the roots lie inside the unit circle in the z-plane, the filter is stable; conversely,
if any root lies on or outside the unit circle, the filter is unstable. As a matter
of terminology, the roots of the characteristic polynomial are the same as the
eigenvalues of:
ŒfFk Kk Hk Fk g (5.15)
Nonetheless, we assume that a similar investigation is achieved for the RKFs. Let’s
examine the case for the RKF with MMNSF. Hence, if we substitute the necessary
equality instead of Kalman filter gain, Eq. (5.15) becomes:
For the RKF with SMNSF:
n h 1 i o
Fk Pk=k1 HkT Hk Pk=k1 HkT C Sk Rk Hk F k (5.16)
5.5 Remark on Stability 89
Fig. 5.12 Estimation result for height h by the RKF with MMNSF in case of fault of zero output
Table 5.7 Comparison of the mean absolute estimation errors of the regular OKF, RKF with
SMNSF, and RKF with MMNSF in case of fault of zero output (for 30 s)
Parameter Absolute values of Absolute values of Absolute values of
error for regular error for RKF with error for RKF with
OKF SMNSF MMNSF
Fault No fault Fault No fault Fault No fault
Δu (m /s) 3.45 0.1313 3.129 1.6054 0.1677 0.096
Δw(m /s) 1.1045 0.0709 0.744 0.3416 0.0621 0.0606
Δq (rad /s ) 0.1388 0.048 0.3018 0.2252 0.0554 0.0529
Δθ (rad ) 0.2646 0.0297 0.3057 0.2664 0.0357 0.0278
Δh(m) 80.984 0.0852 71.2394 17.6275 0.9192 0.0954
Δβ (rad ) 0.0223 0.0236 0.0306 0.019 0.0211 0.0246
Δp (rad / s ) 0.0317 0.0336 0.0245 0.0204 0.0316 0.0321
Δr (rad / s ) 0.0235 0.0215 0.0273 0.0214 0.0187 0.0221
Δφ (rad ) 0.0413 0.0384 0.041 0.0331 0.0404 0.041
90 5 Estimation of the UAV Dynamics in the Presence of Sensor Faults
Parameter 20 s 40 s 45 s 50 s 60 s 90 s
Now we should search for the roots of Eqs. (5.16) and (5.17). The realized
simulations for every malfunction case showed that both RKFs are stable in case
of a measurement malfunction. Nonetheless, the RKF with MMNSF is shown to be
more stable and, usually, only one root reaches near a point of the unit cycle limit.
In Fig. 5.13, the pole-zero map of the RKF with MMNSF is given.
The pole-zero map is obtained for the instantaneous abnormal measurement
malfunction case at the 20th second, where the malfunction is implemented. It is
not very different from that found for the OKF in the same sampling step. During
adaptation with MMNSF, related roots (that belong to the states, the innovation
channel of which is corrected via the scale matrix) moves through the unit cycle limit
but they do not approach the marginal stability border very closely. On the other
hand, for the RKF with SMNSF, where the whole innovation vector is disregarded in
case of malfunction, all of the roots move through the stability border (see Fig. 5.14)
and, usually, one of them takes a nearly marginal (jzj D 1) value, so the stability
characteristics worsen.
The simulation results indicate that the proposed RKFs are not at risk of con-
vergence; namely, they are stable. In order to guarantee the accuracy of simulation
outputs, some even harsher malfunction cases are tested (far greater malfunctions
are implemented). It is seen that the magnitude of the malfunction is not a factor
which may break down the stability of the filter and the roots of the RKFs stay
inside the unit cycle for every condition.
Moreover, the simulation results indicated the better stability characteristics of
the RKF with MMNSF compared to the RKF with SMNSF. For every malfunction
case, the roots of the RKF with MMNSF do not approach the stability limit (unit
cycle) as much those of the RKF with SMNSF.
5.6 Conclusion and Discussion 91
Fig. 5.13 Pole-zero map of the RKF with MMNSF for the case of instantaneous abnormal
measurement (at the 20th second)
In this chapter, RKF algorithms with filter gain correction for the case of mea-
surement malfunctions at the current estimation step have been developed. By the
use of defined variables named scale factors, faulty measurements are taken into
consideration with small weight and the estimations are corrected without affecting
the characteristics of the accurate measurements. In the presented RKFs, the filter
gain correction is performed only in the case of malfunctions in the measurement
system.
RKF algorithms with single and multiple MNSF are proposed. In the first case,
the filter is adapted by using single scale factor as a corrective term on the filter
gain, and, in the second, a scale matrix composed of multiple factors is used to fix
the relevant terms of the Kalman gain matrix individually.
The proposed RKF algorithms with single and multiple MNSF are applied for the
state estimation process of a UAV platform. Algorithms are tested for four different
measurement malfunction scenarios and the results are compared with the outputs
of the OKF for the same cases: instantaneous abnormal measurements, continuous
bias of measurements, measurement noise increment, and fault of zero output. In all
92 5 Estimation of the UAV Dynamics in the Presence of Sensor Faults
Fig. 5.14 Pole-zero map of the RKF with SMNSF for the case of instantaneous abnormal
measurement (at the 20th second)
four scenarios, the OKF becomes faulty while the introduced RKF algorithms stand
robust to the measurement faults. Comparison of the simulation results show that
the performance of the RKF with MMNSF is better than the RKF with SMNSF. The
advantage of the RKF with MMNSF can be especially seen in cases of continuous
bias of measurements and fault of zero sensor output. Nonetheless, simulations for
the stability analyses indicated that both RKF algorithms are stable in case of a
malfunction, but the stability characteristics of the RKF with MMNSF are better
than the RKF with SMNSF.
The proposed approach does not require a priori statistical characteristics of
the faults and can be used for both linear and nonlinear systems. Furthermore,
the presented RKF algorithms are easily implemented in practice and their com-
putational burden is not heavy. These characteristics make the introduced RKF
algorithms extremely important from the point of view of supplying reliable
parameter estimations for the flight control system of an autonomous UAV.
References 93
References
1. White NA, Maybeck PS, DeVilbiss SL (1998) Detection of interference/jamming and spoofing
in a DGPS-aided inertial system. IEEE Trans Aerosp Electron Syst 34(4):1208–1217. doi:10.
1109/7.722708
2. Zhang Y, Li XR (1997) Detection and diagnosis of sensor and actuator failures using interacting
multiple-model estimator. In: Proceedings of the 36th IEEE conference on decision and control,
San Diego, CA, December 1997, vol 5, pp 4475–4480. doi:10.1109/CDC.1997.649671
3. Maybeck PS (1999) Multiple model adaptive algorithms for detecting and compensating sensor
and actuator/surface failures in aircraft flight control systems. Int J Robust Nonlinear Control
9(14):1051–1070
4. Mehra RK (1970) On the identification of variances and adaptive Kalman filtering. IEEE Trans
Autom Control 15(2):175–184
5. Maybeck PS (1982) Stochastic models, estimation, and control, vol I and II. Academic, New
York
6. Salychez OS (1994) Special studies in dynamic estimation procedures with case studies
in inertial surveying, ENGO 699.26 lecture notes. Department of Geomatics Engineering,
University of Calgary, Calgary
7. Mohamed AH, Schwarz KP (1999) Adaptive Kalman filtering for INS/GPS. J Geod
73(4):193–203
8. Wang J, Stewart MP, Tsakiri M (2000) Adaptive Kalman filtering for integration of GPS with
GLONASS and INS. In: Proceedings of Geodesy Beyond 2000: the challenges of the first
decade, Birmingham, July 1999, pp 325–330
9. Sasiadek JZ, Wang Q (1999) Sensor fusion based on fuzzy Kalman filtering for autonomous
robot vehicle. In: Proceedings of the 1999 IEEE international conference on robotics and
automation, Detroit, MI, May 1999, vol 4, pp 2970–2975. doi:10.1109/ROBOT.1999.774048
10. Zhang ST, Wei XY (2003) Fuzzy adaptive Kalman filtering for DR/GPS. In: Proceedings of the
second international conference on machine learning and cybernetics, Xi-an, China, November
2003, vol 5, pp 2634–2637
11. Escamilla-Ambrosio PJ, Mort N (2003) Hybrid Kalman filter-fuzzy logic adaptive multisensor
data fusion architectures. In: Proceedings of the 42nd IEEE conference on decision and control,
Maui, HI, December 2003, vol 5, 5215–5220. doi:10.1109/CDC.2003.1272465
12. Hu C, Chen W, Chen Y, Liu D (2003) Adaptive Kalman filtering for vehicle navigation. J Glob
Position Syst 2(1):42–47
13. Hide C, Moore T, Smith M (2004) Adaptive Kalman filtering algorithms for integrating
GPS and low cost INS. In: Proceedings of the position location and navigation symposium
(PLANS), Monterey, CA, April 2004, pp 227–233
14. Ding W, Wang J, Rizos C, Kinlyside D (2007) Improving adaptive Kalman estimation in
GPS/INS integration. J Navig 60(3):517–529. doi:10.1017/S0373463307004316
15. Jwo DJ, Weng TP (2008) An adaptive sensor fusion method with applications in integrated
navigation. J Navig 61(4):705–721. doi:10.1017/S0373463308004827
16. Kim KH, Lee JG, Park CG (2006) Adaptive two-stage Kalman filter in the presence of
unknown random bias. Int J Adapt Contr Signal Process 20(7):305–319. doi:10.1002/acs.900
17. Hajiyev C (2007) Adaptive filtration algorithm with the filter gain correction applied to
integrated INS/radar altimeter. Proc Inst Mech Eng Part G J Aerosp Eng 221(5):847–885.
doi:10.1243/09544100JAERO173
18. Geng Y, Wang J (2008) Adaptive estimation of multiple fading factors in Kalman filter for
navigation applications. GPS Solutions 12(4):273–279. doi:10.1007/s10291-007-0084-6
19. Hajiyev Ch, Soken HE (2009) Adaptive Kalman filter with the filter gain correction applied
to UAV flight dynamics. In: Proceedings of the 17th Mediterranean conference on control and
automation, Thessaloniki, Greece, 24–26 June 2009
94 5 Estimation of the UAV Dynamics in the Presence of Sensor Faults
20. Soken HE, Hajiyev C (2010) An adaptive sensor fusion method applied to UAV dynamics. In:
Proceedings of the international unmanned vehicles workshop (UVW2010), Istanbul, Turkey,
10–12 June 2010, pp 57–62
21. Hajiyev C, Caliskan F (2003) Fault diagnosis and reconfiguration in flight control systems.
Kluwer Academic Publishers, Boston
22. Hajiyev Ch, Soken HE (2009) Adaptive Kalman filter with multiple fading factors for UAV
state estimation. In: Proceedings of the 7th IFAC international symposium on fault detection,
supervision and safety of technical processes (SafeProcess 2009), Barcelona, Spain, 30 June–3
July, pp 77–82
23. Hajiyev C, Soken HE (2010) Robust estimation of UAV dynamics in the presence of
measurement faults. J Aerosp Eng 25(1):80–89
Chapter 6
Estimation of Unmanned Aerial Vehicle
Dynamics in the Presence of Sensor/Actuator
Faults
6.1 Introduction
When there is an actuator fault in the system, it results in changes in the control
distribution matrix. In a similar manner to the R-adaptation procedure, in case
of an actuator fault, the real error of the innovation covariance will exceed the
theoretical covariance. Thus, the basic premise of the Q-adaptation is also to obtain
an appropriate multiplier matrix for the Q matrix such that the real and theoretical
values of the innovation covariance match. Hence, if a fading matrix (ƒk ) composed
of multiple fading factors is added into the algorithm as [14–17]:
1 X
k
eQ j eQ Tj D Hk Fk Pk1=k1 FkT C ƒk Gk Qk GTk HkT C Rk (6.1)
jDkC1
For a specific case where all the states are measured as .Hk D I/, as in the case
examined here, Eq. (6.2) reduces to:
0 1
1 X k
1
ƒk D @ eQ j eQ T Fk Pk1=k1 Fk T Rk A Gk Qk Gk T (6.3)
jDkC1 j
98 6 Estimation of the UAV Dynamics in the Presence of Sensor/Actuator Faults
where:
Here, ƒii represents the ith diagonal element of the matrix ƒ. Apart from that point,
if there is an actuator fault in the system, ƒk must be included in the estimation
process as:
Due to the scale or fading matrix, the covariance of the estimation error of the RAKF
increases in comparison with the OKF. Therefore, the RAKF algorithm is used only
when the there is a fault and, in all other cases, the procedure is run optimally with
the regular OKF. The fault detection process is applied by checking a statistical
function at each estimation step. At that point, the two following hypotheses may be
introduced:
• 0 ; the system is operating normally
• 1 ; there is a malfunction in the estimation system
In order to detect the faults, the statistical function used is defined as:
1
ˇk D eQ Tk Hk Pk=k1 HkT C Rk eQ k : (6.7)
then the threshold value,
2˛,M , can be found. Hence, when the hypothesis 1 is
correct, the statistical value of ˇ k will be greater than the threshold value
2˛,M , i.e.:
0 W ˇk
2˛;M 8k
(6.9)
1 W ˇk >
2˛;M 9k:
6.4 Numerical Example 99
On the other hand, after detecting the fault in the system, the key point is
detecting the type of fault (either a sensor or an actuator fault). After that, the
appropriate adaptation (R- or Q-adaptation) may be applied [15–17]. The fault
isolation can be realized by an algorithm similar to that proposed for aircraft in
[18]. When a regular Kalman filter is used, the decision statistics change regardless
of whether the fault is in the sensors or in the actuators. If an RKF based on the
Doyle–Stein condition is used, it is easy to distinguish between sensor and actuator
faults. The Kalman filter that satisfies the Doyle–Stein condition is referred in [18],
as the RKF is insensitive to actuator failures and may also be used for fault isolation
in our study. After that, the appropriate adaptation (R- or Q-adaptation) may be
applied. The Doyle–Stein condition is as follows:
Here, K is the filter gain, I is a unit matrix, H is the system measurement matrix, Bs
is the control distribution matrix in continuous time, and s D .sI Fs /1 , where
Fs is the system matrix in continuous time.
The RKF is very useful for isolating sensor and control surface failures, as it is
insensitive to the latter failures. If the Kalman filter process noise covariance matrix
is chosen as
then the filter becomes robust against the actuator faults. Here, Q is the process noise
covariance matrix for the nominal plant, qr is a parameter that approaches to infinity
as the filter’s robustness increases, and V is a positive-definite symmetric matrix.
If a sensor fault occurs, the R-adaptation is realized; else, the procedure is
continued with the Q-adaptation. The algorithm architecture given in Fig. 6.1
summarizes the integrated R- and Q-adaptation procedures. This new filter, where
the robust and adaptive Kalman filters are integrated, can be called the RAKF.
In order to examine the effectiveness of the proposed RAKF, it was applied for the
state estimation of the Zagi UAV. The simulations are realized in 1,000 steps for a
period of 100 s with 0.1 s of sampling time, t.
In case of sensor/actuator fault, the simulations are also carried out for the OKF
so as to compare the results with the RAKF algorithm and understand better the
efficiency of the proposed algorithm.
Nonetheless, for the fault detection process,
2˛,M is taken as 21.7, and this value
comes from the Chi-square distribution when the number of degrees of freedom is
9 and the reliability level is 99 %.
100 6 Estimation of the UAV Dynamics in the Presence of Sensor/Actuator Faults
Fig. 6.1 System architecture for the robust adaptive Kalman filter (RAKF) algorithm
All the figures given in this section consist of three subfigures. The top
figure gives the OKF or RAKF state estimation results and the actual values for
comparison. The middle figure shows the error of the estimation process based on
6.4 Numerical Example 101
the actual state values of the UAV. The bottom figure indicates the variance of the
estimation that is obtained as the diagonal elements of the Kalman filter estimation
error covariance matrix.
During the simulations, for testing the fault detection and isolation procedures
and also the RAKF algorithm, two different faults are implemented to the systems.
Firstly, a sensor fault is formed by adding a constant term to the measurement of the
pitch angle (the constant term is nearly of the same magnitude as the measurement,
so the measurement’s magnitude is doubled for those instants of time), , in
between the 30th and 50th seconds. Furthermore, in order to check that fault
isolation can be achieved successfully and switching from one adaptation stage to
the other is realized properly, as another fault, the actuator fault is simulated by
taking the first column elements of the longitudinal control distribution matrix to
be nearly zero in between the 80th and 90th seconds. The graph for the statistical
values of ˇ k in case of using the OKF is shown in Fig. 6.2. As can be seen, except
for the periods where the sensor and actuator failures occur, ˇ k is lower than the
threshold value.
On the contrary, when one of these failures occurs, ˇ k grows rapidly and exceeds
the threshold value. Hence, the 1 hypothesis is judged to be true and it is shown
that, by using such a procedure, the fault can be detected. The results for the fault
isolation process are given in Fig. 6.3.
Fig. 6.2 Behavior of the statistics ˇ k in case of measurement failure (in between the 30th and 50th
seconds) and actuator failure (in between the 80th and 90th seconds) when the optimal Kalman
filter (OKF) is used
102 6 Estimation of the UAV Dynamics in the Presence of Sensor/Actuator Faults
Fig. 6.3 Behavior of the statistics ˇ k in case of measurement failure (in between the 30th and 50th
seconds) and actuator failure (in between the 80th and 90th seconds) when the robust Kalman filter
(RKF), which is insensitive to actuator failures, is used
As mentioned, the RKF, which is insensitive to actuator failures, is used for this
isolation process. The behavior of the statistics ˇ k in the case of sensor/actuator
failures shows that detecting the actuator failure is not possible when the RKF is
used, since it is insensitive against such faults; on the other hand, the sensor failure
is detected immediately. Hence, the fault isolation can be realized by using such a
procedure. The RKF, which is insensitive to actuator failures and sensitive to sensor
failures, is employed to isolate the sensor and actuator failures in the proposed
RAKF algorithm (see Fig. 6.1).
The results of the first fault scenario realized for testing the RAKF algorithm in
the presence of sensor/actuator failures are given in Figs. 6.4 and 6.5. It is clear that
the RAKF gives superior estimation results compared to the OKF in both cases.
The type of fault is precisely detected and the RAKF tunes itself according to the
appropriate adaptation rule. The RAKF algorithm first detects the fault by using the
fault detection process given with Eqs. (6.7), (6.8), and (6.9). Then, it isolates the
fault using the RKF based on the Doyle–Stein condition, such that it distinguishes
the sensor and actuator faults. As the final step, it adapts the R matrix in between
the 30th and 50th seconds and the Q matrix in between the 80th and 90th seconds.
6.4 Numerical Example 103
Actual Value
0
-10
0 10 20 30 40 50 60 70 80 90 100
5
error(m/s)
-5
-10
0 10 20 30 40 50 60 70 80 90 100
variance [(m/s)2]
0.5
0
0 10 20 30 40 50 60 70 80 90 100
time(sec)
Fig. 6.4 u velocity estimation via the OKF in case of actuator and bias sensor faults
As seen from the figures, it ensures accurate estimation results for the whole process,
while the OKF fails in both the sensor fault and the actuator fault conditions. This
proves us three facts:
– The fault detection process is working correctly and the fault in the system is
detected without any time delay
– By using the RKF based on the Doyle–Stein condition, the type of fault is
distinguished and the RAKF decides to use the right adaptation procedure against
the fault
– Both of the proposed adaptation methods are functioning properly
Nevertheless, the performance of both filters, the OKF and the RAKF, in case of
fault can be understood more clearly by examining Table 6.1.
In the table, the mean absolute values of error, which are calculated for 10 s,
are given for all three conditions: sensor fault, actuator fault, and no fault. The
tabulated results show that the estimations by the OKF worsen for all of the
longitudinal states in case of sensor/actuator fault, while the RAKF retains accurate
104 6 Estimation of the UAV Dynamics in the Presence of Sensor/Actuator Faults
Actual Value
0
-10
0 10 20 30 40 50 60 70 80 90 100
5
error(m/s)
-5
-10
0 10 20 30 40 50 60 70 80 90 100
variance [(m/s)2]
0.5
0
0 10 20 30 40 50 60 70 80 90 100
time(sec)
Fig. 6.5 u velocity estimation via the RAKF in case of actuator and bias sensor faults
Table 6.1 Comparison of the mean absolute values of error for the optimal Kalman filter (OKF)
and the robust adaptive Kalman filter (RAKF) in case of actuator and bias sensor faults
Mean absolute value of error for regular
Parameter OKF Mean absolute value of error for RAKF
Operating Operating
Sensor fault normally Actuator fault Sensor fault normally Actuator fault
u(m/s) 1.8799 0.1891 0.8116 0.0528 0.0578 0.1815
w(m/s) 0.6762 0.0931 1.1939 0.0649 0.0611 0.2678
q(rad/s) 0.0604 0.0621 0.1602 0.0521 0.0542 0.0686
(rad) 0.6497 0.053 0.0796 0.0161 0.0132 0.0325
h(m) 0.5493 0.1059 0.1273 0.0644 0.0636 0.0932
ˇ(rad) 0.0319 0.0302 0.0338 0.0123 0.0155 0.0165
p(rad/s) 0.0514 0.0556 0.0538 0.0282 0.0257 0.0304
r(rad/s) 0.0508 0.045 0.0465 0.0115 0.0146 0.0163
(rad) 0.0758 0.0786 0.0805 0.0167 0.0192 0.0349
6.4 Numerical Example 105
Actual Value
0
-10
0 10 20 30 40 50 60 70 80 90 100
5
error(m/s)
-5
-10
0 10 20 30 40 50 60 70 80 90 100
variance [(m/s)2]
0.5
0
0 10 20 30 40 50 60 70 80 90 100
time(sec)
Fig. 6.6 u velocity estimation via the OKF in case of actuator and noise increment sensor faults
estimation behavior for the whole process. Because of the decoupled dynamics used
for the UAV, any fault in the longitudinal channel (e.g., sensor fault for the sensor
measuring or the actuator fault affecting the longitudinal control derivatives)
does not have an effect on the lateral states as expected. This fact makes the
importance of using multiple adaptive factors for both adaptation procedures more
distinct.
As the second testing scenario for the RAKF algorithm, this time, the sensor
fault is characterized by multiplying the variance of the noise of the measurement
for the u velocity component with a constant term in between the 30th and 50th
seconds, while the actuator fault is the same as the first scenario. As Figs. 6.6 and
6.7 show, again, the OKF outputs involve error, while the RAKF algorithm achieves
estimation of the states accurately. Nonetheless, another consequence of the second
fault scenario is the good performance of the fault isolation scheme; regardless of the
type of sensor fault, the isolation scheme can distinguish the type of fault effectively
and the RAKF works properly.
In order to clarify the results, a comparison of the mean absolute values of error
is given in Table 6.2.
106 6 Estimation of the UAV Dynamics in the Presence of Sensor/Actuator Faults
Actual Value
0
-10
0 10 20 30 40 50 60 70 80 90 100
5
error(m/s)
-5
-10
0 10 20 30 40 50 60 70 80 90 100
variance [(m/s)2]
0.5
0
0 10 20 30 40 50 60 70 80 90 100
time(sec)
Fig. 6.7 u velocity estimation via the RAKF in case of actuator and noise increment sensor faults
Table 6.2 Comparison of the mean absolute values of error for the OKF and RAKF in case of
actuator and noise increment sensor faults
Mean absolute value of error for regular
Parameter OKF Mean absolute value of error for RAKF
Operating Operating
Sensor fault normally Actuator fault Sensor fault normally Actuator fault
u(m/s) 0.6235 0.2317 0.8352 0.0819 0.0635 0.1303
w(m/s) 0.2588 0.1482 1.1963 0.0687 0.0626 0.2621
q(rad/s) 0.0611 0.0673 0.1505 0.055 0.0508 0.0711
(rad) 0.0653 0.0621 0.0688 0.0115 0.0155 0.0376
h(m) 0.1101 0.1 0.1323 0.0826 0.0733 0.0785
ˇ(rad) 0.0355 0.0392 0.0319 0.0163 0.0149 0.0134
p(rad/s) 0.0497 0.0465 0.0553 0.0308 0.026 0.0265
r(rad/s) 0.0456 0.0443 0.0449 0.0135 0.013 0.0112
(rad) 0.08 0.0694 0.0695 0.0262 0.0196 0.0359
References 107
In this chapter, first, the Q-adaptation with multiple adaptive factors for actuator
fault compensation is presented. Then, the RAKF with the R- and Q-adaptations
against sensor/actuator failures is proposed.
The proposed RAKF algorithm with the multiple adaptive factors incorporates
the R- and Q-adaptation procedures and adapts itself against sensor/actuator faults
in order to ensure robustness. The filter first detects the type of fault (either sensor
or actuator) and then isolates it. After that, the appropriate adaptation (R- or Q-
adaptation) is applied such that the estimation characteristic is not deteriorated.
Unlike existing studies, the proposed adaptation methods for both the R- and Q-
adaptations are simplistic and can be applied easily with a small modification to the
Kalman filter.
In this chapter, it has been shown how to incorporate these two different
adaptation processes such that the fault is isolated and, so, the RAKF algorithm
is a complete estimator including fault detection and isolation processes for the
case of sensor/actuator faults. When a regular Kalman filter is used, the decision
statistics change regardless of whether the fault is in the sensors or in the actuators,
whereas when an RKF based on the Doyle–Stein condition is used, it is easy to
distinguish the sensor and control surface/actuator faults. The RAKF is applied for
the state estimation procedure of a UAV and its effectiveness is investigated by
comparisons with the OKF in case of faults. The proposed RAKF algorithm may
play an important role for UAV control systems since it gives accurate estimation
results despite the sensor or actuator faults. Regarding the harsh environment where
UAVs are generally used, encountering any kind of fault has a high probability and,
so, preferring the proposed RAKF algorithm instead of the regular OKF may bring
a significant advantage.
References
7. Hajiyev Ch, Soken HE (2009) Adaptive Kalman filter with the filter gain correction applied
to UAV flight dynamics. In: Proceedings of the 17th Mediterranean conference on control and
automation, Thessaloniki, Greece, 24–26 June 2009
8. Hajiyev Ch, Soken HE (2009) Adaptive Kalman filter with multiple fading factors for UAV
state estimation. In: Proceedings of the 7th IFAC international symposium on fault detection,
supervision and safety of technical processes (SafeProcess 2009), Barcelona, Spain, 30 June–3
July, pp 77–82
9. Hajiyev C, Soken HE (2010) Robust estimation of UAV dynamics in the presence of
measurement faults. J Aerosp Eng 25(1):80–89
10. Soken HE, Hajiyev C (2010) Pico satellite attitude estimation via robust unscented Kalman
filter in the presence of measurement faults. ISA Trans 49(3):249–256. doi:10.1016/j.isatra.
2010.04.001
11. Jwo DJ, Weng TP (2008) An adaptive sensor fusion method with applications in integrated
navigation. J Navig 61(4):705–721. doi:10.1017/S0373463308004827
12. Jwo DJ, Chang SC (2009) Particle swarm optimization for GPS navigation Kalman filter
adaptation. Aircr Eng Aerosp Technol Int J 81:343–352. doi:10.1108/00022660910967336
13. Kim KH, Lee JG, Park CG (2006) Adaptive two-stage Kalman filter in the presence of
unknown random bias. Int J Adapt Control Signal Process 20(7):305–319. doi:10.1002/acs.
900
14. Hajiyev Ch, Soken HE (2011) Fault tolerant estimation of UAV dynamics via robust adaptive
Kalman filter. In: Proceedings of the special international conference on complex systems:
synergy of control, communications and computing (COSY 2011), Ohridi, Republic of
Macedonia, 16–20 September 2011, pp 311–320
15. Hajiyev C, Soken HE (2013) Robust adaptive Kalman filter for estimation of UAV dynamics
in the presence of sensor/actuator faults. Aerosp Sci Technol 28(1):376–383. doi:10.1016/j.ast.
2012.12.003
16. Hajiyev C, Soken HE (2013) Robust self-adaptive Kalman filter with the R and Q adaptations
against sensor/actuator failures. In: Zhang WJ (ed) Self-organization: theories and methods.
Nova Science Publishers, New York
17. Hajiyev C, Soken HE (2015) Fault tolerant estimation of UAV dynamics via robust adaptive
Kalman filter. In: Dimirovski GM (ed) Complex systems: relationships between control,
communications and computing. Springer (in press), Berlin
18. Hajiyev C, Caliskan F (2005) Sensor and control surface/actuator failure detection and isolation
applied to F-16 flight dynamics. Aircr Eng Aerosp Technol Int J 77:152–160
Chapter 7
Fault Detection, Isolation, and Data Fusion for
Unmanned Aerial Vehicle Air Data Systems
7.1 Introduction
For unmanned aerial vehicles (UAVs) to carry out their assigned tasks and complete
missions successfully, appropriate air measurement data should be derived in case
of fault occurrence in sensor measurements. In order to produce appropriate data,
the following should be done: detecting and isolating the faults in the measurement
system, determining the most appropriate data, using other measurement data for
estimation purposes, and deriving the most appropriate data using the available
measurements (decision algorithm). There have been many studies using different
techniques investigating the subject [1–10]. Some of them are given below.
In [1], a method using the extended Kalman filtering (EKF) structure is proposed
and the efficacy of the method is shown using a space vehicle dynamic model which
includes redundant sensors to measure the parameters. Although using multiple
extended Kalman filters together in nonlinear systems for fault detection is a known
method, generally, the initial value of the covariance matrix is chosen by expertise;
however, in this study, a new method for determining the optimal covariance matrix
is proposed.
Two different fusion methods for a UAV that integrate inertial navigation system
(INS) and accurate dynamic model measurements are proposed in [2]. In order
to reduce the growing errors of the INS and obtain better measurements, the
measurements of a dynamic model are integrated using EKF. In that study, the
errors between the INS and flight vehicle model (FVM) measurements are estimated
with the filtering technique using the error model of the system. The first proposed
method is based on comparing speed and Euler angle measurements produced by the
INS and the FVM. In this first method, inertial measurement unit (IMU) acceleration
and angular speed errors, INS speed and Euler angle errors, and FVM speed and
Euler angle errors are determined and used to obtain better measurement data. In the
second method, the estimated values of acceleration and angular speed taken from
the FVM are used for the calibration of the IMU. The errors in the acceleration and
angular speed measurements are estimated using an error model of the FVM.
In [3], a sensor fault detection method using data fusion techniques is inves-
tigated. In this study, the EKF method is used to integrate data. An adaptive
method which updates the sensor noise covariance matrix according to residual
values covariance is used in EKF. This new method is called the adaptive modified
extended Kalman filtering (AMEKF) algorithm. The AMEKF system estimates the
system noise and sensor noise covariance matrices, which can affect normal filter
working conditions, and finds calibration errors.
A weak model-based method for the detection of sensor abrupt faults is
considered in [4]. In that study, some characteristics of the sensors are used to derive
the necessary signal for fault detection. The method is simulated using a small
aircraft and an environment model, which includes sensor noise, wind effect, and
turbulence. In the simulation, angle of attack (AOA) sensor faults are considered.
In [5], a new method is developed for a nonlinear flight model to detect input
and output sensor faults. In the study, a polynomial residual generator structure is
used. This approach, which is robust against faults, can filter the error signal and
decrease the effect of the additive noise and error that occurs due to linearization and
measurement. The input errors in sensors are separated as known input, disturbance,
and faults vectors. The designed system describes each input–output relation by
polynomials and uses residuals to detect and isolate faults. Different fault detection
methods are also used for the same fault conditions to better evaluate the proposed
method. For this purpose, nonlinear geometric approach (NLGA), unknown input
Kalman observers (UIKF), and neural networks (NN) are investigated.
In [6], the NN method is used for fault detection and isolation (FDI), and is
applied to a nonlinear UAV model. In the study, a radial basis function (RBF) type
neural NN is designed for modeling purposes. A novel residual generation technique
is also used to overcome wrong fault error situations. A nonlinear UAV model is
considered and longitudinal equations are used to simulate the fault situation in the
pitch angle gyro.
In [7], a similar system to that used in [4] is investigated. The characteristics of
the sensor signal are again used for FDI. The design is the same as that described in
[4], except for the dynamic dependency filter structure, which changed from linear
autoregressive to pooled autoregressive (PAR). The idea comes from changing the
linear structure of the filter in [4] to obtain better results. The results are again
positive, as in [4].
In [8], interacting multiple model (IMM) and unscented Kalman filter (UKF)
methods are used with nonlinear aircraft dynamics for fault detection purposes.
The differences between using the INS equations and a dynamic model are shown
for the situations without faults. The IMM-UKF structure is proved to estimate the
states of the system in the presence of one faulty sensor. In conclusion, the IMM-
UKF structure works properly if used together for fault detection purposes and is an
effective new way of finding the faults in the system.
In [9], an FDI system is proposed for a nonlinear model. There are widely
used FDI algorithms for linear systems; however, these methods are not sufficient
7.2 Kalman Filter-Based Integrated ADS/GPS Navigation System 111
for nonlinear systems and can cause errors in fault detection. For this reason, an
input estimation structure that estimates the input values using inputs–outputs of the
system and IMU-INS measurement values is proposed in that study.
An observer/Kalman filter identification method for fault diagnosis in a small
autonomous helicopter is proposed in [10]. The helicopter model under nominal
conditions is found using input–output test values in an observer/Kalman filter.
The differences between the designed model and the system outputs let the users
determine the faults in the system. The model of the system, which is the helicopter
model in this case, can be determined by Markov parameters and a Hankel matrix
which has the same rank as the system.
In the given studies, different Kalman filter techniques (EKF, UKF, IMM-UKF,
AMEKF), model-based methods such as using sensor signal characteristics, residual
generators using polynomials, and NN, are used for state estimation and FDI. The
methods vary due to the difficulty of the problem (multiple sensor fault detection,
redundancy of sensors, etc.) where the FDI algorithm is used (aircraft, helicopter,
etc.) and the seriousness of the fault. From these studies, we can understand that the
Kalman filter technique can detect and isolate faults using different dynamic models
and, in many studies, it is also used as a comparison test. In conclusion, this filtering
technique is a useful method in determining faults in sensors.
In the literature survey, it is seen that many of the studies on the subject include
fusing INS and IMU measurements with Global Positioning System (GPS) data.
However, there are few studies on obtaining the best data for flight control from the
air data system (ADS). In this study, the problem of obtaining the most appropriate
data from ADS measurements is considered and ADS data FDI, determining the best
data from the ADS, and estimation of the state subjects are investigated using high-
accuracy GPS/INS measurements and the Kalman filtering technique as a basis.
Here, Veadsx , Veadsy , Veadsz are the speed errors of the ADS, each in Cartesian
coordinates. In the ADS, the true airspeed error is mainly the result of the wind
112 7 Fault Detection, Isolation and Data Fusion for UAV Air Data System
Vadsx
Vadsy
Vadsz V̂adsx
Air Data V̂adsy
System +
− V̂adsz
Vgpsx V̂eadsx
Vgpsy
V̂eadsy
Vgpsz + V̂eadsz
Global Optimal
Positioning Kalman
System
− Filter
Fig. 7.1 Kalman filter-based air data system (ADS)/Global Positioning System (GPS) fusion
scheme
speed. Therefore, choosing the ADS errors as the vector parameters shown in Eq.
(7.1) is for the purpose of obtaining the true airspeed error, which is relatively high,
and, thus, the wind speed.
The uncorrelated model, which is used in the exponential correlation function
for the stagnant processes, will be used [11, 12]. This expression is in discrete and
matrix form. This model is very suitable for simulation purposes [11]:
2 3 2 32 3
Veadsx .k C 1/ 1 ˇVadsx T 0 0 Veadsx .k/
4 Veadsy .k C 1/ 5 D 4 0 1 ˇVadsy T 0 5 4 Veadsy .k/ 5
Veadsz .k C 1/ 0 0 1 ˇVadsz T Veadsz .k/
2 1 2 2
3
ˇVadsx T 2 ˇ Vadsx T wVadsx
6 7
C 4 ˇVadsy T 12 ˇ 2 Vadsy T 2 wVadsy 5 (7.2)
ˇVadsz T 12 ˇ 2 Vadsz T 2 wVadsz
where ˇVadsx ; ˇVadsy ; ˇVadsz are the variables which are inverses of the true airspeed
correlation time, wVadsx ; wVadsy ; wVadsz are the Gauss distributed noises of the ADS
speed measurements, and T is the sampling period.
Using the ADS and GPS speed measurement differences as the measurements
observation vector in the Kalman Filter, the observation vector can be written as:
In Eq. (7.3), Veadsx, Veadsy, Veadsz are the true airspeed measurement errors of
the ADS and, at the same time, are the wind speeds, and vVadsx ; vVadsy ; vVadsz and
vVgpsx ; vVgpsy ; vVgpsz are the zero-mean Gauss noises of the ADS and GPS measure-
ments, respectively. Namely, the difference between the GPS speed measurement
and the ADS true airspeed measurement gives us the ADS speed error and the
wind speed during flight. However, this wind speed information includes the random
noises of both systems. The measurement statements of Eq. (7.3) can be written in
matrix form as:
2 3 2 3 2 3 2 3
z1 .k/ VADSX VGPSX 100 vVx
z.k/ D 4 z2 .k/ 5 D 4 VADSY VGPSY 5 D 4 0 1 0 5x.k/ C 4 vVy 5 (7.4)
z3 .k/ VADSZ VGPSZ 001 vVz
„ ƒ‚ …
H.k/
For integration of the ADS and INS/GPS data, the Kalman filtering equations
below are used:
b
x .k=k/ D b
x .k=k 1/ C K.k/.k/
.k/ D z.k/ H.k/b x .k=k 1/
K.k/ D P .k=k/ H T .k/R1 .k/
(7.5)
P .k=k/ D .I K.k/H.k// P .k=k 1/
P .k=k 1/ D ˆ .k; k 1/ P .k 1=k 1/ ˆ .k; k 1/
CG .k; k 1/ Q .k 1/ G .k; k 1/
where b x .k=k/ is the estimation value, I is the identity matrix, K(k) is the gain value
of the optimal filter, P(k/k) is the estimation error correlation matrix, P .k=k 1/
is the extrapolation error correlation matrix, H(k) is the measurement matrix, (k)
is the innovation process, Q .k 1/ is the system measurement correlation matrix,
G .k; k 1/ is the noise transfer matrix, and ˆ .k; k 1/ is the system transfer
matrix:
2 3
1 ˇVadsx T 0 0
ˆ .k; k 1/ D 4 0 1 ˇVadsy T 0 5 (7.6)
0 0 1 ˇVadsz T
b
x .0=0/ D x.0/
P .0=0/ D P.0/
To start the Kalman filter process, the x(0) and P(0) initial conditions, system
error correlation matrix Q(k), and measurement error correlation matrix R(k) must
first be determined. To perform the filtering, we need system and measurement error
models.
114 7 Fault Detection, Isolation and Data Fusion for UAV Air Data System
In Fig. 7.1, as we can see in the given integration scheme, the Kalman filter
outputs are the ADS speed error estimation values, which are actually the wind
speeds b
Veadsx ; b
Veadsy ; b
Veadsz . These estimated values are subtracted from the ADS
measurements to achieve the speeds b Vadsx ; b
Vadsy ; b
Vadsz . We can use the following
formulas for determining the speed values:
b
Vadsx D Vadsx b
Veadsx
b b
Vadsy D Vadsy Veadsy (7.7)
Vadsz D Vadsz b
b Veadsz
A standard central Kalman filter with central measurements and a covariance matrix
can cause computational overload. In addition, any sensor fault can degrade the
overall performance of the filter. To overcome these disadvantages, a noncentral
Kalman filter is useful. A noncentral Kalman filter is a two-stage data-processing
technique and is composed of one basic filter and one or more local filters [13–15].
Firstly, the local filters process the data optimally to obtain the best result using each
set. Then, the central filter fuses local filter results and produces the most optimal
value. The described method is an effective technique for integrating multiple sensor
measurement sets. The federated Kalman filter has many advantages over the other
noncentral filtering techniques [11, 13–15].
The federated Kalman filter-based navigation sensor fusion scheme, which fuses
INS, GPS, and other navigation sensors, is given in Fig. 7.2.
The aforementioned federated Kalman filtering is a two-level data-processing
technique. Local filters that process data from the sensors are fused in a federated
filter. All the local filters estimate the same variables; thus, a correlation is present
between them. A global optimal estimation and error variance can be described by
the following equations:
1 1
Pm D P1 1
1 C P2 C C PN (7.8)
Xem D Pm P1 1 1
1 X e1 C P 2 X e2 C C P N X eN
In the formulas, Pi is the covariance matrix of the ith filter and Xei is the
estimation value of the ith filter.
7.3 Federated Kalman Filter-Based Integrated ADS and GPS/INS Data 115
In this study, INS/GPS measurements, which have high accuracy, and ADS
measurements, which have low accuracy but high frequency, are integrated using
the Kalman filtering technique in order to obtain high-accuracy measurement data
at high frequency. The designed system based on the indirect Kalman filtering
technique is shown to be successful in finding the wind speed data using the known
error or UAV dynamics of the system and determined statistical values.
In applications, generally, ADS and GPS or GPS/INS are integrated on the basis
of a Kalman filter. A federated Kalman filter-based navigation sensor fusion scheme
with FDI block to fuse the reliable data coming from the three ADS measurement
sets and integrated GPS/INS (faulty data are determined and isolated) is shown in
Fig. 7.3.
In the figure, it can be seen that the federated filter is composed of local filters
that process ADS navigation measurements and GPS/INS measurements (position
and velocity), together with the ADS data, FDI block for faulty data determination
and isolation, and a main filter that fuses the estimates of the local filters after the
elimination of faulty data. The main best variable estimates and error variance of
the filters can be found using the equations given below [7, 9]:
1 1
Pm D P1 1
1 C P2 C P3 (7.9)
Xem D Pm P1 1 1
1 X e1 C P 2 X e2 C P 3 X e3
116 7 Fault Detection, Isolation and Data Fusion for UAV Air Data System
Fig. 7.3 Kalman filter-based ADS and GPS/inertial navigation system (INS) data fusion with fault
detection and isolation (FDI)
The aforementioned algorithm can be used to derive the best data from the ADS
measurements.
X
k
ˇ.k/ D Q T .j/.j/
Q (7.10)
jDkMC1
7.4 Sensor FDI Algorithms 117
Q
where .j/ is the normalized innovation sequence of the Kalman filter and M is the
width of the sliding window.
The Kalman filter normalized innovation can be calculated as follows:
1=2
Q
.k/ D H.k/P .k=k 1/ H T .k/ C R.k/ .k/ (7.11)
then the threshold value,
2˛,Ms , can be found. Hence, when the hypothesis H1 is true,
the statistical value of ˇ(k) will be greater than the threshold value
2˛,Ms , i.e.:
H0 W ˇ.k/
2˛;Ms ; 8k
(7.13)
H1 W ˇ.k/ >
2˛;Ms ; 9k
If the fault is detected, then it is necessary to determine what sensor set is faulty. For
this purpose, the s-dimensional sequence Q is transformed into s one-dimensional
sequences to isolate the faulty sensor, and for each one-dimensional sequence
Q i .i D 1; 2; : : : ; s/, the corresponding monitoring algorithm is run. The statistics of
the faulty sensor set are assumed to be affected much more than those of the other
sensors. Let the statistics be denoted as i (k). When max fi .k/=i D 1; 2; : : : ; sg D
m .k/ for i ¤ j and i .k/ ¤ j .k/, it is judged that the sensor set has failed.
.
Let the statistics that is a rate of sample and theoretical variances, b
i2
i2 , be
used to verify the variances of one-dimensional innovation sequences i .k/; i D Q
1; 2; : : : ; s. When Q i N .0;
i /, it is known that [18]:
i
2M1 ; 8i; i D 1; 2; : : : ; s (7.14)
i2
where i D .M 1/b
i2 .
2
As
i D 1 for the normalized innovation sequence, it follows that:
Using Eq. (7.15), it can be proved that any change in the mean of the normalized
innovation sequence can be detected [19]. When a fault affecting the mean or
variance of the innovation sequence occurs in the system, the statistics i exceeds
118 7 Fault Detection, Isolation and Data Fusion for UAV Air Data System
the threshold value
2˛;M1 , depending on the level of significance ˛ and degrees of
freedom .M 1/.
Three ADS sets are taken into consideration. The sets are composed of values taken
from an AOA sensor, sideslip angle sensor (AOS), and speed values determined
from true airspeed sensor (TAS) measurements. There are three AOA, two AOS, and
three TAS sensors. The first set is built using AOA sensor 1 (AOA1), AOS sensor 1
(AOS1), and TAS sensor 1 (TAS1). The second set includes AOA sensor 2 (AOA2),
AOS sensor 2 (AOS2), and TAS sensor 2 (TAS2). The third set is composed of AOA
sensor 3 (AOA3), AOS2 (as there are two different AOS sensors), and TAS sensor
3 (TAS3).
Three different cases are investigated in the simulations [20]:
1. TAS1 bias D 10 knots (C10 knots disturbance on the first set’s true airspeed data
after 80 s)
2. AOA3 bias D 5ı (constant C5ı bias on the third set’s angle of attack data after
80 s)
3. AOS1 bias D5ı (constant C5ı bias on the first set’s sideslip angle data after 80 s)
The faulty data are eliminated using the statistical test as described earlier. After
that, data from another set that are known to be true are used in place of the faulty
data. For example, when the TAS1 data is determined as faulty, TAS2 will be used
instead to make set 1 working. Of course, as there are only two sideslip angle
sensors, when one AOS value is wrong, there is only one other value to replace it.
In the second simulation, the AOA3 value is faulty and it is replaced by the AOA1
value as soon as the system detects the fault.
These sets are simulated with systems that isolate the faults and eliminate the
faulty data (Figs. 7.4, 7.5, 7.6, and 7.7) and the systems that do not perform isolation
(Figs. 7.8, 7.9, 7.10, and 7.11). In the figures, “beta” denotes the sideslip angle and
“KTAS” refers to the knots true airspeed.
The case for the system with FDI when the TAS1 bias equal to 10 knots occurs at the
80th second is investigated. Using the statistical tests, faulty data can be determined
and isolated. The variables AOA and TAS estimated using the filters are given in
Figs. 7.4 and 7.5, respectively.
7.5 Simulation Results for Indirect Kalman Filter-Based ADS and GPS/INS. . . 119
Fig. 7.4 Angle of attack (AOA) estimation values with fault isolation [when true airspeed sensor
1 (TAS1) bias D 10 knots]
Fig. 7.5 TAS estimation values with fault isolation (when TAS1 bias D 10 knots)
120 7 Fault Detection, Isolation and Data Fusion for UAV Air Data System
Fig. 7.6 Wind velocity estimation values in the x-direction with fault isolation (when TAS1
bias D 10 knots)
In Fig. 7.6, the ADS speed error estimation in the x-direction (wind velocity
estimation in the x-direction) and real wind velocity in the x-direction are shown
together.
It can be seen that the wind velocity in the x-direction is accurately estimated by
the federated Kalman filter, even though the first and second sets do not give good
results at different time intervals. We can also tell that some of the measurement
values are eliminated by the system that includes the federated filter technique.
Secondly, the case where the bias in the third set’s angle of attack value is C5ı
after the 80th second is considered. In these simulations, it must be remembered that
the isolation algorithm is also working and that the faulty data are eliminated.
Again, for this set, fault isolation is carried out and AOA and TAS estimation
values are shown in Figs. 7.7 and 7.8, respectively.
As seen in Figs. 7.7 and 7.8, the estimated values of AOA and TAS via the
federated filter are close to the real values.
Thirdly, the case where the bias on the first set’s sideslip angle value is C5ı after
the 80th second is considered. In these simulations, it must be remembered that the
isolation algorithm is also working and that the faulty data are eliminated. Wind
speed values in the y-direction for this case are presented in Fig. 7.9.
As seen from Fig. 7.9, the wind velocity in the y-direction is accurately estimated
by the federated Kalman filter. Similar simulation results are obtained for the rest of
the parameters of the UAV ADS.
7.5 Simulation Results for Indirect Kalman Filter-Based ADS and GPS/INS. . . 121
Fig. 7.7 AOA estimation values with fault isolation (when AOA3 bias D C5ı )
Fig. 7.8 TAS estimation values with fault isolation (when AOA3 bias D C5ı )
122 7 Fault Detection, Isolation and Data Fusion for UAV Air Data System
Fig. 7.9 Wind speed values in the y-direction with fault isolation (when AOS1 bias D C5ı )
The federated filter is investigated for the same three cases. If isolation is not made,
the federated filter can give faulty results [20]. Again, firstly, the results when the
TAS1 bias is equal to 10 knots is given. Graphs that include the AOA, TAS, and
x-direction wind speed estimation values are given in Figs. 7.10, 7.11, and 7.12,
respectively.
As can be seen from Figs. 7.11 and 7.12, the TAS and wind speed estimation
values in the x-direction of the federated Kalman filter diverge. The cause is the
data coming from the first estimation set. The filter takes the values from the first set
into account; thus, the federated values are not very good for this reason (no faulty
data isolation is made).The federated filter estimated values of AOA are not affected
significantly by the TAS1 fault, as shown in Fig. 7.10.
Next, the case where the bias in AOA in set 3 is C5ı is considered. The estimated
AOA and TAS values for this case are given in Figs. 7.13 and 7.14, respectively.
As seen in Fig. 7.13, the estimated AOA values of the federated Kalman filter
also diverge. The values with errors, which are not isolated and used in calculations,
cause this divergence. The federated filter estimated values of TAS are not affected
significantly by the AOA3 fault, as shown in Fig. 7.14.
Wind speed estimation results in the y-direction when the AOS1 bias D 5ı and
fault isolation is absent are shown in Fig. 7.15.
7.5 Simulation Results for Indirect Kalman Filter-Based ADS and GPS/INS. . . 123
Fig. 7.10 AOA estimation values without fault isolation (when TAS1 bias D 10)
Fig. 7.11 TAS estimation values without fault isolation (when TAS1 bias D 10 knots)
124 7 Fault Detection, Isolation and Data Fusion for UAV Air Data System
Fig. 7.12 Wind velocity estimation values in the x-direction without fault isolation (when TAS1
bias D 10)
Fig. 7.13 AOA estimation values without fault isolation (AOA3 bias D C5ı )
7.5 Simulation Results for Indirect Kalman Filter-Based ADS and GPS/INS. . . 125
Fig. 7.14 TAS estimation values without fault isolation (AOA3 bias D C5ı )
Fig. 7.15 Wind speed estimation results in the y-direction without fault isolation (when AOS1
bias D 5ı )
126 7 Fault Detection, Isolation and Data Fusion for UAV Air Data System
As can be seen from the presented graphs, wind speed estimation values in the
y-direction of the federated Kalman filter diverge.
In UAVs, having highly accurate ADS parameters is needed for the purposes of
control. To obtain the required speed data, fusing the measurements of the system
is a must. In this study, GPS/INS measurements, which have high accuracy, and
ADS measurements, which have low accuracy but high frequency, are integrated
using the Kalman filtering technique in order to obtain high-accuracy measurement
data at high frequency. The designed system based on the indirect Kalman filtering
technique is shown to be successful in finding the wind speed data using the known
error or UAV dynamics of the system and determined statistical values.
In this chapter, the problem of obtaining the most appropriate data from ADS
measurements is considered and ADS data FDI, determining the best data from
the ADS, and estimation of the state subjects are investigated using high-accuracy
GPS/INS measurements and the Kalman filtering technique as a basis.
FDI algorithms are developed and diagnostic tests performed for faulty cases.
The mentioned tests are based on the statistical properties of the normalized
innovation process of the Kalman filter.
To fuse the data originating from the different measurement sets, a federated
Kalman filter is presented. The mentioned algorithms can be used to derive the
best data from the ADS measurements. The best data are produced using the FDI
algorithm and leaving the faulty data out before fusing with the federated filter. As a
result, the air data parameters are estimated by the federated Kalman filter very well
and the estimated values converge to the real values.
References
1. Jayaram S (2010) A new fast converging Kalman filter for sensor fault detection and isolation.
Sensor Rev 30(3):219–224
2. Bryson M, Sukkarieh S (2004) Vehicle model aided inertial navigation for a UAV using low-
cost sensors. In: Proceedings of the 2004 Australasian conference on robotics and automation,
Canberra, Australia, December 2004
3. Mosallae M, Salahshoor K (2008) Sensor fault detection using adaptive modified extended
Kalman filter based on data fusion technique. In: Proceedings of the 4th international
conference on information and automation for sustainability (ICIAFS08), Colombo, Sri Lanka,
December 2008
4. Samara PA, Sakellariou JS, Fouskitakis GN, Fassois SD (2003) Detection of sensor abrupt
faults in aircraft control systems. In: Proceedings of the 2003 IEEE conference on control
applications (CCA 2003), Istanbul, Turkey, June 2003
5. Bonfe M, Castaldi P, Geri W, Simani S (2006) Fault detection and isolation for on-board
sensors of a general aviation aircraft. Int J Adapt Control Signal Process 20:381–408
References 127
6. Samy I, Postlethwaite I, Gu DW (2009) Sensor fault detection and accommodation using neural
networks with application to a non-linear unmanned air vehicle model. In: Proceedings of the
IMechE Part G: J Aerosp Eng, vol 224
7. Samara PA, Fouskitakis GN, Sakellariou JS, Fassois D (2008) A statistical method for the
detection of sensor abrupt faults in aircraft control systems. IEEE Trans Control Syst Technol
16(4):789–798
8. Cork L, Walker R (2007) Sensor fault detection for UAVs using a nonlinear dynamic model
and the IMM-UKF algorithm. IEEE Inf Decis Control, pp 230–235
9. Marzat J, Piet-Lahanier H, Damongeot F, Walter E (2010) Fault diagnosis for nonlinear aircraft
based on control-induced redundancy. In: Proceedings of the conference on control and fault-
tolerant systems (SysTol), Nice, France, October 2010
10. Heredia G, Ollero A (2009) Sensor fault detection in small autonomous helicopters using
observer/Kalman filter identification. In: Proceedings of the 2009 IEEE international confer-
ence on mechatronics, Malaga, Spain, April 2009
11. Hajiyev C (2010) Experimental data processing methods and engineering applications. Nobel
Publication and Distribution Inc., Ankara (in Turkish)
12. Grishin YP, Kazarinov YM (1985) Fault-tolerant dynamic systems. Radio i Svyaz, Moscow (in
Russian)
13. Carlson NA (1990) Federated square root filter for decentralized parallel processors. IEEE
Trans Aerosp Electron Syst 26:517–525
14. Carlson NA, Berarducci MP (1994) Federated Kalman filter simulation results. Navigation
41:297–321
15. Hajiyev C, Tutucu A (2001) INS/GPS integration using parallel Kalman filtering. In: Proceed-
ings of the 15th IFAC symposium on automatic control in aerospace, Bologna/Forli, Italy,
pp 453–458
16. Willsky AS (1976) A survey of design methods for failure detection in dynamic systems.
Automatica 12(6):601–611
17. Hajiyev C, Caliskan F (2003) Fault diagnosis and reconfiguration in flight control systems.
Kluwer Academic Publishers, Boston
18. Hajiyev C (2006) Innovation approach based measurement error self-correction in dynamic
systems. Measurement 39:585–593
19. Hajiyev C, Caliskan F (2005) Sensor and control surface/actuator failure detection and isolation
applied to F-16 flight dynamic. Aircr Eng Aerosp Technol 77(2):152–160
20. Hajiyev Ch, Vural SY, Kargin V (2012) Fault detection, isolation and federated filtration in
UAV air data system. In: Proceedings of the national aerospace conference, 12–14 September
2012. Turkish Air Force Academy, Istanbul, UHUK-2012-043, 16 (in Turkish)
Chapter 8
Stability Analysis for Unmanned Aerial Vehicles
8.1 Trimming
The equations of motion for unmanned aerial vehicles (UAVs) are given in Chap. 2.
In this chapter, the equations of motion are further examined and their linearization
is discussed in detail for the stability analysis. In this procedure, the first step is to
investigate the trim conditions.
xP D f .x; u/ ; (8.1)
where u represents the system control inputs and x represents the states of the
system. The trim condition for the system can be found as:
f .b u/ D 0:
x;b (8.2)
If the system is in equilibrium with the given states and control inputs, the
variables do not change. Stable (trimmed) flight condition for a UAV is an example
of such a situation. In this case, we can say that a subspace of the UAV states is in
equilibrium. In the literature, this case is called the trim condition. Wings-level,
fixed-altitude turn, and climb–descent maneuvers are different examples of trim
conditions.
The states of a UAV can be calculated for different trim conditions using a
minimization function. The linearization of the equations is necessary for controller
design.
We can find the linearized equations of motion for a UAV starting with the
equations given in Chap. 2.
Around the equilibrium point, we have:
b
xP D f .b
x;b
u/ (8.3)
P
xPQ D xP b
x; (8.5)
xPQ D f .x; u/ f .b u/ ;
x;b (8.6)
xPQ D f .x C xQ xQ ; u C uQ uQ / f .b u/ ;
x;b (8.7)
xPQ D f .b u C uQ / f .b
x C xQ ;b u/ :
x;b (8.8)
@f .b
x;b
u/ @f .b
x;b
u/
xPQ Š xQ C uQ : (8.10)
@x @u
@f @f
In order to linearize the equations, we need to solve for @x and @u at the
equilibrium point.
The longitudinal and lateral motions can be handled separately. The variables that
we choose to describe the longitudinal and lateral motions of the UAV are h, V, ˛,
u, w, q, and ˇ, v, p, r, ·, , respectively.
8.2 Derivation of the Transfer Functions 131
The longitudinal equations of motion for the UAV can be found by taking
the lateral variables as zero and using definitions of ˛ D arctan wu and Vt D
p
u2 C w2 , where Vt is the true airspeed.
The differentiation matrix for the states can be found as:
2
3
@Pu @Pu @Pu @Pu @Pu
; ; ; ;
6 @u 7
6
@w
@q
@
@h 7
6 7
6 @wP @wP @wP @wP @wP 7
6 ; ; ; ; 7
6 @u 7
6
@w
@q
@
@h 7
6 7
@f 6 @Pq @Pq @Pq @Pq @Pq 7
6 ; ; ; ; 7
D 6 @u @w @q @ @h 7 (8.11)
@x 6 ! ! ! ! !7
6 7
6 @P @P @P @P @P 7
6 ; ; ; ; 7
6 7
6 @u @w @q @ @h 7
6 ! ! ! ! ! 7
6 7
4 @hP @hP @hP @hP @hP 5
; ; ; ;
@u @w @q @ @h
Similarly, the differentiation matrix for the control inputs can be derived and the
longitudinal state and control matrices—Alon and Blon , respectively—can be found.
The equations for the longitudinal stability coefficients are given below:
V 2 cxıe S
Xıe D (8.15)
2m
Sprop cxp ıa h
Xıt D (8.16)
m
uS Scz˛ w Scczq uq
Zu D q C Œczo C cz˛ ˛ C czıe ıe C (8.17)
m 2m mV
wS Scz˛ u Scczq uq
Zw D Œczo C cz˛ ˛ C czıe ıe C (8.18)
m 2m mV
Sczq Vc
Zq D u C (8.19)
2m
132 8 Stability Analysis for UAV
Sczıe V 2
Zq D (8.20)
2m
Scmıe V 2 c
Mıe D (8.24)
2Jy
The derivation procedure for the lateral equations of motions and, therefore,
the lateral state and control matrices—Alat and Blat , respectively—are similar. The
equations for the lateral stability coefficients are given below:
SVair 2
Yıq D cyıa (8.28)
2m
SVair 2
Yır D cyır (8.29)
2m
vb2 S Sbv
Lv D cpp p C cpr r C cpo C cpˇ ˇ C cpıa ıa C cpır ır
8mVair 2m
Sbcpˇ p 2
C u C w2 (8.30)
4m
SV 2 air b
Lp D 1 q C cpp (8.31)
8m
SV 2 air b
Lr D 2 q C cpr (8.32)
8m
8.2 Derivation of the Transfer Functions 133
SbVair 2
Lıq D cpıa (8.33)
4m
SbVair 2
Lır D cpır (8.34)
4m
vb2 S Sbv
Nv D crp p C crr r C cro C crˇ ˇ C crıa ıa C crır ır
8mVair 2m
Sbcrˇ p 2
C u C w2 (8.35)
4m
SV 2 air b
Np D 3 q C crp (8.36)
8m
SV 2 air b
Nr D 4 q C crr (8.37)
8m
SbVair 2
Nıq D crıa (8.38)
4m
SbVair 2
Nır D crır (8.39)
4m
The next step is to give the linearized state equations for the UAV. The equations
in the state space form are:
0 1 2 30 1 2 3
uP .Xu/ .Xw/ .Xq/ .g cos .// 0 u Xıe Xıt
B C 6 7B C 6 7
B wP C 6 .Zu/ .Zw/ .Zq/ .g sin .// 0 7 B C 6 07 " #
B C 6 7 B w C 6 Zıe 7 ı
B C 6 7B C 6 7 e
B qP C D 6 Mu Mw Mq 0 0 7 B q C C 6 Mıe 07
B C 6 7B C 6 7 ıt
B P C 6 0 0 1 0 7 B C
05@ A 4 0 6 0 7
@ A 4 5
hP sin ./ cos ./ 0 t 0 h 0 0
(8.40)
Here, ten states are used to describe the lateral and longitudinal motions of the UAV.
The heading angle equation is also added to the system.
Finally, the linearized equations of motion for the Zagi UAV are found as below.
These equations will be used when we design the various controller algorithms that
are presented in the next chapters of this book. The values given in the matrices Alon ,
Alat , Blon , and Blat are for the trim conditions. The longitudinal equations are:
xP D Ax C Bu (8.43)
2 3
0:3356 1:3181 1:9276 9:6610 0
6 7
6 1:7916 3:9003 9:8215 1:7035 07
6 7
6 7
Alon D 6 0:7020 3:5375 11:3920 0 07 (8.44)
6 7
6 0 0 1:0000 0 07
4 5
0:1736 0:9848 0 17:4865 0
2 3
0:7436 6:8728
6 7
6 3:7855 0 7
6 7
6 7
Blon D 6 47:9170 0 7 (8.45)
6 7
6 0 0 7
4 5
0 0
ıe
ulon D (8.46)
ıt
T
Xlon D u w q h (8.47)
ıa
ulat D (8.50)
ır
T
xlat D v p r ' ‰ (8.51)
The effect of elevator and thrust inputs on the states can be represented with the
transfer functions that can be derived using the Laplace transformation as:
Y.s/
D CŒsI A1 B C D (8.52)
U.s/
Y.s/ Cadj .sI A/ B jsI Aj
D (8.53)
U.s/ jsI Aj
where Y(s) is the output and U(s) is the control input. The determinant jsI Aj
gives the characteristic equation. MATLAB is used to calculate all the transfer
functions. As an example, if the output variable is chosen as u, the forward velocity,
then we can calculate ıeu and ıtu . First, the “state_space_long D ss(Alon, Blon,
[1,0,0,0,0], 0)” command gives us the state space model for the system and then
the “tf(state_space_long)” command can be used to derive the equations for the
forward velocity [1–3]:
In a similar way, the output matrix C (given as [1,0,0,0,0] in the “ss” function)
can be changed to obtain the desired transfer functions and the effect of control
inputs on the rest of the longitudinal states can be found:
q 4:825s2 C 68:38s
D 4 ; (8.59)
ıt s C 15:63s C 88:02s2 C 62:64s C 87:23
3
4:825s C 68:38
D 4 ; (8.61)
ıt s C 15:63s C 88:02s2 C 62:64s C 87:23
3
If the same process is repeated for the lateral equations of motion, we can obtain the
transfer functions for the lateral motion as:
In the next sections, we investigate the characteristic equations for the stability
analysis.
As we can see, there are two equations that describe the longitudinal motion
of the UAV. These two equations represent two different modes: the phugoid mode,
which is easy to control but takes a long time to dampen, and the short-period mode,
which makes the plane show fluctuation motions with high frequency as the name
refers and has a high damping ratio.
In the short-period mode, the changes in the attack angle and pitch angle are
significant, whereas the change in the forward velocity is negligible. On the other
hand, in the phugoid mode, the changes in u and pitch angle are significant, and
this mode can be described as the change in the kinetic and potential energies of the
aircraft. Next, we investigate the roots of the characteristic equation [2–5].
The roots of the equation are 7.5215 C 4.6367i, 7.5215 4.6367i,
0.2935 C 1.0155i, and 0.2935 1.0155i. There are two roots which are close to
the imaginary axis and describe the phugoid motion and two roots which are far
away from the imaginary axis and describe the short-period motion of the UAV. The
angular frequency (w), damping factor (—), and period (T) values are calculated and
given in Table 8.1. The UAV is longitudinally stable in the given condition because
of the location of the roots of the characteristic equation.
As can be seen from the calculations and the values given in the table, the roots
are on the left of the imaginary axis for the longitudinal motion of the UAV and the
damping factor and natural frequency values are in the desired zone.
The characteristic equation of the lateral motion is investigated and the stability
condition is checked. The characteristic equation is as follows:
2
s C 2:5638s C 20:9734 .s 0:0435/ .s C 2:1658/ D 0 (8.75)
As can be seen, the lateral characteristic equation is also of the fourth order. It
represents three different modes, which are the roll, Dutch roll, and spiral modes [2,
4, 6].
The roll mode is mainly a pure rolling motion. If the aircraft is not stable, this
motion will lead to an increase in rolling in the direction in which the aircraft
turns. The roll mode is represented by one root and gives a first-order nonoscillatory
response; it is generally stable at low attack angles.
The spiral mode also gives a first-order response and it is represented by one first-
order root. This mode includes both rolling and yawing motions. If it is not stable,
the rolling motion increases while yawing, and the aircraft starts a spiral motion.
However, instability in the spiral mode is tolerable under some specific conditions.
The Dutch roll mode is triggered by perturbations in the sideslip angle and
includes yaw, roll, and sideslip changes combined. If the directional stability is not
high enough to lead the aircraft into the spiral mode, the Dutch roll mode occurs
and the aircraft shows oscillatory yaw and roll motions.
The roots that symbolize each motion are found from the characteristic equation
and given in Table 8.2.
The results show that the spiral mode root is positive and this mode is unstable.
The roll mode is stable because of the root at 2.1658. The Dutch roll mode
Table 8.2 Roll mode, Dutch roll mode and spiral mode characteristic values
Natural
Motion frequency Damping factor Roots Time constant
Roll N/A N/A 2.1658 0.4616 s
mode
Spiral N/A N/A 0.0435 (unstable) 16 s (for double amplitude)
mode
Dutch 4.58 rad/s 0.27 1.282 C 4.397i 0.808 s
roll mode 1.282 4.397i
References 139
damping factor is 0.27 and its natural frequency is 4.58 rad/s. For stable modes,
the times to get back to 36.8 % of the amplitude are given. For the spiral mode, the
time to double amplitude is calculated and given in the table.
8.5 Conclusion
The linearized equations of motion of UAVs are used to find the transfer functions
that relate the actuator inputs to the state outputs. The UAV is analyzed for
lateral and longitudinal stability by means of these transfer functions. First, the
characteristic equations are obtained using the transfer functions and the stability
condition is investigated regarding the roots of the polynomials. The short-period
and phugoid modes for longitudinal motion and the roll, Dutch roll, and spiral
modes for lateral motion are discussed. It is shown that the UAV is stable in the
analyzed flight condition.
References
1. Tewari A (2002) Modern control design with Matlab and Simulink. Wiley, Chichester
2. Pamadi BN (2003) Performance, stability, dynamics, and control of airplanes. AIAA education
series, Reston
3. Etkin B, Reid LD (1996) Dynamics of flight: stability and control. Wiley, New York
4. Yechout TR, Morris SL, Bossert DE, Hallgren WF (2003) Introduction to aircraft flight
mechanics. AIAA education series, Reston
5. Blakelock JH (1991) Automatic control of aircraft and missiles. Wiley, New York
6. Nelson RC (1998) Flight stability and automatic control. McGraw-Hill, New York
Chapter 9
Classic Controller Design for Unmanned Aerial
Vehicles
where c(t) is the output, e(t) is the error for making the actual state reach the desired
condition as e.t/ D xd .t/ xa .t/, and K is the proportional gain. Here, xd (t) is the
desired state condition and xa (t) is the actual (or observed) state.
Integral gain takes the integral of the error input and can be described by:
Zt
1
c.t/ D e.t/dt: (9.2)
i
0
Using the Laplace transformation, the transfer function can be represented as:
1
c.s/ D e.s/: (9.3)
i s
de.t/
c.t/ D d (9.4)
dt
PI control fuses two control effects and makes the steady-state error zero. The
integral action adds a root at zero to the system transfer function and removes
the error at the steady state. On the other hand, the integral action also affects the
overall system response, so P-type control is used to obtain the required response,
compensating for the unwanted effects of the integral gain.
Another possible combination is to use the derivative effect (D) together with the
proportional gain. Then, the controller is called a PD controller and its response can
be given by:
de.t/
c.t/ D Ke.t/ C Kd (9.8)
dt
where the transfer function is:
C.s/
D K .1 C d s/ : (9.9)
E.s/
In case we use the PD-type controller, the controlled system acts quickly because
of the differential effect; however, the steady-state error remains. If we just want
systems to reach the desired state in a short time, then this type of controller may be
used in some cases.
The PID-type controller has all three control type effects and it gives the response
according to the control law of:
Zt
K de.t/
C.t/ D Ke.t/ C e.t/dt C Kd : (9.10)
i dt
0
The PID-type controller makes the system respond quickly and removes the
steady-state error. By choosing the appropriate gains, we can obtain the required
response [1].
In our design, PID effects are used together to control the pitch angle, speed,
heading, and altitude. A comparison is made to choose the best controller for the
system. The controller coefficients are determined using the root locus method (by
choosing the gain), carrying out simulations, and using optimization methods such
as determining the signal response specifications through MATLAB.
144 9 Classic Controller Design for the UAV
u
w
deltae1
q
theta
h
deltat
Vt
all
Longitudinal Dynamics2
Qref de Theta
ue 10 -K-
-K-
s+10
Krg1 Scope2-de
Kamp Elevator Actuator
s
s+1
Washout filter
1
u
2 3
x' = Ax+Bu w
1 q
y = Cx+Du
deltae1
State-Space1 4
theta
2 5
deltat 6 h
all
If we use the transfer function that relates ı e to q and include the elevator actuator
function, we can choose the appropriate Kq gain for the inner loop using the root
locus analysis. In the outer loop, a washout filter is added to the system to exclude
the unnecessary signal inputs. For controlling , a P-type controller is used. K can
also be chosen using the root locus analysis.
The feedback for the pitch angular rate and pitch angle is usually done together
to obtain the damping factor and angular frequency of the required values in the
inner loop for the short-period mode and in the outer loop for the phugoid mode.
However, as we can see in the linearized longitudinal equations, the values for the
short-period mode are already in the required zone and feedback is not needed in
the inner loop. We can check this from Military Specification MIL-F-8785C or
normal aircraft standard FAR. Also, there are studies on the subject claiming that,
for small UAVs, different standards should be used to evaluate the characteristics of
the system [2–4]:
Here, the elevator actuator function, which has the time constant of 0.1 s, is modeled
as:
10
TFelevator D (9.13)
s C 10
The root locus of the system for negative values of gain is shown in Fig. 9.3.
With the K feedback, we can observe from the root locus that, for gain values of
approximately 0.5, the short-period mode damping values do not decrease much
(0.665), while the phugoid characteristics are improved (to a damping factor of
0.7). In other words, feedback improves the phugoid characteristics. The transfer
function of the system, which includes the washout filter, and its step response is as
follows (Fig. 9.4):
146 9 Classic Controller Design for the UAV
L
D (9.14)
ref Lref
L D 2; 492s8 C 7; 619s7 C 94; 320s6 C 584; 500s5 C 18; 160; 00s4 C 3; 062; 000s3
C3; 225; 000s2 C 2; 244; 000s3 C 762; 500
(9.15)
Lref D 103 0:001s11 C 0:05226s10 C 1:197s9 C 15:8s8 C 131:2s7 C 695:4s6
C 2; 279s5 C 4; 372s4 C 5; 669s3 C 4881s2 C 2; 768s C 760:9
(9.16)
As we can see, the system is stable and the steady-state error is very low.
For the feedback loop, an equivalent feedback system transfer function can be
given as:
G.s/
TFfeedback D (9.17)
1 C G.s/ .H.s/ 1/
where G(s) is the multiplication of the elevator actuator transfer function (/ı e ) and
K , and H(s) is the washout filter transfer function [5].
The position constant where the zn are the zeros and the pn are the poles can be
given as:
k zn
Kp D (9.18)
pn
The steady-state error can be calculated using the position constant as:
1
e1 D (9.19)
1 C Kp
Using the steady-state error and position constant formulas, we can choose the
gain for the feedback loop to reach a low steady-state error value. In our case, the
steady-state error that we obtain with the root locus-based design is 0.003, which is
sufficiently small. Besides, we will build an outer loop controller for the altitude in
the next section. Thus, we do not change the inner loop controller structure and we
also do not use the method given by Eqs. (9.17), (9.18), and (9.19).
angle to the altitude is found. The transfer function from the output pitch angle to the
reference input is calculated and then combined with h/, which can be calculated
by multiplying h/ı e and /ı e . As a result, we find the transfer function of h/ ref .
Then, a root locus study is made to investigate the system:
The critical gain value for stability is 0.299, as we can see in Fig. 9.5. A PID
controller can be used for such a system; however, this system is already stable
and it has a root at zero, which gives the integrator effect. Thus, instead of a PID
controller, a PD controller is designed. The PD controller gives a quick response
and is appropriate for the requirements.
Fig. 9.5 Root locus used in the design of the altitude controller (outer loop)
9.2 Classical Controller for the Longitudinal Motion 149
Firstly, the coefficient for the PID controller, which was planned to be used, was
calculated using the Ziegler–Nichols technique, such that the integral of the error
is minimum [6]. The critical values of the gain and angular frequency are Khpu D
0:299 and !u D 2:35rad=s. According to the Ziegler–Nichols rules, the P gain
must be 0.6 times the critical gain, and the integral and differential gains should be
calculated accordingly [6, 7]:
2
Tu D D 2:72s (9.23)
!u
0:6 Khpu
Khi D D 0:134 (9.26)
0:5 Thu
The PID values, calculated using the Ziegler–Nichols method, are as given above.
The gains are used in the system together with the saturators, which will keep the
reference pitch angle at the required levels. The results are satisfactory. Both the
pitch angle and the altitude reach the desired values in a short period of time.
In order to better control the speed and suppress the overshoot of the response,
the gains are recalculated. The responses of the controlled and uncontrolled systems
are given in Figs. 9.6, 9.7, and 9.8. Both the results found using the Ziegler–Nichols
gains values and tryouts are given. The improvement in the response of the system
can be seen when the gains are found by tryouts. For the tryouts, this response is
obtained by taking Khp D 0:58, Khi D 0:25, and Khi D 0:2.
The response of the altitude system is improved. However, the integral effect is
already contained in the system equations, and we prefer to use a PD-type controller,
not a PID-type controller. The gains of the PD-type controller are determined using
the Ziegler–Nichols values as a starting point. The proportional gain value Khp is
chosen as 0.1794, and tryouts are done using the previously determined Khd . The
final value of Khd that is used in the PD controller is 0.132. In this case, a zero at
1.33 is added to the system and a nonvibratory quick response is obtained. The
step response of the system with the PD controller is given in Fig. 9.9.
The response of the designed altitude controller for the case where the reference
input is 100 m and the speed is held constant at 20 m/s using the speed controller
is given in Fig. 9.10. Moreover, the response of the speed controller, which tries to
hold the speed at 20 m/s, is given in Fig. 9.11. As we can see, the PD controller is
working efficiently as an altitude controller for our system. Nonetheless, in cases
where other disturbances are present, a PID-type controller may be used instead.
In Fig. 9.12, the Simulink diagram for the final structure of the altitude controller
is given. As can be seen, the control is obtained using the K in the inner loop for
pitch angle control and the PD-type controller in the outer loop for altitude control
150 9 Classic Controller Design for the UAV
Fig. 9.7 Response of the altitude control system with a PID controller using Ziegler–Nichols gains
9.2 Classical Controller for the Longitudinal Motion 151
Fig. 9.8 Response of the altitude control system with a PID controller using the gains found in
tryouts
Fig. 9.10 Response of the PD-type altitude controller to a 100-m reference input
Fig. 9.11 Response of the system when the speed is trying to be held at 20 m/s (while the altitude
controller is working with the reference input)
u
w
deltae1
q
theta
Scope5-q
h
deltat
Vt
all
Uzunlamasina Denklemler
du/dt -K-
Step3 Derivative1 Kd
9.2 Classical Controller for the Longitudinal Motion
-K- s
Kp s+1
Washout filter
-K-
Krg
Scope8-theta
and saturators. The saturators are used for limiting the elevator angle and reference
pitch angle values. Although not shown in the diagram, the ıt input is determined
using the speed controller.
A controller that controls the speed using the thrust input, ıt, is designed. A P-type
controller is used. Both the results of the step response and the simulations show that
a P-type speed controller which uses ıt as the input is satisfactory. In order to find
the gain value K, the root locus method is used. The Vt /ıt transfer function does not
normally exist in the previously found functions; however, we can obtain it using
the equation below:
P C WW
VP t Vt D U U P C V VP (9.27)
We can find the Vt /ıt transfer function if we assume that, for a given linearization
point, Vut and Vwt are almost constant and v is zero. The other values can be determined
using the calculated values for the speed values, etc.
A sixth row to include the Vt transfer function is added to the system. The new
Alon matrix, denoted Alon _ vt , can be given as:
2 3
0:3356 1:3181 1:9276 9:6610 0 0
6 1:7916 3:9003 9:8215 1:7035 0 0 7
6 7
6 7
6 0:7020 3:5375 11:3920 0 0 07
Alon_vt D 6 7 (9.28)
60 0 1:0000 0 0 07
6 7
4 0:1736 0:9848 0 17:4865 0 0 5
1:234 3:5261 9:1667 3:8967 0 0
As explained previously, the measured value can be chosen as Vt and then the
related transfer function can be calculated. The transfer function for the thrust
actuator is modeled with a first-order transfer function with a time constant of 0.5 s:
The open-loop transfer function can be obtained by combining these two transfer
functions:
The root locus for the obtained Vt /ıt transfer function is as in Fig. 9.13. As can be
seen, the system is stable for gains smaller than 0.126 (K < 0:126). Firstly, the PID
gains proposed by the Ziegler method are checked. The determined critical values
are Kpu D 0.126 and ! u D 1.07 rad/s. The calculation for the gains are as follows:
2
Tu D D 5:8721s; (9.32)
!u
0:6 Kpu
Ki D D 0:0555: (9.35)
0:5 Tu
In Fig. 9.14, we see that the overshoot values and settling time are not good.
Hence, we designed a simpler controller using only the proportional gain value. The
system is already stable because of the roots of the transfer function. The response of
the system when Kvt D 0:013 is given in Fig. 9.15. We may obtain better responses
by increasing the gain value; however, the change in ıt to obtain the better responses
affects the overall system response. The response of the system when the gain value
is chosen as 0.025 is given in Fig. 9.16.
It is clear that the response quickens. However, the response of the overall control
system is also affected. In Figs. 9.17 and 9.18, the results for the altitude change
156 9 Classic Controller Design for the UAV
Fig. 9.14 PID controller response for the gains calculated using the Ziegler–Nichols method
Fig. 9.15 Response of the speed controller system to a 15-m/s reference input with only P gain
(Kvt D 0:013)
9.2 Classical Controller for the Longitudinal Motion 157
Fig. 9.16 Response of the speed controller system to a 15-m/s reference input with only P gain
(Kvt D 0:025)
Fig. 9.17 Change in altitude when the speed controller changes the UAV’s speed by 15 m/s
(Kvt D 0:015)
158 9 Classic Controller Design for the UAV
Fig. 9.18 Change in altitude when the speed controller changes the UAV’s speed by 15 m/s
(Kvt D 0:025)
when the speed value is set to a certain value are given. It can be understood that,
when the gain used for controlling Vt is increased, the change in altitude is also
increased. Changes of around 10 m in altitude over a short period, as we see for the
case Kvt D 0:025, is an undesired response for mini-UAV operations (Fig. 9.18).
Thus, a gain value that gives a quick response for controlling Vt and results in lower
altitude changes must be used.
In Fig. 9.19, the step response for the controllerwith Kvt D 0:015 is given. This is
the final gain value that we choose for the speed controller. For speed changes below
10 m/s, it may be increased to Kvt D 0:02, depending on the desired response.
The overall structure for the speed controller is given in Fig. 9.20.
The equations that we derived for the lateral motion were given in Chap. 8. In
addition to the classical states, we also used the heading equation, which was given
as P D r sec . In this section, we design a control system using the derived
equations and the matrices Alat and Blat .
For the lateral classical controller design, we prefer a method that consists of
inner loops for increasing the stability. A yaw damper is used to reduce the rudder
9.3 Classical Controller for the Lateral Motion 159
u
w
deltae1
q
theta
h
deltat
Vt
all
Longitudinal Dynamics1
dt
ue 2
-K-
s+2
KVt Throttle Actuator Scope-dt
reference value to zero when the controller for the heading is working. Changes in
the roll rate (p), roll angle (·), and heading angle (‰) are used to determine the
required input to the aileron. This method is a classical and effective method for
designing the lateral controller. However, during the design process, the coupling
160 9 Classic Controller Design for the UAV
between the aileron and rudder should be kept in mind to overcome the difficulties;
washout filter, aileron and rudder equations should be handled together, and actuator
equations should be included in the state equations while the controller is designed
(Fig. 9.21).
We can start the design with the p loop that has the least effect on the controller
[8]. p/ı a can be obtained using the lateral equations. The aileron actuator transfer
20
function is chosen as 20Cs .
Combining the given transfer functions, we can find p/ua (roll angle
change/control input):
The root locus drawn for this system is given in Fig. 9.22.
9.3 Classical Controller for the Lateral Motion 161
What we see from the root locus is that the root of the roll mode equation must
be moved to the left on the real axis to increase the stability of the system and
the damping factor of roll angle change. The roots of the characteristic equation
previously determined are:
From the root locus, we understand that the root may be moved to 4.03 by
choosing a gain value of 0.165. At the same time, this gain value faintly decreases
the Dutch roll mode damping factor. We also need to check the value of the spiral
mode root. The spiral mode root is unstable; however, when the gain is changed, the
root becomes 0.09 and the time for double amplitude can be found as:
log.2/
TD D 7:5s; (9.38)
0:09
which is good enough to ensure a normal working condition. Thus, the gain value
Kp can be chosen as 0.165.
The second loop to be closed is r/ı r (change in yaw angle by rudder input).
162 9 Classic Controller Design for the UAV
Starting from this step of the lateral controller design, we need to define output
matrix C and gain matrix K by taking the coupling effect in lateral motion and the
washout filter used in the r/ı r loop into account.
The state matrix (Alat ) of the lateral system is:
2 3
1:0502 1:9276 9:8215 9:6610 0
6 1:2213 1:9155 1:0096 0 07
6 7
6 7
Alat D 6 1:7255 0:0919 1:7198 0 07 (9.39)
6 7
40 1:0000 0:1763 0 05
0 0 1:0154 0 0
Here, ı a , aileron input, ı r , rudder input, and xw , washout, are added into the lateral
state matrix. Ta , Tr , and Tw are the time constants for the aileron actuator, rudder
actuator, and washout filter, respectively. The new lateral state matrix is:
2 3
1:0502 1:9276 9:8215 9:6610 0 0 1:8218 0
6 1:2213 1:9155 1:0096 0 0 8:3479 10:8560 0 7
6 7
6 1:7255 0:0919 1:7198 0 0 4:2400 2:1272 0 7
6 7
6 7
60 1:0000 0:1763 0 0 0 0 0 7
XD6 7 (9.41)
60 0 1:0154 0 0 0 0 0 7
6 7
60 0 0 0 0 20:0000 0 0 7
6 7
40 0 0 0 0 0 20:0000 0 5
0 0 1:0000 0 0 0 0 5
The values in the new B matrix are the rudder and aileron actuator constants. The
control inputs and states in this case are:
ua
ulat D (9.43)
ur
T
xlat D v p r ' ‰ ıa ır xw (9.44)
At this point, if we can define the necessary C and K matrices, then the feedback
can be achieved by calculating:
The matrix C includes feedback from xw , p, and ®. In B, one input from the
aileron and one input from the rudder are present. K must be determined at each
step using the appropriate feedback channel.
For example, for the first step, the gain matrix is:
0:1650 0 0
KD ; (9.46)
0 00
where xw , p, and ® are chosen as the outputs. In this way, for the first step, the
following matrix can be calculated:
2 3
00 00 00 00
60 0 00 00 0 07
6 7
60 0 00 00 0 07
6 7
6 7
60 0 00 00 0 07
BKC D 6 7 (9.48)
60 0 00 00 0 07
6 7
6 0 3:3000 0 0 00 0 07
6 7
40 0 00 00 0 05
00 00 00 00
The system must be updated using the newly calculated A matrix and finding r/ı r
by means of the xw state. The new lateral state matrix with the first feedback is:
164 9 Classic Controller Design for the UAV
2 3
1:0502 1:9276 9:8215 9:6610 0 0 1:8218 0
6 1:2213 1:9155 1:0096 0 0 8:3479 10:8560 0 7
6 7
6 1:7255 0:0919 1:7198 0 0 4:2400 2:1272 0 7
6 7
6 7
60 1:0000 0:1763 0 00 0 0 7
X1 D 6 7 (9.49)
60 0 1:0154 0 00 0 0 7
6 7
60 3:3000 0 0 0 20:0000 0 0 7
6 7
40 0 0 0 00 20:0000 0 5
0 0 1:0000 0 00 0 5
Fig. 9.23 Root locus diagram for the yaw angle loop
9.3 Classical Controller for the Lateral Motion 165
negative gains, a gain value of 0.405 makes the damping factor 0.43 and does not
affect the spiral mode root much. Therefore, Kr can be chosen as 0.4. After this
step, we must find ®/ı a , which is the equation that is necessary for ® feedback.
Similar steps to those introduced in the previous section can be used to find the
equation that relates ® to ı a and ua :
Fig. 9.24 Root locus diagram for the roll angle feedback loop
166 9 Classic Controller Design for the UAV
A gain value of K' 1 means that the spiral mode root becomes more stable
and the Dutch roll mode damping factor also stays high. In this case, the roll mode
root stays at a value of 5.5. This choice of K gives us a nonvibratory response. Also,
choosing a high proportional gain value decreases the steady-state error value.
In order to design the heading controller, we should use a last gain matrix to
determine the equation that relates heading to the aileron input. The last gain matrix
we found is:
0:1650 0 1
Klast D (9.53)
0 0:4 0
Recalculating the lateral A matrix and choosing the appropriate output, we can find:
The open-loop roots of the last loop are: 19.117 C 1.92i, 19.117 1.92i,
1.5285 C 4.79i, 1.5285 4.79i, 5.6534, 1.3728 C 2.4857i, and 1.3728
2.4857i.
The proportional gain value is chosen as 0.5 (Fig. 9.25). With this value, the roll
mode root value is improved for better response and the spiral mode is stabilized.
Also, the Dutch roll mode damping factor becomes 0.3. The response of the system
is given in Figs. 9.26, 9.27, 9.28, 9.29, and 9.30.
In the heading controller, a P-type controller is chosen because of its simple
design and quick response. As we can see from the figures, the controller directs the
aircraft to the reference heading value quickly, without significant change in the roll
angle.
In this system, the rudder and aileron actuator time constants are taken as 0.05 s
20
and their transfer functions are positive and negative sC20 , respectively.
The washout filter in the rudder loop has a time constant of T D 0.2 s, and its
s
transfer function is sC5 .
9.3 Classical Controller for the Lateral Motion 167
Fig. 9.26 Step response of the system to unit heading reference input
168 9 Classic Controller Design for the UAV
0.35
X: 26
Y: 0.35
0.3
0.25
Heading angle (rad)
0.2
0.15
0.1
0.05
0
0 50 100 150 200 250
t(s)
Fig. 9.27 Response of the system to a 0.35 radians reference heading input
0.08
0.07
0.06
0.05
r (rad/s)
0.04
0.03
0.02
0.01 X: 9
Y: 0.0004563
0
−0.01
0 50 100 150 200 250
t(s)
0.09
0.08
0.07
0.06
Phi – roll angle (rad)
0.05
0.04
0.03
0.02
0.01
0
0 50 100 150 200 250
t(s)
0.03
0.02
P change in roll angle(rad/s)
0.01
X: 23
Y: −6.287e-005
0
−0.01
−0.02
−0.03
0 50 100 150 200 250
t(s)
9.4 Conclusion
In this chapter, a flight control system is designed for UAVs using classical control
methods. Separate discussions for controlling different parameters, such as altitude
and heading, are given, and controller design procedures are explained in detail.
This method is the most common method for UAV control, and good results for
controlling both the lateral and longitudinal motion parameters also prove that.
References
1. Golnaraghi F, Kuo BC (2009) Automatic control systems, 9th edn. Wiley, Englewood Cliffs
2. Yechout TR, Morris SL, Bossert DE, Hallgren WF (2003) Introduction to aircraft flight
mechanics. AIAA education series, Reston
3. Williams W (2003) UAV handling qualities : : : ..you must be joking. Aerospace Sciences
Corporation Pty. Ltd., Elizabeth, Australia. Available online at: http://aerosciences.com.au/
hidden/UAV%20Handling%20Qualities%20Paper%20v1.pdf
4. Capello E, Guglieri G, Marguerettaz P, Quagliotti F (2012) Preliminary assessment of flying and
handling qualities for mini-UAVs. J Intell Robot Syst 65(1–4):43–61. doi:10.1007/s10846-011-
9565-5
5. Pamadi BN (2003) Performance, stability, dynamics, and control of airplanes. AIAA education
series, Reston
6. Nelson RC (1998) Flight stability and automatic control. McGraw-Hill, New York
7. Haugen F (2004) PID control. Tapir Academic Press, Trondheim
8. Stevens BL, Lewis FL (1992) Aircraft control and simulation. Wiley, Chichester
Chapter 10
Linear Quadratic Regulator Controller Design
10.1 Introduction
The first steps to design a control system for an UAV is defining the dynamic
model and finding the aerodynamic coefficients. The nonlinear model can then
be linearized to design a linear controller for the UAV [1]. The characteristic
values of different motions that show how the aircraft normally behaves must
also be investigated at this stage. Finally, the controller and filter can be designed
considering the required specifications.
In this chapter, the optimal control technique linear quadratic regulator (LQR)
is chosen to design a controller for the longitudinal motion of a small fixed-wing
type UAV. The Kalman filter technique is then applied to see how the controller
is affected by disturbances. The effectiveness of the controller with and without
the Kalman filter is also shown through simulations. Linear quadratic control is
an optimal control technique that is used for controlling the aircraft. This control
technique aims to decrease the energy that is used to control the aircraft. This
technique can be applied together with a filtering technique (Kalman filter, etc.)
in cases where some of the states are not available for measurement or when the
measurements are noisy.
Studies that include LQR control design for different types of aircraft have
been carried out in the literature [2]. In [3], an LQR method is used to control
the trajectory and mission paths of the autonomous helicopter. Nonlinear motion
dynamics are linearized at certain operating points and a linear model is obtained
by Taylor series expansion. By using the LQR methodology, the attitude of the
autonomous Puma helicopter is controlled. The study by Öner et al. [4] shows that
the LQR controller is quite effective in the vertical flight mode for all possible
yaw angles. In [5], the design procedure for a gain-scheduled LQR controller for
an autonomous airship is presented. Two types of control subsystems (lateral and
longitudinal) have been designed from the nonlinear six degrees of freedom (6DOF)
airship model to fulfill different goals (yaw as well as speed and position control).
In references [2–5], the LQR controller is not applied together with a filtering
technique.
In some studies, filtering techniques are used to analyze the effects of LQR design
and increase the effectiveness of the controller. The paper by Santoso et al. [6]
depicts the application of linear quadratic optimal control to the longitudinal flight
motion of a UAV which has elevon control only. The LQR controller is followed by
a Kalman filter-based estimator for unmeasurable states. The LQR controller is then
combined with the Kalman estimator using the separation principle to investigate the
feasibility of altitude control. The simulation results show improvements compared
with classical design counterparts in the sense that the combined approach offers
more design flexibility and is able to tolerate noisy environments.
Some advanced techniques, including robust and adaptive control techniques,
are also used in UAV control. For example, in one study by Johnson et al., the H-
infinity technique is used to design a velocity and altitude controller that follows a
determined model [7]. However, in this study, these advanced techniques are only
investigated to design a better controller and are not used.
In this chapter, an altitude controller is designed for a small UAV using the
optimal LQR control method with the Kalman estimator. The effectiveness of the
controller with the gains found by the LQR method is investigated under the effect
of disturbances and with the Kalman filter.
In order to design an optimal controller for the UAV, we use the linear quadratic
controller method in this chapter. This method depends on minimizing criteria using
the feedback u D kx.t/. The performance index that is used for minimizing both
the control effort and the states can be given as [8]:
Z1
1
JD xT Qx C uT Ru dt (10.1)
2
0
kxk D xT .t/x.t/ and x is an n 1 state vector. The given norm that depends
on the states can be used as a measure of the system response. The index also
includes u inputs to control the system. The Q and R matrices in J are called the
weight matrices. The weight matrices allow us to define the effect of each state and
control input in the controller response when we use the feedback for designing the
controller.
10.2 Linear Quadratic Optimal Controller 173
The weight given to the ith element of the matrix Q is a measure of the control effort
to control the related state. The bigger the element, the bigger the control effort used
for that state.
The matrix R is a p p, real, symmetric, and positive-definite matrix. Choosing
a positive-definite R is a must for an optimal solution.
If we can find the feedback u D kx.t/ that minimizes the J index for the system
given by xP D AxCBu, then we can solve the optimal control problem. In order to do
this, we need to relate the Lyapunov stability criteria with the quadratic performance
index.
The Lyapunov stability criteria deal with the behavior of the system response. When
the system leaves a steady state, whether it returns to a stable point or not can be
measured using the Lyapunov criteria. If a nonvibratory system response takes the
system very close to a point where it is stable for every x, then the system is called
Lyapunov stable. If the nonvibratory response of a system from a point x goes near
point xe (equilibrium point) for every x, then the system is called Lyapunov stable.
If the system response goes to xe for every x, then the system is called Lyapunov
asymptotically stable.
In order to control this criterion, the norm values kxk D xT .t/x.t/ that produce
a hypercubic space near the equilibrium point can be used. If the system response
stays in an area defined as S, it shows that the system is stable. In conclusion, using
the norm value, the total energy in the system can be calculated and the stability
condition can be checked. This method is called the Lyapunov method.
Generally, Lyapunov-type scalar functions can be defined as:
V.x/ D xT Px (10.3)
Lyapunov function can be built for the system and the function has partial derivatives
for every x variable, the stability condition can be checked using these partial
derivatives.
The positive definiteness requirements are:
V.0/ D 0 (10.5)
The following conditions are checked in order to understand whether the system
is stable or not:
P
(a) Asymptotically stable if V.x/ < 0 and x ¤ 0
(b) Generally asymptotically stable if the requirements in (a) are met and V.x/ >
1 when kxk > 1
Thus, a quadratic Lyapunov function can be defined and used for checking
system stability. This leads us to the linear quadratic optimal control method.
Z1
1
JD xT Qx C xT K T RKx dt (10.6)
2
0
The derivative of the Lyapunov function can be found, keeping in mind that
x(1) D 0:
1
VP .x.t// D xT Qx C xT K T RKx (10.8)
2
The Lyapunov function can also be written as:
1 T
V.x/ D x Px (10.9)
2
10.2 Linear Quadratic Optimal Controller 175
P 1 T
V.x/ D x Px C xT PPx (10.10)
2
Combining the equations found for the derivative of V(x), the following equation
can be solved:
The value of V(x) at 0 is equal to the performance index value for the system
trajectory starting from 0. Thus, the J index can be calculated as:
1 T
JD x .0/Px.0/ (10.12)
2
The values of matrix P that minimize the performance index J can then be calculated
as:
@ xT .0/Px.0/
D0 (10.13)
@kij
If the values of P that are optimal make the matrix P positive-definite, this
guarantees a stable loop defined by the Lyapunov function and we can calculate
the optimal gain. We may obtain the same result using the Riccati equation, which
is:
R D T ; (10.14)
from the definition of the matrix R. Using the given equation for R and deriving P
to find the optimal values, we obtain the optimal control input:
@ 1 T 1
K T BT P K T BT P D 0 (10.15)
@kij
K D R1 BT P (10.17)
The Riccati matrix equation can be found by using the previously determined value
of K [8].
AT P C PA PBR1 BT P C Q D 0 (10.19)
For designing a controller, both the first formula and the given Riccati equation can
be used. If the Riccati equation is to be used, the longitudinal matrices Alon , Blon and
lateral matrices Alat , Blat must first be determined. Then, the equation for finding the
matrix P can be determined or MATLAB functions can be used to find the required
values.
As an advantage, the LQR design can be done easily. However, a disadvantage
of this method is expecting that the values of the states are known exactly at every
moment. However, in general, these values are not known exactly and a filtering
technique is necessary. In this study, Kalman filtering is used to estimate the states
that are needed in the optimal controller.
Using the LQR method, a longitudinal controller for altitude and speed and a lateral
controller for heading is designed. The structure of the system allows us to use the
gain values found using the Riccati equation directly. The formulas given below
describe the motion of the vehicle around the stable point. Here, the xQ represents the
divergence from the required point and xd represents the reference input:
xQ D x xd (10.20)
As we can understand, if Axd D xdP holds using only u D kx.t/ without adding
anything into the state matrices, it will be enough to control the system optimally.
With choosing xd D [0; 0; 0; 0; hd], where hd is the reference altitude input, we can
10.3 Altitude and Speed Controller Design Using the LQR Method 177
reach the reference point only using state feedback in the given system. In short,
in cases where we have to control the state which has no direct effect on the other
states (where the A matrix columns related to the state are zeros), we can use the
given method directly.
The scheme of the controller is given in Fig. 10.1. In the scheme, hd values are
multiplied by K1 D [0;0;0;0;1] and added to u D kLQR x.t/ to find the real control
input that is used in the state model. The longitudinal states h, u, w, q, and are
scoped.
The controller is used to take the altitude value to the reference while holding the
other values constant. As we can see in Figs. 10.2, 10.3, 10.4, 10.5, 10.6, and 10.7,
the changes in the values other than the altitude goes to zero and the states hold their
values.
The weight matrices (Q, Q1 , and R) used in the controller are given below. The
change in U observed while the altitude value is changing can be controlled by
changing the corresponding weight matrix value. For this purpose, two different Q
matrices (Q and Q1 ) are used and the results are given. In R, the value for the elevator
is chosen to be high in order to use the elevator effort more effectively:
100 0
RD (10.25)
0 1
2 3
1000 0
60 1 0 0 0 7
6 7
6 7
Q D 60 0 1 0 0 7 (10.26)
6 7
4 0 0 0 100000 0 5
0000 10
2 3
10 0 0 0 0
60 1 0 0 0 7
6 7
6 7
Q1 D 6 0 0 1 0 0 7 (10.27)
6 7
4 0 0 0 100000 0 5
0 000 10
simout1
u_tilde To Workspace1
1
Matrix
Gain1 Kazanc Matrisi
x' = Ax+Bu
20 K*u K* u
y = Cx+Du u_tilde
h_d Saturation1 Plant
w_tilde simout2
delta_e To Workspace2
delta_t
q_tilde
-K-
Gain theta_tilde
10
simout
To Workspace
Fig. 10.2 LQR altitude controller response for a 20-m hd input (with Q)
2.5
2
Pitch angle (degrees)
1.5
0.5
X: 19
Y: 0.0002329
0
-0.5
0 5 10 15 20 25 30 35 40
t(s)
Fig. 10.5 LQR altitude controller response for a 20-m hd input (with Q1 )
10.3 Altitude and Speed Controller Design Using the LQR Method 181
4
Pitch angle (degrees)
1
X: 25
Y: 0.002847
0
-1
0 5 10 15 20 25 30 35 40
t(s)
Fig. 10.9 Response of the LQR speed controller to a 10-m/s reference input
The designed speed controller takes the aircraft to the reference speed while
holding the other states constant. The changes in different states when a 10-m/s
change in U is commanded are given in Figs. 10.9, 10.10, and 10.11.
The integrator is added to the system as the first state. The error of the speed is
fed to the system as a sixth state in the longitudinal equations. The new A and B
matrices built by adding the integrator are given below:
2 3
0 1:0000 0 0 0 0
60 0:3356 1:3181 1:9276 9:6610 07
6 7
6 7
60 1:7916 3:9003 9:8215 1:7035 07
Alon D6 7 (10.28)
60 0:7020 3:5375 11:3920 0 07
6 7
40 0 0 1:0000 0 05
0 0:1736 0:9848 0 17:4865 0
2 3
0 0
6 0:7436 6:8728 7
6 7
6 7
6 3:7855 0 7
Buzu D 6 7 (10.29)
6 47:9170 0 7
6 7
40 0 5
0 0
184 10 LQR Controller Design
Fig. 10.10 Change in altitude when the LQR speed controller takes the system to 10 m/s
-0.02
-0.04
-0.06
Pitch angle (rad)
-0.08
-0.1
-0.12
-0.14
-0.16
X: 18
Y: -0.1832
-0.18
-0.2
0 10 20 30 40 50 60 70 80 90 100
t(s)
Fig. 10.11 Change in pitch angle when the LQR speed controller takes the system to 10 m/s
10.4 LQR-Type Heading Controller 185
It is observed in the tryouts that the reference speed value can be reached quickly
without affecting the other states (altitude, pitch angle, etc.) when we choose the
diagonal elements of Qlon as 1000, 10, 1, 1, 100, and 30000. The diagonal elements
of Rlon are chosen as 100 and 100 to give the same effort to both elevator and thrust
control inputs:
2 3
1000 0 0 0 0 0
60 10 0 0 0 0 7
6 7
6 7
60 0 100 0 7
Quzu D6 7 (10.30)
60 0 010 0 7
6 7
40 0 0 0 100 0 5
0 0 0 0 0 30000
100 0
Ruzu D (10.31)
0 100
For the lateral controller, an integrator may also be added to the system; however, the
structure of the system tells us that, for controlling the heading using the reference
input xd D [0; 0; 0; 0; d], the normal lateral matrices are sufficient. For the case
where heading is the input, the following equation holds:
Thus, as stated, we can use the lateral state and control distribution matrices to
calculate the result of the Riccati equation and the LQR gain values. The Q and R
matrices for the lateral LQR controller are chosen as follows:
2 3
10 0 0 0
6 0 10 0 0 0 7
6 7
6 7
Q D 6 0 0 100 0 0 7 (10.33)
6 7
4 0 0 0 10000 0 5
00 0 0 1000
100 0
RD (10.34)
0 100
The scheme of the lateral LQR heading controller is given in Fig. 10.12. As we
can understand from the results (Figs. 10.13, 10.14, and 10.15), the LQR controller
is working optimally and takes the UAV to the desired heading with little changes
occurring in ı a and ı r . The vibrations in states may be lessened by choosing
appropriate Q and R weight matrices.
186
Matrix Matrix
Gain1 Gain
x' = Ax+Bu
0.35 K*u K*u
y = Cx+Du beta_tilde
Psi_d Saturation1 Plant
p_tilde
delta_a
delta_r
r_tilde
-K-
Gain Phi
Psi
10
0.4
0.35
X: 26
Y: 0.35
0.3
Heading angle(rad)
0.25
0.2
0.15
0.1
0.05
0
0 50 100 150 200 250 300
t(s)
0.16
0.14
0.12
r(change in yaw angle(rad/s))
0.1
0.08
0.06
0.04
0.02
X: 37
Y: -2.034e-007
0
-0.02
0 50 100 150 200 250 300
t(s)
Fig. 10.14 Change in r (yaw angle change, rad/s) when a d D 0.35 radians input is given
188 10 LQR Controller Design
5
Phi(roll angle(rad))
1
X: 30
Y: 0.0004451
0
-1
0 50 100 150 200 250 300
t(s)
Fig. 10.15 Change in roll angle (·) when a d D 0.35 radians input is given
In this section, it is assumed that the measured states are correct at any time and the
designed LQR is fed with them to control the UAV. However, noise is always present
and the system may not give accurate results under normal operating conditions. For
this purpose, a Kalman filter that has the capability to estimate the real values of the
states is incorporated in the study. The results of the system with a Kalman filter
and LQR and only LQR are then compared to show the effectiveness of the Kalman
filter. The system model is discretized using the Euler method. Firstly, the nominal
model is tested with 20 m altitude and 20.02ı heading inputs, with the optimal
gain values calculated previously. Then, the system with noise is simulated. The
response of the system with the controller under the disturbances is tested with and
without the Kalman filter and the effectiveness of the designed LQR and Kalman
filter system is shown [9].
The Kalman filter uses state equations (state space matrices) and initial values
to calculate the residue and gain values and to estimate the real signal value. The
steps of the Kalman filter can be explained using the linear discretized state and
measurement equations:
In the state equation, X(k) is the state vector of the system, A is the system
transition matrix, u(k) is the input vector, B is the control distribution matrix, w(k)
is the random Gaussian noise vector (system noise) with zero mean and known
covariance structure, and G is the transition matrix of the system noise. In the
measurement equation, y(k) is the measurement vector, H is the measurement
matrix, and v(k) is the measurement noise vector with zero mean and known
covariance structure.
There is no correlation between the system noise w(k) and the measurement noise
v(k). The covariance matrices for the w(k) and v(k) vectors are given by:
E w.k/wT .j/ D Q.k/ı.kj/
E v.k/v T .j/ D R.k/ı.kj/
Here, E is the expected value operator and ı(kj) is the Kronecker symbol.
The optimum linear Kalman filter that estimates the state vector of the system
(10.35) is expressed using the following recursive equations system:
Equation of the extrapolation value:
The Kalman filter tries to estimate the real signal from the signal with distur-
bance, which has a Gaussian distribution, using the described steps and decreasing
the value between two signals [10].
The model of longitudinal motion can be discretized using the Euler approach for
a sampling time of dt. Thus, the new discretized A and B matrices to be used in the
filtering approach can be found by A1 D A dt C I and B1 D B dt.
The discretized UAV model can be given as follows:
2 30 1
0:997 0:013 0:02 0:097 0 u.k/
6 0:018 0:961 0:098 0:017 0 7 B w.k/ C
6 7B C
6 7B C
X .k C 1/ D 6 0:007 0:035 0:886 0 0 7 B q.k/ C
6 7B C
40 0 0:01 1 0 5 @ .k/ A
0:0017 0:01 0 0:175 1 h.k/
2 3
0:007 0:07
6 0:04 0 7
6 7
6 7 ıe .k/
C 6 0:48 0 7 C Gw.k/ (10.43)
6 7 ıt .k/
40 0 5
0 0
Change in altitude
25
LQR without disturbances
LQR(with disturbances)
20 LQR(with Kalman filter
15
altitude,h
10
-5
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
step
Change in altitude
20.15
LQR without disturbances
LQR(with disturbances)
20.1 LQR(with Kalman filter
20.05
altitude,h
20
19.95
19.9
19.85
Fig. 10.19 Change in h in the LQR controller (zoomed, 5,000 steps D 500 s)
10.5 LQR Controller with the Kalman Estimator 193
0.2
error,kalman-normal
0.1
-0.1
-0.2
-0.3
-0.4
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
step
Fig. 10.20 Differences between Kalman filter estimation and normal values for altitude
0.5
error,kalman-normal
-0.5
-1
-1.5
-2
-2.5
-3
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
step
Fig. 10.21 Differences between Kalman filter estimation and normal values for forward velocity
194 10 LQR Controller Design
0.4
error,disturbance-normal
0.2
-0.2
-0.4
-0.6
-0.8
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
step
Fig. 10.22 Differences between readings with disturbances and normal values for altitude
3
error,disturbances-normal
-1
-2
-3
-4
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
step
Fig. 10.23 Differences between readings with disturbances and normal values for forward
velocity
10.5 LQR Controller with the Kalman Estimator 195
0.035
0.03
0.025
0.02
0.015
0.01
0.005
0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000
Fig. 10.24 Change in the fifth diagonal element of the covariance matrix (P1 )
A similar study for the lateral controller is accomplished. For this system, heading
reference is the input. If we choose Xd D [0; 0; 0; 0; d] and add the necessary LQR
gain input, we can build a model for the system. The MATLAB code for the system
is similar to that for the longitudinal case.
The change in r (yaw angle change, rad/s) and heading angle are given in the
graphs. The simulation results are given in Figs. 10.25, 10.26, 10.27, 10.28, 10.29,
10.30, 10.31, and 10.32. The results show the effectiveness of the Kalman filter.
The differences between the normal and noisy values for heading and yaw angles
are also given.
The presented simulation results (Figs. 10.25, 10.26, 10.27, 10.28, 10.29, 10.30,
10.31, and 10.32) show that, in case of disturbances, using a Kalman filter to
estimate the values of the lateral states clearly increases the effectiveness of the
LQR controller.
196 10 LQR Controller Design
Heading angle
25
20
15
heading,psi
0
0 50 100 150 200 250 300 350 400 450 500
step
Heading angle
20.3
20.2
LQR without disturbances
20.1 LQR(with disturbances)
LQR(with Kalman filter
20
heading,psi
19.9
19.8
19.7
19.6
19.5
19.4
390 395 400 405 410 415 420 425 430 435
step
Fig. 10.26 Heading controller response to a 20ı input (zoomed, 500 steps D 50 s)
10.5 LQR Controller with the Kalman Estimator 197
25
change in roll angle,r
20
15
10
-5
0 50 100 150 200 250 300 350 400 450 500
step
0.1
0.05
errors,kalman-normal
-0.05
-0.1
-0.15
-0.2
0 50 100 150 200 250 300 350 400 450 500
step
Fig. 10.28 Differences between Kalman filter estimations and normal system for r
198 10 LQR Controller Design
0.2
0.15
errors,kalman-normal
0.1
0.05
-0.05
-0.1
-0.15
-0.2
0 50 100 150 200 250 300 350 400 450 500
step
Fig. 10.29 Differences between Kalman filter estimations and normal system for heading angle
2
errors,disturbance-normal
1.5
0.5
-0.5
-1
0 50 100 150 200 250 300 350 400 450 500
adim
Fig. 10.30 Differences between readings with disturbances and normal system for r
10.5 LQR Controller with the Kalman Estimator 199
0.4
errors,disturbance-normal
0.2
-0.2
-0.4
-0.6
-0.8
0 50 100 150 200 250 300 350 400 450 500
step
Fig. 10.31 Differences between readings with disturbances and normal system for heading
0.025
0.02
0.015
0.01
0.005
0
0 50 100 150 200 250 300 350 400 450 500
Fig. 10.32 Change in the diagonal element (fifth element) of the P2 covariance matrix
200 10 LQR Controller Design
Optimal controllers using the LQR method to control the longitudinal and lateral
flight dynamics of small UAVs have been designed and examined in this chapter.
The effects of disturbances on measurements are taken into account and the
effectiveness of the Kalman filter in obtaining correct measurements and achieving
the desired control level are shown using the controller designed for the system.
A Kalman estimator is designed and added to the system and the response of the
controller with and without disturbances and the Kalman filter is tested. The results
reveal the effectiveness of the Kalman filter and the LQR controller.
The simulations show that, in case of disturbances, using a Kalman filter to
estimate the values of the states clearly increases the effectiveness of the LQR
controller. In conclusion, it is found that the LQR method with a Kalman estimator
is effective in controlling the longitudinal and lateral flight dynamics of UAVs and
can be used for such applications.
References
1. Jang JN (2003) Longitudinal stability augmentation system design for the dragon fly UAV
using a single GPS receiver. In: Proceedings of the AIAA guidance, navigation, and control
conference and exhibit, Austin, TX, August 2003. AIAA paper 2003-5592
2. Kinoshita T, Imado F (2006) A study on the optimal flight control for an autonomous UAV.
In: Proceedings of the IEEE 2006 international conference on mechatronics and automation
(ICMA), Luoyang, China, June 2006, 43(38)
3. Franko S (2009) LQR-based trajectory control of full envelope, autonomous helicopter. In:
Proceedings of the world congress on engineering (WCE 2009), London, UK, July 1–3 2009,
vol I
4. Öner KT, Çetinsoy E, S{r{moğlu E, Hancer C, Ayken T, Ünel M (2009) LQR and SMC
stabilization of a new unmanned aerial vehicle. World Acad Sci Eng Technol 34:373–378
5. Masar I, Stöhr E (2011) Gain-scheduled LQR-control for an autonomous airship. In: Proceed-
ings of the 18th international conference on process control, Tatranská Lomnica, Slovakia,
14–17 June 2011, pp 197–204
6. Santoso F, Liu M, Egan GK (2007) Linear quadratic optimal control synthesis for an UAV. In:
Proceedings of the 12th Australian international aerospace congress (AIAC12), Melbourne,
Australia, 16–22 March 2007
7. Johnson MD, Calise AJ, Johnson EN (2003) Evaluation of an adaptive method for launch vehi-
cle flight control. In: Proceedings of the AIAA guidance, navigation, and control conference
and exhibit, Austin, TX, August 2003
8. Anderson BDO, Moore JB (1990) Optimal control: linear quadratic methods. Prentice Hall,
Englewood Cliffs
9. Hajiyev C, Vural SY (2010) LQR controller with Kalman estimator applied to UAV longitu-
dinal dynamics. International workshop on unmanned air vehicles UVW2010, Istanbul, 10–12
June 2010
10. Hajiyev C (1999) Radio navigation. Istanbul Technical University, Istanbul (in Turkish)
11. Sage AP, Melsa JL (1971) Estimation theory with applications to communications and control.
McGraw-Hill, New York
Chapter 11
Fuzzy Logic-Based Controller Design
The accuracy of the proposition is a value between 0 and 1, and depends on the
relation between x and A.
For example, the relationship between throttle and vehicle speed can be defined
by a fuzzy rule as IF throttle openness is high THEN vehicle speed is high. A
constructed fuzzy model can be named by the fuzzy rules included in it. The output
structure in a fuzzy proposition determines the fuzzy rule type. Four different types
are present:
1. Mamdani-type fuzzy rules
2. Singleton-type fuzzy rules
3. Takagi–Sugeno-type fuzzy rules
4. Tsukamoto-type fuzzy rules
This type of rule includes semiqualified knowledge about the variables, and they
can be given as:
IF x Ai THEN y Bi .
In this rule, x is the linguistic input and the Ai s are antecedent linguistics terms.
Likewise, y is the linguistic output and the Bi s are consequent linguistic terms. The
linguistic variables x and y and the linguistic terms Ai and Bi in the rules are fuzzy
sets that are defined in their own definition sector.
For example, the membership functions for the antecedent and consequent fuzzy
sets can be given as:
x 2 X Rp (11.1)
y 2 Y Rq (11.2)
11.1 Fuzzy Logic-Based Systems 203
The fuzzy sets Ai are fuzzy spaces in the antecedent space, where related
consequent propositions are real. The linguistic variables Ai and Bi are generally
chosen between the terms little, very little, moderate, big, very big, etc. If we define
these sets as A and B, we can define the Ai and Bi as subsets Ai 2 A and Bi 2 B.
For example, a heater system with constant gas flow where the heat received
depends on the given O2 input may be described as a fuzzy model:
The input is the O2 flow and the output is heating power. The sets for input and
output linguistic terms can be given as:
The relationship between the inputs and outputs can be given as:
Rule 1: IF O2 flow is LOW THEN the heating power is low
Rule 2: IF O2 flow is MODERATE THEN the heating power is high
Rule 3: IF O2 flow is HIGH THEN the heating power is low
The linguistic terms are defined in the fuzzy system by membership functions.
For the exemplified inputs and outputs, the membership functions given in Fig. 11.2
can be used.
0 1 2 3 0 25 50 75 100
Singleton-type fuzzy rules are a special type of Mamdani rules. The output fuzzy
set is a single set.
The singular membership function can be given as:
1; x D x
.x/ D (11.7)
0; IF NOT
In Mamdani-type IF–THEN rules, the antecedent and consequent parts define the
system with fuzzy propositions. In the consequent part of Takagi–Sugeno (T-S) type
fuzzy models, a mathematical function is present. Thus, this type of model is built
with both linguistic and mathematical expressions.
The rules in T-S models are as given below:
IF x is Ai THEN yi D fi (x), i D 1, 2, : : : .., r.
Unlike the Mamdani-type model, the input x is exact. All the f (x) functions are
of the same type; however, the parameters within them are different. For simplicity,
a linear function can be used as an output function.
0.3
0.3
1
0
min 0 max
0.3 and 0.5–0.3
0.3
1
0
min 0 max
0.5
0.5
0
min 0 max
Ri .x; y/ D Ai .x/ ^ Bi .y/ D min .Ai .x/; Bi .y// (11.8)
The minimum operator is applied for all possible x–y combinations in the
Cartesian space of X and Y. The fuzzy R relation that defines the model is given
by the unification of r fuzzy rules that are in the Mamdani rule base:
r
R D [ Ki (11.9)
iD1
If we choose the maximum operator for unification, the membership function for
the R relation becomes:
R .x; y/ D max .Ai .x/ ^ Bi .y// D max .min .Ai .x/; Bi .y/// (11.10)
lir lir
For an example, we can take the case that is shown in Fig. 11.3, where one rule
is presented:
IF “altitude is low (0.3 shows low altitude membership value)” and “vertical
velocity (w) is moderately low (0.5 shows low mod. membership value)” THEN
elevator angle output is “low mod”
We can use min.–max. relations to obtain the efficacy factor of the output. In this
case, the lowest value of the inputs is taken and the area of the membership function
for low or moderate elevator angle output is calculated. Mostly, more than one rule
is effective and, in this case, every single output value can be used to find an average
value:
206 11 Fuzzy Logic-Based Controller Design
As observed in the other controller designs, the classical methods require knowledge
of the system model. However, a fuzzy logic controller may be designed using
previously known input–output relations and expert knowledge, without actually
having the system model.
A Mamdani-type controller is based on using knowledge-based rules that are
useful to take the system to previously determined requirements. Mamdani proposed
that, starting from requirements and expert knowledge, a controller based on IF–
THEN rules and membership functions for input and outputs can be designed. In
this study, Mamdani-type controllers are used [3].
In this case, the error and the derivative of the error are chosen as inputs. This
type of controller is similar to a proportional-derivative (PD) controller. For altitude,
velocity, and heading controllers, a PD-type fuzzy logic controller is designed.
The PD controller scheme used in this study is shown in Fig. 11.4. Error (e(t))
and the change in error (e(t) e(t 1)) is used to calculate the output u. Triangular-
type membership functions and PD-type inputs are used in the controller.
The rule base for both Vt–dt and h–™ref are as shown in Table 11.1. For an
inner loop controller, the one that was previously used in the classical method is
used (P controller) and the outer loop controllers are designed based on fuzzy logic
rules.
The rules are defined as below:
-K-
-K- u fcn y
K1
K3
FLC PLANT
-K-
1
K2
z
Unit Delay
1. IF the error in altitude is negatively high and the change in the error of altitude is
negatively high THEN the required output for pitch angle is negatively high
A total of 25 rules are defined. In these rules, the value of the error itself is
considered to be important. Then, the coefficients K1 and K2 may be changed
for different operating conditions. Also, K3 can be changed for the required
output. Decreasing the input and output coefficients means widened membership
functions for the input functions and centered membership functions for the output
membership functions [2].
An area-centered structure is used as an output defuzzifier.
Firstly, longitudinal fuzzy controllers are designed. For every controller, error and
error change membership functions are defined. For these calculations, MATLAB
software is used [4]. The membership functions for the altitude controller are given
in Figs. 11.5, 11.6, 11.7, and 11.8.
The membership functions for the change in error are similar; however, the range
is between 100 and C100. The error ranges in altitude, change of error in altitude,
and the output pitch angle are 400/C400 m, 100/C100, and C0.1745/0.1745
radians (10ı ), respectively. Seven output membership functions are defined.
The same rules and functions are used for the velocity controller. The range
of error and change in error membership functions in the velocity controller are
20/C20 m/s and 20/C20, respectively. For the dt output, the range is between
1 and C1. In the dt membership functions, zero is showing the actual dt input.
When the controller can be designed using expert knowledge and considering
different conditions, no input–output relationship may be necessary. However,
optimizing the system requires a detailed study on membership functions and, in
this study, the coefficients are used to develop a better system.
The results show that a Mamdani-type PD controller with homogenously dis-
tributed membership functions gives good results when the coefficients are deter-
mined carefully. The advantages of the fuzzy logic controller are that it can
be designed without knowing the dynamics of the system and it can be easily
changed using coefficients instead of the membership functions themselves. Using
208 11 Fuzzy Logic-Based Controller Design
Fig. 11.6 Change in error membership functions for the altitude controller
11.2 Fuzzy Controllers 209
In the lateral controller design (heading controller), the rules and membership
functions types are the same as for the longitudinal controller. The input and
output functions are defined in radians. The range for error and the change in error
Fig. 11.9 Results from the fuzzy logic controller using a 100-m altitude input
11.2 Fuzzy Controllers 211
10
6
de(degrees)
X: 49.2
Y: 1.661
2
-2
0 20 40 60 80 100 120 140
t(s)
Fig. 11.11 Change in elevator angle when the altitude input is 100 m
212 11 Fuzzy Logic-Based Controller Design
0.3
0.25
0.2
Pitch angle (radian)
0.15
0.1
X: 50
0.05 Y: 0.02881
-0.05
0 20 40 60 80 100 120 140
t(s)
Fig. 11.12 Change in pitch angle when the altitude input is 100 m (in radians)
are 5/C5 radians and 1/C1, respectively. For the output functions, 20/C20ı
(0.349/C0.349 radians) is chosen as the range. Again, the output coefficient is
used to scale the membership functions.
The membership functions are given in Figs. 11.18 and 11.19, and the heading
controller surface diagram is presented in Fig. 11.20.
The rules for the lateral controller are similar to those used for the longitudinal
controller [3].
The same types of membership functions are also used for the change in error.
11.2 Fuzzy Controllers 213
Fig. 11.13 K3 output coefficient calculation using the signal constraint method
Fig. 11.14 14-m/s velocity input result of the fuzzy logic controller
214 11 Fuzzy Logic-Based Controller Design
12
10
4
de(degrees)
2 X: 41.2
Y: -0.2067
0
-2
-4
-6
-8
0 20 40 60 80 100 120 140
t(s)
0.25
0.2
0.15
Pitch angle(radians)
0.1
0.05
X: 58
Y: -0.0007956
0
-0.05
-0.1
0 20 40 60 80 100 120 140
t(s)
The results of the lateral fuzzy logic controller are given in Figs. 11.21, 11.22,
and 11.23.
The heading controller and lateral fuzzy logic controller schemes are shown in
Figs. 11.24 and 11.25, respectively.
The results show that the lateral controller also works well if good coefficients
can be determined.
0.4
0.35
X: 70.6
Y: 0.3498
0.3
Heading angle(radians)
0.25
0.2
0.15
0.1
0.05
0
0 50 100 150 200 250
t(s)
Fig. 11.21 Response of the fuzzy logic controller to a 0.35 radians (20.02ı ) input
218 11 Fuzzy Logic-Based Controller Design
0.2
0.15
0.1
P(roll angle change-rad/s)
0.05
X: 36.6
Y: -0.000292
0
-0.05
-0.1
-0.15
-0.2
0 50 100 150 200 250
t(s)
0.015
0.01
Rudder input(dr-rad)
0.005
X: 46
Y: 1.234e-005
0
-0.005
-0.01
-0.015
0 50 100 150 200 250
t(s)
Fig. 11.23 Change in dr (rudder input) for a 0.35 radians heading input
11.4 A Comparison of Flight Controllers for Unmanned Aerial Vehicles (UAVs) 219
Kfl
Phi angle
+
− Kfl +
−
Heading Ref1
Saturation
Heading Controller
1
z
Scope1-Heading error Unit Delay
Heading feedback
−K
Kphi
+ + −K
+ − −
−
Saturation1 Kp
Heading Ref1 Saturation
Heading controller
μa
20
Aileron actuator
s+20
X>
deltaa R
1 [
z −20 Phl
Unit Delay s+20 deltar Psl
Scope2- roll error
Saturation2 all >
Rudder actuator
Lateral Dynamics
+
− Kr
Rudderref s
−K
s+5
Washout filter
Three types of flight control systems (classical control methods, optimal control
methods, and fuzzy logic technique) for a small UAV are compared in this section.
Classical control system design is the most common system for these kinds of
vehicles, and it has been proven to be a good way of controlling aircraft. On the
other hand, there are several advantages of other control systems. Linear quadratic
regulator (LQR) controllers can be designed quickly and fuzzy logic controllers
can be designed without having previous knowledge of the system dynamics. The
obtained results may be utilized in the design of UAV flight control systems. The
response values of the controllers are given in Tables 11.3, 11.4, and 11.5.
220 11 Fuzzy Logic-Based Controller Design
Table 11.4 Linear quadratic LQR controller Settling time Rise time
regulator (LQR) controller
response data Altitude (h) 6s 3.05 s
Total vel. (Vt ) 7.82 s 1.6 s
Heading ( ) 14.3 s 6.7625 s
Table 11.5 Fuzzy logic Fuzzy logic controller Settling time Rise time
controller response data
Altitude (h) K3 D 10 37 s 15 s
Total vel. (Vt ) K3 D 30 37 s 15.32 s
Heading ( ) 60 s 26 s
All the given values are read from simulations for similar inputs. The PID
controller works well as altitude and heading controllers, but the effect of altitude
and velocity on each other makes it necessary to use a slow-response controller for
the velocity. The LQR controller responses are quick; however, the control inputs
are high, even though the cost matrices are adjusted carefully. Besides, the values
that are fed to the LQR controller must be filtered, as previously described in the
Kalman filtering section. The fuzzy logic controller works well, even though no
optimization technique is used to optimize the system and it was designed without
dynamic model knowledge. In the fuzzy logic controller, a careful arrangement of
membership functions may lead to better results.
In conclusion, it is found that all three types of controllers are suitable for use
as small UAV controllers. It is observed that, although the PID controller responses
are good, an LQR controller and a fuzzy logic controller with membership function
optimization can also be used as a controller in UAV systems.
In this chapter, fuzzy logic-based longitudinal and lateral controllers have been
designed and different UAV flight control methods were compared. The fuzzy logic
controller works well, even though no optimization technique is used to optimize
the system and it was designed without dynamic model knowledge. In the fuzzy
logic controller, a careful arrangement of membership functions may lead to better
results. For fuzzy logic controllers, in general, simulations are used to determine if
the system is stable with the designed controller or not. In this case, the simulations
showed that the considered system is stable.
References 221
Comparison of different UAV flight control methods shows that the PID con-
troller responses are good and, also, LQR and fuzzy logic controllers with mem-
bership function optimization are effective in controlling the UAV. Designing better
controllers without neglecting the effect of each motion on one another, simulating
the system under different inputs and conditions, and trying to design more efficient
and robust controllers may be potential research topics.
References
1. Klir GJ, Yuan B (1995) Fuzzy sets and fuzzy logic theory and applications. Prentice Hall, Upper
Saddle River
2. Passino KM, Yurkovich S (1998) Fuzzy control. Addison-Wesley, Menlo Park
3. Vural SY, Hajiyev C (2014) A comparison of longitudinal controllers for autonomous UAV. Int
J Sustain Aviat 1(1):58–71
4. Hines JW (1997) MATLAB supplement to fuzzy and neural approaches in engineering. Wiley,
New York
About the Authors
Halil Ersin Soken was born in Istanbul, Turkey. He graduated “Summa cum Laude”
from the Astronautical Engineering Department of Istanbul Technical University in
2007. In 2008, he completed his second degree in Aeronautical Engineering. He
earned his Ph.D. in Space and Astronautical Science from The Graduate University
for Advanced Studies (Japan) in 2013. During his Ph.D. studies, he was a visiting
researcher at Aalborg University (Denmark) for 3 months.
Dr. Söken worked as a research assistant at Istanbul Technical University (2010)
and Japan Aerospace Exploration Agency – JAXA (2010–2013). Currently, he is
a research associate at the Institute of Space and Astronautical Science (ISAS) of
JAXA.
His research interest is Kalman filter-based estimation techniques and satellite
attitude determination and control, specifically for small satellites.
S. Yenal Vural was born in Adana, Turkey. He is a Ph.D. student at the Aeronautical
Engineering Department at Istanbul Technical University (ITU). He received his
B.Sc. degree in Mechanical Engineering from ITU in 2002 and his M.Sc. degree in
Aerospace Engineering from the same university in 2008.
His research interests include unmanned aerial vehicles (UAVs), flight control,
Kalman filtering, and fault-tolerant flight control.
Index
Innovation (cont.) L
based KF adaptation, 67 Lapse rate, 37, 38
channel, 90 Lateral
covariance, 54, 55, 75, 97 characteristic, 137, 138
matrix, 75 classical controller, 158–169
process, 113, 126 controller, vii, 141, 159, 162, 185, 195,
sequence, 58, 67, 73, 74, 75, 78, 116, 210, 212, 215, 216, 220
117, 189 dynamics, 60, 79, 219
Innovation-based adaptive estimation (IAE), equations, 18, 19, 61, 132, 134, 136, 141,
67, 71, 72 160
Input flight dynamics, 200
axes, 47 fuzzy logic controller, 210–215, 219
error, 110 LQR controller, 185, 186
function, 207, 208 LQR heading controller, 185
Kalman observers, 110 matrices, 176, 185
–output gradient, 47 motion, 130, 136, 138, 139, 158–169
–output relationship, 110, 206, 207 stability, vi, 132, 138–139
–output test, 111 state matrix, 162, 163
values, 111, 177 states, 105, 132, 195
vector, 60, 189 variables, 131
Integrated Latitude, 27, 42, 43
ADS and INS/GPS, 113 Lift, 6
ADS/GPS navigation system, 111–114 Likelihood of system failures, 95
INS/GPS system, 68, 72, 111 Linear
R and Q-adaptation, 99 acceleration, 27
system, 111 controller, 171
Interacting multiple model (IMM), 110, discrete dynamic system, 52
111 discrete Kalman filter, 52–59
filters, 52, 58
filter theory, 52
K Kalman filter, v, 85, 189
Kalman quadratic
estimation, 62, 197, 198 estimation (LQE), 51
estimator, 172, 188–200 optimal control, 172–176
filter, v, vi, vii, 25, 44, 51–70, 71–78, 80, regulator (LQR), vi, vii, 171–200, 219,
85, 87, 88, 95, 96, 97, 99, 100, 101, 220, 221
102, 104, 107, 109–126, 171, 172, system, 51, 92, 110
188, 189, 190, 193, 195, 197, 198, Linearization, 52, 110, 129, 130, 154
200 Linearized
filter application, 76 equations, v, 18–23, 130, 134, 139
filter gain, 73, 88 longitudinal equations, 18, 145
filtering algorithm, 110 Local filter, 72, 114, 115
filtering equations, 113 Longitude, 42, 43
filtering structure, 109 Longitudinal
filtering technique, vi, 71, 111, 115, 126, channel, 105
176, 190, 220 characteristic equation, 137
filter innovation sequence, 75 control, 101, 105
filter normalized innovation, 117 control derivatives, 105
gain, 54–57, 74, 75, 76, 78, 91 control distribution matrix, 101
gain matrix, 76, 91 controller, 144, 176, 210
observers, 110 dynamics, 144, 159, 190
Kronecker equations, 18, 20, 110, 131, 134, 144, 183
delta function, 61 flight dynamics, 200
symbol, 52, 189 fuzzy controllers, 207
Index 229
LQR controller, 190 noise covariance matrix, 61, 72, 73, 75, 96
LQR with Kalman estimator, 190–195 noise scale factor (MNSF), vi, 73–76,
matrices, 176 78–80, 96
motion, 137, 138, 139, 144–158, 170, 171, noise vector, 52, 189
190 residual, 54
stability, 131, 137–139 Medium altitude long endurance (MALE), 5,
states, 103, 131, 135, 177, 190 44
Loop transfer function, 154 Membership function, vii, 202–208, 212, 216,
Loss of accuracy, 45, 46 220
LQR Microelectromechanical systems (MEMS)
altitude controller, 176–180 inertial sensors, 27–29
control design, 171 technology, 27
gain, 185, 195 Minimum
heading controller, 185 error variance, 51
Lyapunov mean square error estimator, 56
asymptotically stable system, 173, operator, 204, 205
174 MMAE. See Multiple model-based adaptive
criteria, 173–174 estimation (MMAE)
function, 173–175 MMNSF. See Multiple, measurement noise
method, 173 scale facto (MMNSF)
stability criteria, 173–174 MNSF. See Measurement, noise scale factor
stability function, 174 (MNSF)
Model-based adaptive estimation, 66, 71, 72
Moment of inertia, 13, 14
M Monitoring, 5, 6, 47, 117
Mach number, 30–32 Moving
Magnetic measurements, 40–41 frame, 15
MALE. See Medium altitude long endurance window, 72, 75
(MALE) Multi-input/output feedback control, 95
Malfunction, vi, 45, 69, 71, 73–76, 81, 82, 85, Multiple
90–92, 96, 98 adaptive factor, 97–98, 105, 107
Mamdani fading factor, 96
method, 204 measurement noise scale facto (MMNSF),
rule base, 205 73, 75–92
rules, 204, 205 model, 66, 71, 72, 110
type controller, 206 sensor fault detection, 111
type fuzzy rule, 202–203 Multiple model-based adaptive estimation
type IF-THEN rule, 204 (MMAE), 66, 71, 72
type PD controller, 207
Mathematical model, 45, 52, 58, 79
Maximum eigenvalue, 78 N
Mean Navigation
absolute values of error, 103–106 sensor, 114, 115
chord length, 20 systems, v, 25–47, 111
of the innovation, 95, 116, 117 NLGA. See Nonlinear equation, geometric
Measurement approach (NLGA)
channel, 62, 71, 73, 78, 79, 80, 82 Noise increment, 45, 46, 76, 81–82, 84–87, 91,
data, vi, 45, 64, 109, 115, 126 105, 106
equation, 37, 46, 52, 188, 189 Nonlinear equation
faults, v, vi, 45–47, 62, 71, 76, 81, 82, 92 estimation, 96
matrix, 54, 60, 68, 99, 113, 189 fuzzy logic based system, 215
noise, 46, 52, 54, 56, 57, 60, 61, 66, 67, 69, geometric approach (NLGA), 110
72, 73, 76, 81, 84–87, 91, 95, 96, motion dynamics, 171
189 system, 92, 109–111, 129
230 Index
Regulator Simultaneous
design, 171–200 adaptation procedure, 96
Relative altitude, 35 localization and mapping (SLAM), 44
Residual Single
based adaptive estimation (RAE), 68–69, adaptive factor, 72, 96
71, 72 fading factor, 96
sequence, 72 measurement noise scale factor (SMNSF),
series, 66–68 73–80, 82–90, 92
vector, 69, 72 Singleton type fuzzy rules, 202, 204
Riccati equation, 175, 176, 185 SLAM. See Simultaneous, localization and
Rigid body, v, 9–15, 29–30 mapping (SLAM)
Ring laser gyro (RLG), 28, 29 Sliding window, 67–69, 117
Robust Small perturbation theory, 16–18
adaptive Kalman filter (RAKF), v, vi, SMNSF. See Single, measurement noise scale
96–107 factor (SMNSF)
algorithm, 74, 76 Smoothing, 51, 53
controller, 221 Speed
Robust Kalman filter (RKF), vi, 73–76, 78–92, controller, 144, 149, 154–159, 176–185
99, 102, 103, 107 -over-ground, 38–40
algorithm, vi, 73, 76–78, 81, 91, 92 Spiral mode, 138, 139, 161, 165, 166
Robustness, 78, 97, 107 SPS. See Standard positioning system (SPS)
Roll Stability
angle, 15, 38, 159–161, 165–166, 188 analysis of robust the Kalman filter, 87
mode, 138–139, 161, 164, 166 analysis of the fuzzy controllers, 215–219
motion, 138 augmentation, 95
rate, 160–161 characteristics, 90, 92
Root locus, 141, 143, 145–148, 154, 155, 160, coefficient, 18, 20, 131, 132
161, 164, 167 condition, 59, 137–139, 174
Rudder actuator condition of filter, 59
deflection, 19, 60 of conventional digital filters, 59
criteria, 173
derivatives, 18–20
S of the filter, v, 59, 90
Satellite radio navigation, v, 41–44 function, 174
Scalar measure, 78, 82, 85 limit, 90
Scale factor, 28, 47, 73, 74, 78, 80, 81, 82, 85, of the optimal discrete Kalman filter, 59–60
87, 91 Standard positioning system (SPS), 42
Self-tuning circuit, 64 State
Sensor correction series, 68
bias, 46 equation, 52, 130–135, 160, 188, 189
data fusion, 72 estimates, 47, 54, 96
drift, 46 estimation, v, vi, 47, 51, 55, 56, 60–61, 69,
failures, 102 71, 73, 74, 77, 78, 91, 97, 99, 100,
faults, 45, 46, 47, 71–92, 96, 99, 101–106, 107, 111
110, 111, 114, 116 feedback, 177
fusion, 51, 114, 115 matrices, 162, 163, 176
Sensor/actuator model, 177
failures, vi, 102, 107 prediction, 55
faults, 72, 95–107 space, v, 18, 22, 60, 133, 135, 144, 188
Sequential state estimation, 51 space form, 18, 22, 60, 133, 144
Signal vector, 51, 52, 59–61, 117, 172, 189
constraint method, 212, 213 Static
input, 145 acceleration, 30
processing techniques, 95 pressure, 31, 32
232 Index
Steady state error, 142, 143, 147, 166 True airspeed, 30, 31, 32, 33, 111, 112, 113,
Strapdown system, 26 118, 131
Suboptimality of the filter, 66 True airspeed sensor (TAS), 118–123, 125
Suborbital, 5 TSK fuzzy logic system. See Takagi–Sugeno
Subway control systems, 201 (TSK) fuzzy logic system
Supersonic, 5 Tsukamoto type fuzzy rules, 202
Surface
failure, 99
radar, 33–34 U
System UAV. See Unmanned aerial vehicle (UAV)
dynamics and control, 53 UIKF. See Unknown input Kalman observers
dynamics matrix, 60 (UIKF)
failures, 95 Unit
identification, 51, 223 circle, 60, 88
measurement matrix, 99 Unknown input Kalman observers (UIKF),
model, 64, 177, 188, 206 110
noise, 52, 57, 60, 67, 68, 72, 110, 189 Unmanned aerial vehicle (UAV)
process and measurement noises, 60 dynamics, 47, 51–69, 71–92, 115, 126
states, 67, 111, 177 flight control, vi, vii, 47, 92, 219–221, 224
theory, 59
transfer function, 143, 147
transformation matrix, 68 V
uncertainties, 71 Variance
of the estimation, 77, 101, 104
of Kalman estimation, 62
T of the noise, 81, 105
TACAN, 25 Velocity
Takagi–Sugeno (TSK) fuzzy logic system, 201 controller, 182, 207–210
TAS. See True airspeed sensor (TAS) Vertical velocity, 27, 205
Taylor series, 18, 130, 171 Vision based system, 44
Temperature gradient, 38 VOR/DME, 25
Theoretical error, 74, 75
Three-axis magnetometer, 40
Threshold, 75, 98, 101, 117, 118 W
Thrust actuator, 154 Wash-out filter, 144, 145, 147, 153, 160, 162,
Time 164, 166, 219
constant, 138, 145, 154, 162, 166 Weight matrices, 172, 177, 185
delay, 103 White noise, 54, 67, 190
domain, 59 Wiener filter, 58
to double amplitude, 139
to half amplitude, 137
Total Y
pressure, 31, 32 Yaw
velocity, 27, 38 angle, 15, 17, 161, 164, 171, 187, 195
Trace of damper, 158, 162–165
the covariance matrix, 56
scale matrix, 78, 80, 84, 87, 90
Transfer Z
function, vi, 130–137, 139, 141–143, 145, Zagi UAV, 19–23, 61, 76, 99, 134
147–148, 154, 155, 160, 164, 166 z-domain, 59
matrix, 52, 57, 113 Zero-mean, 52–54, 61, 113, 189
matrix of the system, 52, 57 Ziegler–Nichols method, 149, 156
Transformation matrix, 68 z-transform, 59, 87