M12-T1 — Error Budgets and Covariance Analysis

Learning objectives

After this lesson, you should be able to:

• Compute the sensitivity of the landing miss to each error source from the closed-loop equilibrium of a linear loop.
• Combine independent sources by variance, and combine biased or correlated sources correctly.
• Obtain the stationary variance of a state driven by white noise from the Lyapunov equation.
• Rank the sources by their share of the variance and say which reduction repairs a failing budget.

The intuition

M1-T4: From Mission Requirement to GNC Budget turned the requirement into three numbers and added them as squares. That treated each of navigation, guidance and control as one lump. A real team cannot improve a lump; it improves a sensor, a filter tuning or a gain. The next question is therefore which error, inside navigation, produces how much of the miss.

The answer depends on the loop, not only on the sensor. A tilt error of one degree sounds small. But Hopper's controller reacts to the tilt it believes it has, and a controller that believes it is tilted tries to correct it by thrusting sideways: in the loop of module 4 each radian of pitch error becomes 11 m of landing offset, while a metre of position error becomes exactly one metre. So the tilt sensor needs to be 11 times better, per unit, than the position fix to cost the same. The number that links a source to the miss is its sensitivity, and an error budget is a table of sensitivities multiplied by source sizes.

Three kinds of error source

Constant for a flight. A navigation error that does not change during the few seconds of descent, a target-point error, a wind that blows steadily. The loop settles at an equilibrium that depends on it, and the miss is that equilibrium. Across many flights the value changes at random, with a spread: this is the case a sensitivity handles.

Varying in time. Turbulence is a gust that keeps changing. Its effect on the miss is the stationary spread of the loop driven by random noise, found from a covariance equation rather than an equilibrium.

The same every time. A bias, such as an error in the surveyed coordinates of the pad. It does not average out; it moves the centre of the scatter and so the whole distribution slides off the target.

The picture is the reason to do the analysis. Before the calculation the tilt and rate errors look the most dangerous, because attitude is what makes Hopper unstable. After it, they are a few per cent of the variance: the feedback is so stiff in pitch that small tilt errors are corrected within a second and leave little trace on the landing point. The largest contribution is the plainest one, the position fix, and the second is the wind that the loop (which has no integral state) cannot reject.

The same analysis applies to any vehicle whose loop is linear near its operating point; M10-T2: The Autopilot by Successive Loop Closure would give the aircraft a different table, and the method would not change.

Variance adds across independent sources, but each source is first multiplied by its sensitivity through the loop. The source that matters is the one with the largest product, which is often not the one that looks most dangerous.

Definitions

Error source $e_i$

A quantity with a spread, such as a navigation error in one state, a target offset or a wind acceleration. Zero mean unless stated.

Sensitivity $s_i$

The landing offset, in metres, caused by one unit of the source held constant: $\partial m/\partial e_i$ at the equilibrium of the closed loop.

Closed-loop matrix $A_{cl}$

$A - BK$: the dynamics of the loop with the gain $K$ applied.

Stationary covariance $P$

The covariance of the state when white noise has been acting for long enough that it no longer changes.

Bias

An error with a nonzero mean, the same in every flight. It shifts the scatter instead of widening it.

For a source held constant, the loop settles where its derivative is zero. A source that enters the loop through the input or the measurement does so through a matrix $G_i$:

The sources are combined by variance, with the sensitivities as weights. For correlated sources the covariance matrix of the sources, $\Sigma$, replaces the sum of squares:

For noise that varies in time, the stationary covariance of $\dot x = A_{cl}x + w$ with white noise of intensity $Q$ solves a linear equation:

Derivation: the Lyapunov equation and the effect of correlation

Both results follow from the definition of covariance. The Lyapunov equation first.

1. $$\dot P = \frac{d}{dt}\,\mathbb{E}[xx^{\mathsf T}] = \mathbb{E}[\dot x\,x^{\mathsf T}] + \mathbb{E}[x\,\dot x^{\mathsf T}]$$

   The covariance of a zero-mean state is the average of $xx^{\mathsf T}$; differentiate with the product rule.

2. $$\dot P = A_{cl}P + P A_{cl}^{\mathsf T} + Q$$

   Substitute $\dot x = A_{cl}x + w$. The noise contributes $Q$, its intensity, because $w$ is white.

3. $$A_{cl}P + P A_{cl}^{\mathsf T} + Q = 0$$

   When the loop is stable the covariance stops changing, so its derivative is zero. This is a set of linear equations in the entries of P.

Now the variance of a sum with correlation.

1. $$m = \sum_i s_i e_i \quad\Rightarrow\quad \sigma_m^2 = \mathbb{E}\Bigl[\bigl(\textstyle\sum_i s_ie_i\bigr)^2\Bigr]$$

   The miss is the sensitivity-weighted sum of the sources; its variance is the average of the square.

