ACS6124 Multisensor and Decision Systems
Part I – Multisensor Systems
Assessment
Lecturer: Yuanbo Nie
email: y.nie@sheffield.ac.uk
2022/2023
This is the first of two assignments for the module ACS6124, with this one assessing the learning
outcomes of the Multisensor Systems component of the module.
Assignment Code: ACS6124-001
Assignment weighting: 50% of the module assessment marks.
Assignment released: 28 Feb 2023 (Tuesday, Week 4, Semester 2)
Assignment due: 23:59, 24 April 2023 (Monday, Week 9, Semester 2)
Late submission penalty will be applied
Assignment format: A report of 15 pages maximum (excluding the Appendix)
Accompanied Matlab code to be submitted and evaluated
Submission mode: Report and code must be submitted electronically via Blackboard
1 Assignment Processes
Report format details
The report should be in a pdf form using a margin (top, bottom, left and right) of no less than
20mm, and text should be size 12 point (minimum 11pt) with 1.5 line spacing.
Submission mode details
This report and the code are to be submitted via Turnitin using the Assessment tab on Blackboard
by 23:59, Monday 24th April 2023 (Week 9 of Semester 2). Please note that you will be allowed
only a single submission and so the first submission made will be the one assessed.
Penalties for Late Submission
Work that is submitted after this deadline (without medical or other similar documented evidence
unless agreed) will incur a penalty for late submission. The usual late submission penalty of 5%
reduction in the mark for every working day (or part thereof) that the assignment is late and a mark
of 0 for submission more than 5 days late. For more information http://www.shef.ac.uk/ssid/
assessment/grades-results/submission-marking
1
Unfair Means
This is an individual assignment. You should not discuss the assignment with other students or
work together with other students in its completion. The assignment must be wholly your own
work. References must be provided to any other work that is used as part of the assignment. Any
suspicion of the use of unfair means will be investigated and may lead to penalties. For more
information see: http://www.shef.ac.uk/ssid/unfair-means
Extenuating Circumstances
If you have any medical or special circumstances that you believe may affect your performance on
the assignment then you should raise these with the Module Leader at the earliest opportunity.
You will also need to submit an extenuating circumstances form. More information at: http:
//www.shef.ac.uk/ssid/forms/circs
Support
Please email me if you have any questions about the assignment.
Feedback
Written feedback will be provided on Blackboard within 15 working days, in line with Department
guidelines.
Assignment technical requirements
Write a report fulfilling the requirements of the Task 1 to 3 of the Assignment. The report should
follow the specified guidelines and meet the objectives provided for each task. Task 4 will be graded
through the evaluation of the submitted code.
Link to Labs
The workload of this assignment has been designed based on the assumption that the lab exercises
have all been successfully completed.
• Tasks requested in this assignment are extensions of lab activities to a real-world engineering
problem. For example
– In Lab B, scripts for running extended Kalman Filter and iterated extended Kalman
Filter have been coded up for scalar nonlinear dynamical systems. In this assignment,
they simply need to be modified to work with multiple states, inputs and outputs.
– In Lab C, script for fault detection using the CUSUM method has been created. It can
considered as an option to achieve various fault detection tasks in this assignment with
appropriate modifications.
• Critical thinking developed in lab activities can be helpful in addressing the challenges faced
in this assignment.
• Completion of this assignment requires a similar level of Matlab skills as in lab activities.
2
2 Assignment: Integrated Navigation for Aircraft
This project concerns the accurate estimation of aircraft position, airspeed body components and
attitude with an integrated system consisting of inertial measurement units (IMUs), the global
positioning system (GPS), and air-data systems. A schematic drawing of the system architecture is
shown in Figure 1.