2. $$= \sum_i s_i^2\,\mathbb{E}[e_i^2] + 2\sum_{i<j}s_is_j\,\mathbb{E}[e_ie_j]$$

   Expand. The cross terms are the covariances of the sources.

3. $$\mathbb{E}[e_ie_j] = \rho_{ij}\sigma_i\sigma_j$$

   The definition of the correlation coefficient. With $\rho = 0$ the cross terms vanish and the rule of module 1 returns.

4. $$\sigma_m^2 = (s_1\sigma_1)^2 + (s_2\sigma_2)^2 + 2\rho\,(s_1\sigma_1)(s_2\sigma_2)$$

   For two sources: a positive correlation between two contributions of the same sign adds to the variance, a negative one removes from it. Two contributions of 0.75 m and 0.20 m with $\rho = 0.8$ add to 0.92 m, not 0.78 m.

A budget is a sensitivity-weighted sum of variances. Bias adds to the mean, not the variance; correlation adds cross terms, which can make the total either larger or smaller than the independent sum.

A worked example

Hopper's lateral loop of module 4 (gain $K = (0.150,\ 0.272,\ 1.658,\ 0.386)$, hover, 0.99 m per axis allowed by the 3 m, 99% requirement) with seven error sources.

1. Requirement. A probability of 99% inside 3 m is 3.035 sigma, so the spread per axis must be 0.989 m or less.
2. Sensitivities. A constant error in the position estimate moves the landing point 1.00 m per metre; velocity 1.81 m per m/s; pitch angle 11.05 m per rad; pitch rate 2.57 m per rad/s; a constant wind acceleration 1.14 m per m/s². With the source sizes 0.75 m, 0.12 m/s, 0.018 rad, 0.006 rad/s and 0.45 m/s², they give 0.750, 0.218, 0.199, 0.015 and 0.512 m.
3. Turbulence. The Lyapunov equation for gust noise of intensity 0.30 gives a stationary position spread of 0.306 m.
4. Total. With a 0.45 m target-point error, the variances add to a spread of 1.099 m against 0.989 m allowed. The probability inside 3 m is 97.6%: the requirement fails.
5. Ranking. Navigation position holds 46.5% of the variance, the gust bias 21.7%, the target point 16.8%. Halving the position error gives 0.887 m and 99.7%; halving both attitude terms barely helps (1.086 m, 97.8%).
6. Bias and correlation. A 0.4 m pad-survey bias on one axis lowers 97.6% to 97.0%. Position and pitch navigation errors with a correlation of 0.8 combine to 0.92 m, not the 0.78 m an independent sum gives.

import math
import gnckit as g

GR, LAM, B = 9.81, 0.7, 19.62
A = [[0, 1, 0, 0], [0, 0, GR, 0], [0, 0, 0, 1], [0, 0, LAM ** 2, 0]]
Bu = [0.0, GR, 0.0, B]                     # input: gimbal angle (rad); it tilts the thrust and turns the body
K = [0.150, 0.272, 1.658, 0.386]           # the lateral LQR of module 4
Acl = [[A[i][j] - Bu[i] * K[j] for j in range(4)] for i in range(4)]

R, P = 3.0, 0.99
k99 = math.sqrt(-2 * math.log(1 - P))
sig_req = R / k99
print(f"requirement: {R:.1f} m with probability {P:.0%} -> {k99:.3f} sigma -> {sig_req:.3f} m per axis")

# Sensitivity of the landing position to each error source: the equilibrium of the closed loop with the source held constant.
def offset(rhs):
    return g.mat_solve(Acl, rhs)[0]

names = ['position', 'velocity', 'pitch angle', 'pitch rate']
units = ['m', 'm/s', 'rad', 'rad/s']
sig_n = [0.75, 0.12, 0.018, 0.006]
sens = []
for i in range(4):
    rhs = [-Bu[r] * K[i] for r in range(4)]        # estimation error e_i enters through the feedback u = -K (x + e)
    sens.append(abs(offset(rhs)))
print("navigation error sources (one sigma) and their sensitivity (m of landing offset per unit of error):")
nav_var = []
for i in range(4):
    c = sens[i] * sig_n[i]
    nav_var.append(c * c)
    print(f"  {names[i]:11s} sigma {sig_n[i]:6.3f} {units[i]:5s} sensitivity {sens[i]:7.3f} m per {units[i]:5s} -> {abs(c):.3f} m")

sig_g = 0.45
wind_gain = abs(offset([0, -1.0, 0, 0]))                # landing offset per unit constant wind acceleration (m per m/s^2)
sig_d = 0.45
print(f"guidance: target-point error {sig_g:.2f} m; wind: {wind_gain:.3f} m per m/s^2, a gust bias of {sig_d:.2f} m/s^2 gives {wind_gain * sig_d:.3f} m")

# Dynamic wind: white acceleration noise of intensity q gives a stationary covariance from the Lyapunov equation.
q = 0.30                                            # (m/s^2)^2 per rad/s, one-sided intensity of the gust acceleration
Q = [[0.0] * 4 for _ in range(4)]; Q[1][1] = q
Pc = g.lyap(Acl, Q)
sig_dyn = math.sqrt(Pc[0][0])
print(f"wind turbulence: stationary position spread from the Lyapunov equation {sig_dyn:.3f} m")

lines = [('navigation: position', nav_var[0]), ('navigation: velocity', nav_var[1]), ('navigation: pitch angle', nav_var[2]), ('navigation: pitch rate', nav_var[3]),
         ('guidance: target point', sig_g ** 2), ('wind: gust bias', (wind_gain * sig_d) ** 2), ('wind: turbulence', sig_dyn ** 2)]
tot = sum(v for _, v in lines)
print("budget (variance share of each source):")
for n_, v in lines:
    print(f"  {n_:25s} {math.sqrt(v):.3f} m  {100 * v / tot:5.1f}% of the variance")
sig_tot = math.sqrt(tot)
p_in = 1 - math.exp(-R * R / (2 * tot))
print(f"total per axis {sig_tot:.3f} m against the {sig_req:.3f} m allowed; probability inside 3 m {p_in:.4f}, outside {1 - p_in:.4f}")

# what-if: halve the largest contribution
vmax = max(lines, key=lambda t: t[1])
tot2 = tot - 0.75 * vmax[1]
print(f"halving the largest source ({vmax[0]}): total {math.sqrt(tot2):.3f} m, probability inside {1 - math.exp(-R * R / (2 * tot2)):.4f}")
# the largest source is not always the cheapest one to fix: what if all the navigation attitude terms are halved
tot3 = tot - 0.75 * (nav_var[2] + nav_var[3])
print(f"halving the two attitude terms: total {math.sqrt(tot3):.3f} m, probability inside {1 - math.exp(-R * R / (2 * tot3)):.4f}")

# a bias adds linearly: a pad survey error of 0.4 m on one axis
bias = 0.4
p_b = 0
# probability inside the circle when the centre is offset by the bias: Monte Carlo free, by numerical integration of the 2-D Gaussian
n = 400
acc = 0.0
h = 2 * R / n
for ix in range(n):
    x = -R + (ix + 0.5) * h
    for iy in range(n):
        y = -R + (iy + 0.5) * h
        if x * x + y * y <= R * R:
            acc += math.exp(-(((x - bias) ** 2 + y ** 2) / (2 * tot))) * h * h / (2 * math.pi * tot)
print(f"a 0.4 m pad-survey bias on one axis: probability inside 3 m {acc:.4f} (without the bias {p_in:.4f})")

# correlated navigation errors: position and pitch errors from the same filter with correlation rho
rho = 0.8
s1, s3 = sens[0] * sig_n[0], sens[2] * sig_n[2]
v_corr = s1 * s1 + s3 * s3 + 2 * rho * s1 * s3
print(f"position and pitch navigation errors with correlation {rho}: their combined spread {math.sqrt(v_corr):.3f} m against {math.sqrt(s1 * s1 + s3 * s3):.3f} m if treated as independent")

requirement: 3.0 m with probability 99% -> 3.035 sigma -> 0.989 m per axis
navigation error sources (one sigma) and their sensitivity (m of landing offset per unit of error):
  position    sigma  0.750 m     sensitivity   1.000 m per m     -> 0.750 m
  velocity    sigma  0.120 m/s   sensitivity   1.813 m per m/s   -> 0.218 m
  pitch angle sigma  0.018 rad   sensitivity  11.053 m per rad   -> 0.199 m
  pitch rate  sigma  0.006 rad/s sensitivity   2.573 m per rad/s -> 0.015 m
guidance: target-point error 0.45 m; wind: 1.138 m per m/s^2, a gust bias of 0.45 m/s^2 gives 0.512 m
wind turbulence: stationary position spread from the Lyapunov equation 0.306 m
budget (variance share of each source):
  navigation: position      0.750 m   46.5% of the variance
  navigation: velocity      0.218 m    3.9% of the variance
  navigation: pitch angle   0.199 m    3.3% of the variance
  navigation: pitch rate    0.015 m    0.0% of the variance
  guidance: target point    0.450 m   16.8% of the variance
  wind: gust bias           0.512 m   21.7% of the variance
  wind: turbulence          0.306 m    7.8% of the variance