State
Estimation
and Fault
Detection
Algorithm
Inertial Measurement Unit
Multi-Antenna GPS
Air-data System
Acceleration
Rotational Rates
Position & Altitude
Ground Speed
Components
Attitude
True Airspeed
Angle of Attack
Side-slip Angle
Position
Altitude
Airspeed Body
Components
Roll Angle
Pitch Angle
Yaw Angle
Windspeed
Components
Sensor Biases
& Faults
Figure 1: Overview of the Sensor Fusion Architecture.
The IMU uses rate gyros for aircraft attitude and heading angle estimation. However, these sensors
have a major drawback of long-term instability caused by uncertainties (biases, drifts, scale factors,
misalignments, nonlinearities and noises). Attitude and heading estimation with only rate gyros
yields time accumulated errors.
The integration of rate gyros with other long-term stable sensors has the capacity of calibrating
gyros in flight. The integrated navigation system has the following advantages:
• long-term stable measurements,
• sensor system uncertainties can be estimated and removed,
• high data rate for flight control and display,
• sensor fault detection and re-configuration capability,
• flexible integration structures.
The measurements available are as follows: the IMU will measure the acceleration (Ax, Ay, Az in
m/s2) and rotational rate (roll rate p, pitch rate q, yaw rate r in rad/s) of the aircraft, the GPS
receiver will provide aircraft position (xE , yE , zE in m), ground speed earth components (uE , vE ,
wE in m/s) and attitude (roll angle ϕ, pitch angle θ, yaw angle ψ in rad), and the air-data sensors
will offer processed true airspeed (VTAS in m/s), angle of attack (α in rad), and side slip angle (β
in rad).
3
2.1 Dynamics Model
For the Kalman filter design, we would like to express the system dynamics in the following general
format
x˙(t) = f(x(t), c(t), t) + g(x(t), ω(t)), (1)
d(t) = h(x(t), t), (2)
dm(tk) = d(tk) + ν(tk), k = 1, 2, . . . (3)
with kinematic equation (1), observation equation (2) and measurement model (3). x is the state
variable, c is the input variable, d is the output variable with dm the measured outputs, ω is the
system noise and ν is the output measurement noise. Note that the notation used here is different
from the lecture slides to avoid confusions with problem-specific variable names.
2.1.1 Kinematic Equations
For aircraft trajectory in a short range, the simplified kinematic equations using flat non-rotating
earth assumption are
x˙E =(u cos θ + (v sinϕ+ w cosϕ) sin θ) cosψ − (v cosϕ− w sinϕ) sinψ + VwxE ,
y˙E =(u cos θ + (v sinϕ+ w cosϕ) sin θ) sinψ + (v cosϕ− w sinϕ) cosψ + VwyE ,
z˙E =− u sin θ + (v sinϕ+ w cosϕ) cos θ + VwzE ,
u˙ =Ax − g sin θ + rv − qw,
v˙ =Ay + g cos θ sinϕ+ pw − ru,
w˙ =Az + g cos θ cosϕ+ qu− pv,
ϕ˙ =p+ q sinϕ tan θ + r cosϕ tan θ,
θ˙ =q cosϕ− r sinϕ,
ψ˙ =q
sinϕ
cos θ
+ r
cosϕ
cos θ
,
with VwxE , VwyE , VwzE wind components in the earth reference frame. They are assumed to be
constant, i.e. V˙wxE = 0, V˙wyE = 0 and V˙wzE = 0 but they are also unknowns hence need to be
estimated by the algorithm. u, v and w are airspeed components in the aircraft’s body reference
frame which also can not be directly measured.
In classical aircraft dynamic models, the acceleration and angular rates are additionally expressed
as functions of flight control inputs, aerodynamic coefficients and other state and output variables.
To avoid working with such a model with complicated aerodynamic coefficients, for this assignment
we will use a multi-sensor integration structure where the IMU measurements directly serve as
the system input and the GPS and air-data measurements serve as the outputs. Note that this
is sufficient to estimate all the state variables. This sensor fusion scheme is often used in flight
tests where the aerodynamic coefficients are not yet known in good precision, as well as in fault
tolerant flight control designs where structure failures in flight could lead to changes in aerodynamic
characteristics. Therefore, we have state vector
x = [xE , yE , zE , u, v, w, ϕ, θ, ψ, VwxE , VwyE , VwzE ]
⊤,
input vector c = [Ax, Ay, Az, p, q, r]
⊤, and output vector
d = [xGPS , yGPS , zGPS , uGPS , vGPS , wGPS , ϕGPS , θGPS , ψGPS , VTAS , α, β]
⊤.
4
2.1.2 Observation Equations
A GPS receiver with multiple antennas offers 9 observables. They are the position coordinates,
ground speed components, and attitude angles:
xGPS =xE ,
yGPS =yE ,
zGPS =zE ,
uGPS =uE = (u cos θ + (v sinϕ+ w cosϕ) sin θ) cosψ − (v cosϕ− w sinϕ) sinψ + VwxE ,
vGPS =vE = (u cos θ + (v sinϕ+ w cosϕ) sin θ) sinψ + (v cosϕ− w sinϕ) cosψ + VwyE ,
wGPS =wE = −u sin θ + (v sinϕ+ w cosϕ) cos θ + VwzE ,
ϕGPS =ϕ,
θGPS =θ,
ψGPS =ψ,
Airdata sensors give three variables: the true airspeed vTAS , angle of attack α and side slip angle
β:
VTAS =
√
u2 + v2 + w2,
α =tan−1
w
u
,
β =tan−1
v√
u2 + w2
,
Together they form the observation model (2).
2.1.3 Measurement Equations
In the kinematic equation (1) and measurement model (3), we see the appearance of system noise
ω and output measurement noise ν. For this assignment, we will model the process noises related
to input measurements, i.e. for the input measurements from the IMU we have
Axm = Ax + ωAx , Aym = Ay + ωAy , Azm = Az + ωAz ,
pm = p+ ωp, qm = q + ωq, rm = r + ωr,
with ω = [ωAx , ωAy , ωAz , ωp, ωq, ωr]
⊤ and
E{ω(tk)} = 0, E{ω(ti)ω⊤(tj)} = Qδij , Q = diag(σ2Ax , σ2Ay , σ2Az , σ2p, σ2q , σ2r).
For outputs from the GPS and air-data sensors, consider
xGPSm = xGPS + νx, yGPSm = yGPS + νy, zGPSm = zGPS + νz,
uGPSm = uGPS + νu, vGPSm = vGPS + νv, wGPSm = wGPS + νw,
ϕGPSm = ϕGPS + νϕ, θGPSm = θGPS + νθ, ψGPSm = ψGPS + νψ,
VTASm = VTAS + νVTAS , αm = α+ να, βm = β + νβ ,
with measurement noises ν = [νx, νy, νz, νu, νv, νw, νϕ, νθ, νψ, νVTAS , να, νβ ]
⊤ and
E{ν(tk)} = 0, E{ν(ti)ν⊤(tj)} = Rδij ,
R = diag(σ2xE , σ
2
yE , σ
2
zE , σ
2
u, σ
2
v , σ
2
w, σ
2
ϕ, σ
2
θ , σ
2
ψ, σ
2
VTAS , σ
2
α, σ
2
β).
5
Therefore, we have and input measurement vector cm = [Axm , Aym , Azm , pm, qm, rm]
⊤ and output
measurement vector
dm = [xGPSm , yGPSm , zGPSm , uGPSm , vGPSm , wGPSm , ϕGPSm , θGPSm , ψGPSm , VTASm , αm, βm]
⊤.
The standard deviations (σ) of the input and output noises are assumed to be constant and known
in nominal flight: 0.01 m/s−2 for accelerations (i.e. for σAx , σAy , σAz ), 0.01 deg/s for angular rates
(i.e. for σp, σq, σr), 5 meters for GPS position measurements (i.e. for σxE , σyE ), 10 meters for GPS
altitude measurements (i.e. for σzE ), 0.1 m/s for velocity measurements (i.e. for σu, σv, σw, σVTAS ),
and 0.1 degrees for angle measurements (i.e. for σϕ, σθ, σψ, σα, σβ). When using these numbers,
note that for dynamics equations, all angles should be in radians and all angular rates should be in
rad/s. Also, the input and output noises are assumed to be uncorrelated, i.e. E{ω(ti)ν⊤(tj)} = 0.
2.2 Assignment Tasks
Task 1: Kalman Filter Design
You will need to design an integrated GPS/IMU/Air-data navigation system by constructing a nav-
igation model with 12 states (position, velocity, attitude and wind components), 6 inputs and 12
measurements, and extend your scalar extended Kalman filter (EKF) or iterated extended Kalman
filter (IEKF) algorithm from the lab to work with vectors. The choice between EKF and IEKF is
at your own discretion.
For this task you will work with dataTask1.mat which contains input measurements c k, the outputs
d k and the time variables t and dt. c k and d k are matrices with rows the recordings of different
time instances and the columns the variables are in the same ordering of cm and dm presented earlier.
You will realize that in comparison to your lab exercise, the expected value for the initial condition
xˆ0|0 = E{x(t0)} and the error covariance P0|0 are not explicitly given. You could use engineering
judgement to configure this. For example, for state variables that can be directly measured (xE , yE ,
zE , ϕ, θ, ψ), simply assume the initial value to be the same as the noisy measurements at the first
time instance.
• Task 1.1: For the use of Kalman filter, you will need to derive the Jacobian matrices for our
nonlinear system. Present the expressions for your derived Jacobian matrices F = ∂f∂x , G =
∂g
∂ω
and H = ∂h∂x in the report with f and g from (1) and h from (2). Hint: substitute IMU
measurement equations into the kinematic equations, then group different terms appropriately
into f(x(t), c(t), t) and g(x(t), ω(t)). Equations for the kinematic model and observation model
are can be found in usefulEquations.m, and you can use Matlab’s Symbolic Math Toolbox
to help you compute the required matrices. [10%]
• Task 1.2: In clear figures, plot the trajectories of the estimated 12 states and comment on your
findings. [10%]
Task 2: Removing Sensor Biases
All sensors will have errors even in nominal operations, typically in the form of scale factors, biases
and noises. In Task 1 we have observed the handling of noisy measurements with the Kalman filter.
In this part of the assignment, we will extend the algorithm to deal with sensor biases.
In particular, the IMU is known to develop gyroscopic bias over time. dataTask2.mat contains the
sensor measurement for the very same flight trajectory, with the addition of IMU sensor biases.
• Task 2.1: Re-run Task 1 using dataTask2.mat and compare the trajectories of the estimated
states with the ones obtained from the previous task. Show the comparison in a clear figure
(i.e. trajectories from both tasks shown together in the same figure) for the variables with
6
the most pronounced differences (around 5 variables would be a good selection). Explain the
reasons behind the degradation in the estimation performance. [10%]
In Kalman filter design, a common practice to handle unknown biases is to simply include them as
variables to be estimated. Let’s consider the case of constant bias where we have
Axm = Ax + bAx + ωAx , Aym = Ay + bAy + ωAy , Azm = Az + bAz + ωAz ,
pm = p+ bp + ωp, qm = q + bq + ωq, rm = r + br + ωr,
b˙Ax = 0, b˙Ay = 0, b˙Az = 0,
b˙p = 0, b˙q = 0, b˙r = 0.
• Task 2.2: Implement the integrated navigation scheme using the updated navigation model
with 18 states (position, velocity, attitude, wind components and sensor biases), 6 inputs
and 12 measurements. Verify that the influences of unknown biases have been successfully
mitigated. In a clear figure, present your estimation of the bias terms through the Kalman
filter scheme. [10%]
Task 3: Sensor Fault Detection and Diagnose
The integrated navigation framework with Kalman filter can also be used for sensor fault detection
and diagnosis (FDD). We will consider additive faults which can be viewed as a generalization of
the bias case, but with bias terms no longer constant.
• Task 3.1: Run the state estimation with the updated navigation model using data from
dataTask3 1.mat and identify the instances when sensor faults occur (there may be more
than 1 instance). In a clear figure, present the history of the innovation of the EKF/IEKF.
Comment on your observation and provide some rationales behind the use of Kalman filter
innovation/residual for fault detection. [10%]
Kalman filter-based FDD often model the faults as a random walk process. For our problem, we
will have
bAx(tk+1) = bAx(tk) + ωbAx (tk), bAy (tk+1) = bAy (tk) + ωbAy (tk), bAz (tk+1) = bAz (tk) + ωbAz (tk),
bp(tk+1) = bp(tk) + ωbp(tk), bq(tk+1) = bq(tk) + ωbq (tk), br(tk+1) = br(tk) + ωbr (tk).
• Task 3.2: Update the navigation model using 6 additional noise terms with the faults (bias)
modelled as in (11). Now you should have
ω = [ωAx , ωAy , ωAz , ωp, ωq, ωr, ωbAx , ωbAy , ωbAz , ωbp , ωbq , ωbr ]
⊤.
The standard deviation of the additional noise terms needs to be set to relatively high values
to capture the fault signals. Iteratively improve your design and verify that the influences of
faults have been successfully mitigated. Check the history of the innovation of the EKF/IEKF
again and compare it with what you have observed in Task 3.1. Briefly explain the differences. [10%]
• Task 3.3: Consider dataTask2.mat to be the measurements from nominal operations and
dataTask3 1.mat to be the one with various sensor faults, under the same noise and bias
magnitudes. Develop a fault detection scheme to identify the time instances where faults
occur for measurements in Ax, Ay, Az, p, q and r respectively. Justify your design choice of
such a scheme in a paragraph. Also, present your estimation of the fault signals through the
Kalman filter scheme in a clear figure, together with the time instances where the detection
scheme raises alarms. Hint: when generating the figure, scale the axes to clearly show the
features of the fault, rather than the Kalman filter convergence process. [10%]
7
Beyond input measurements, malfunctioning could also occur in sensors for measuring the outputs.
For example, the pitot tube for airspeed measurement and the angle of attack sensor could freeze
and provide erroneous readings due to ice (e.g. in the cases of Air France flight 447 and XL Airways
Germany flight 888T). A similar effect could occur in the event of a data-injection attack, where
false signals replace the sensor measurements.
• Task 3.4: dataTask3 2.mat contains the data with the angle of attack sensor stuck mid-flight.
Enhance your algorithm design to be able to detect such events, as well as to be able to
continue to provide reliable state estimation afterwards. To simply this task for you, you can
assume it is known beforehand that the fault come from the angle of attack sensor. You have
the full freedom to design such an FDD scheme but you should not further increase the number
of state variables in the Kalman filter. Justify your design choice in a paragraph and show
evidences of successful fault detection and mitigation. [10%]
Task 4: Estimation and FDD Performance Competition
In this part of the assignment, you will submit your design of the state estimation and FDD algorithm
and its performance will be measured against the work of your fellow peers. Use the Blackboard
template SubmissionTemplate.m and follow the instructions inside to prepare the script. You may
use dataTask4.mat to further test your design, which contains a different flight tape with different
occurrences of failures.
WARNING: It is your responsibility to ensure the proper functioning of your script, as run-time
errors can lead to a mark of 0% for this task. You can use testYourScript.m to test your script on
different flight tapes, to help you identify any syntax error or output variable dimension mismatches.
After submission, the performance will be evaluated as follows based on another flight trajectory:
• Estimation Performance The estimation performance will be evaluated based on the weighted
mean squared error (WMSE) between the estimated states and the true states. The assigned
weights are shown as below:
States xE yE zE u v w ϕ θ ψ VwxE VwyE VwzE
Weights 0.2 0.2 0.1 10 10 10 5000 5000 5000 10 10 10
[6%]
• Detection Performance The detection performance evaluation contains two parts. First, the
estimation of the input measurement bias/fault signals will be evaluated based on the weighted
mean squared error. The weights assigned are shown as below:
States bAx bAy bAz bp bq br
Weights 1 1 1 50 50 50
Next, the time instances for the first detection of faults in measurements are evaluated, in-
dependently for each input variable of Ax, Ay, Az, p, q and r, and output variable α. Note
that input measurements will contain some constant biases in nominal operations hence these
biases should not lead to fault detection. Triggering of alarm before the occurrence of faults
or after the end of faults will be marked as 0% for this particular detection task. [14%]
8
3 Marking Criteria
This assignment is marked out of 100% and contributes to 50% of the overall module assessment.
Among this 100%, 20% is associated with the performance of the submitted Matlab script as part
of Task 4, and 80% is from the marking of the report.
The full range (0-80%) will be used when marking the report. Most part of this assignment require
you to present the results in clear figures. Therefore, the marking will have a strong emphasis on
the clarity and well thoughtfulness of these figures. Please aim to generate clear figures and write
concise discussions, keep in mind that the 15 pages allowance is a limit rather than a target. Con-
tents beyond a maximum of 15 pages (include all texts and figures) will be disregarded.
Additionally, penalties up to 5% for EACH violation of the following requirements will be applied:
• minimum 11pt font size for the main text
• minimum 10pt font size texts in the figures
• page margin (top, bottom, left and right) no less than 20mm, strictly applied for both texts
and figures
• proper figure references (e.g. not appears to be ??)
• proper figure captions
WARNING: You may lose up to 25% of your Part I assignment grade due to formatting issues.
DOUBLE CHECK before you submit. Note that most LATEX templates do not allow arbitrary font
sizes (e.g. 11.5pt).
Marking Criterion / Comments Marks
Task 1.1
• Correctly derive the expression for matrix F / 3
• Correctly derive the expression for matrix G / 4
• Correctly derive the expression for matrix H / 3
Task 1.2
• Clear figures with the trajectory of the estimated states / 4
• The estimation of state trajectories in the figures are correct / 6
Task 2.1
• A clear figure showing key differences / 5
• Explanations of the degradation in estimation performance / 5
Task 2.2
• A clear figure with the estimation of bias terms / 4
• The estimation of bias terms in the figures are corrects / 6
Task 3.1
• A figure showing the history of KF innovation / 5
• Rationales behind the use of KF innovation for fault detection / 5
Task 3.2
• Show evidences that influences of faults have been successfully mitigated / 4
• Comparison of innovations of KF with brief explanations / 6
Task 3.3
• Justification for the choice of fault detection scheme / 4
• A clear figure with the estimation of fault signals and detection times / 6
Task 3.4
• Justification for the choice of fault detection scheme / 4
• A clear figure showing evidences of successful fault detection and mitigation / 6
9