total per axis 1.099 m against the 0.989 m allowed; probability inside 3 m 0.9759, outside 0.0241
halving the largest source (navigation: position): total 0.887 m, probability inside 0.9967
halving the two attitude terms: total 1.086 m, probability inside 0.9780
a 0.4 m pad-survey bias on one axis: probability inside 3 m 0.9698 (without the bias 0.9759)
position and pitch navigation errors with correlation 0.8: their combined spread 0.917 m against 0.776 m if treated as independent

In robotics

An analysis of this kind sits behind every serious budget.

• Allocation by effort. Because the variance goes as the square of size, the cheapest way to repair a failing budget is almost always to reduce the largest term. In the example one change, a better position fix, repairs the whole budget; no change in tilt sensing could.
• Sensors are bought by their sensitivity. A gyro grade, a beacon geometry or an altimeter is worth buying when the sensitivity of the miss to it is large. In Hopper's loop the answer for pitch is no, because the loop is stiff; for an aircraft landing at a long range it might be yes.
• Finding what the loop cannot reject. The second largest term, the gust bias, is a constant disturbance that the loop has no integral to reject. The fix is the integral state of M4-T6: Tracking, Integral Action and Feedforward, which removes the term altogether rather than shrinking it.
• Covariance analysis in the filter. The navigation errors in the table are the output of the filter of M7-T2: Coupling GNSS and INS: Loose, Tight and Deep; their correlation, which the filter's covariance matrix reports, is the input to the cross terms above.
• Where it stops. The analysis assumes a linear loop and Gaussian sources. Saturation, delays and a beacon fault break both assumptions, and for those the next lesson, M12-T2: Monte Carlo Dispersion Analysis, replaces the algebra with a simulated campaign.

The covariance propagation used here is the same one that appears in Linear Algebra M9-T4: Covariance Propagation and Probability M4-T4: Propagating Uncertainty Through Transformations.

Common misconceptions

Key takeaways

1. The miss caused by a constant error source is its size times a sensitivity, which is an equilibrium of the closed loop.
2. Independent sources combine by variance; correlated sources need the cross terms; a bias moves the centre, not the width.
3. Time-varying noise gives a stationary covariance from the Lyapunov equation of the closed loop.
4. In the example the sum is 1.10 m per axis against 0.99 m allowed: 97.6% inside 3 m instead of 99%.
5. The ranking, not the total, says what to fix: halving the 0.75 m position error repairs the budget (99.7%), halving the attitude terms does not.
6. The analysis assumes a linear loop and Gaussian sources, which is why a simulated campaign follows.

References and further reading

Each source is marked free or paid. Free to read does not mean free to reuse: the course resources page lists every source with its licence.

• A. Gelb (ed.), Applied Optimal Estimation, MIT Press, 1974, the chapters on covariance analysis and error budgets. Paid.
• B. D. O. Anderson and J. B. Moore, Optimal Filtering, Prentice-Hall, 1979, for the Lyapunov equation and stationary covariance; out of print, available from libraries. Free to borrow.
• P. S. Maybeck, Stochastic Models, Estimation, and Control, vol. 1, Academic Press, 1979, an introduction to covariance propagation. Paid (reprinted by Navtech).
• MIT 16.30 Feedback Control Systems (OpenCourseWare) , for the closed-loop equilibrium and the steady-state error. Free; licence CC BY-NC-SA 4.0.

Knowledge check

Three questions. Write an answer, then reveal the explanation. Your answers stay in your browser.

What’s next

Related content

Nothing below is required to finish this lesson.

• Lab
  M12-L1: How Many Runs Are Enough?

  A requirement says 99% of landings within 3 m. Build the error budget, run a dispersion campaign, and state with what confidence the requirement is met.

Other courses on this site treat neighbouring topics from their own angle:

• Probability M1-T3: Simulating Randomness: Sampling and Monte Carlo
• Probability M4-T4: Propagating Uncertainty Through Transformations
• Linear Algebra M9-T4: Covariance Propagation
• UAV M8-T9: Fault detection, isolation and recovery
• UAV M8-T7: FMEA
• UAV M8-T3: Verification and validation
• UAV M8-T4: SIL, HIL and ground testing
• UAV M6-T14: Software-in-the-Loop simulation
• UAV M7-T7: Jitter and deadline monitoring
• UAV M8-T1: System requirements engineering
• UAV M8-T5: Flight-test planning
• Navigation M10-T2: Engineering Trade-offs and Safety Cases
• Robot Learning M11-T2: Seeds, Variance and Statistics

This course leaves the following to the pages that already teach them:

• Requirements engineering, FMEA and fault-tree analysis as methods
• SIL and HIL test benches, real-time scheduling and time synchronization
• Flight-test planning and incident investigation