M12-T3 — A Stability-Margin Policy

Learning objectives

After this lesson, you should be able to:

• State a margin policy with nominal and dispersed requirements for gain, phase and disk margin.
• Compute the lower and upper gain margins of an unstable plant and the delay margin of a loop.
• Verify margins over a grid of dispersed parameters and identify the corner that governs.
• Confirm frequency-domain margins in the time domain, and turn a violated margin into a requirement on a component.

The intuition

Hopper's pitch loop is the reference of module 4: a PD controller on an unstable axis, natural frequency 4 rad/s, damping 0.7. By itself the loop has a phase margin of 55.5° with the 40 rad/s gimbal servo. Is that enough? The answer depends on how far the real loop can be from the model. The servo may be slower in cold weather, the flight computer adds delay, the mass changes the authority by 40%, and the position of the centre of mass changes the divergence rate. A margin that is comfortable at the nominal point may be gone at the corner of that range.

What a policy contains

Requirements at two levels. A nominal requirement (here 45° of phase margin, 6 dB of gain margin each way, a disk margin of 0.5) applies to the design as modelled. A smaller dispersed requirement (30°, 4.5 dB, 0.35) applies at every combination in the stated envelope, because the model is not the vehicle. The two numbers are commonly paired in the aircraft specifications that inspired them; the governing document of a real programme must be consulted, and the values here are illustrative.

An envelope. The list of parameters to vary and their ranges: authority 11.77 to 19.62 (mass), servo 25 to 60 rad/s, delay 10 to 60 ms, divergence rate 0.6 to 0.9 per second. The policy says that the requirement holds everywhere in the envelope, not just at corners that were tried first.

Where the loop is broken. Margins depend on the point at which the loop is cut. The usual choices are at the actuator input (to test the effect of the actuator and the controller) and at the sensor output; a policy requires both, since errors enter in both places.

A confirmation in the time domain. Margins come from a linear model. A simulation with the nonlinear actuator (rate and position limits) checks that the loop behaves the way the margins say.

Unstable plants have two gain margins

For a stable plant a loop loses stability when the gain is raised. Hopper's pitch axis is unstable, and the controller has to supply stiffness just to hold it up, so it also loses stability when the gain falls too low: below 0.0297 times the design gain (30.5 dB) the feedback can no longer overcome the unstable term. That lower margin is large here, and the upper margin is set by the servo and the delay. Both are written into the policy, and both must be checked.

A margin policy is a statement of how much is required, over what range of conditions, at which loop breaking points. Each of those three can be wrong while the others are right.

Definitions

Gain margin

The factor by which the loop gain can change before instability, in dB. An unstable plant has a lower and an upper one.

Phase margin

The extra phase lag, at the gain crossover frequency $\omega_c$, that would make the loop unstable.

Delay margin

The pure delay that has the same effect as the phase margin: $\tau_m = \mathrm{PM}/\omega_c$, with PM in radians.

Disk margin $\alpha$

The radius of the largest disk of simultaneous gain and phase changes that the loop tolerates, as in the robustness lesson of module 4.

Envelope

The set of parameter values over which the requirements must hold.

Nominal and dispersed requirement

The margins required at the design point and at every point of the envelope.

For the pitch loop with controller $K(s) = K_P + K_Ds$, plant $b/(s^2-\lambda^2)$, a servo $\omega_a/(s+\omega_a)$ and a delay $\tau$, the loop is:

At the gain crossover frequency $\omega_c$ (where $|L| = 1$):

For the lower gain margin, consider the closed-loop characteristic polynomial with the gain multiplied by $k$:

Derivation: phase margin, delay margin, and the lower gain margin

The three results come from reading the loop on the imaginary axis.

1. $$L(j\omega) = (K_P + jK_D\omega)\cdot\frac{b}{-\omega^2-\lambda^2}\cdot\frac{\omega_a}{j\omega+\omega_a}\cdot e^{-j\omega\tau}$$

   Substitute $s = j\omega$. The plant factor is real and negative: $-b/(\omega^2+\lambda^2)$.

2. $$\angle L = -180^\circ + \arctan\frac{K_D\omega}{K_P} - \arctan\frac{\omega}{\omega_a} - \omega\tau$$

   The phases add. The negative real plant contributes −180°; the PD zero adds a lead; the servo pole and the delay add lags that grow with frequency.

3. $$\mathrm{PM} = 180^\circ + \angle L(j\omega_c) = \arctan\frac{K_D\omega_c}{K_P} - \arctan\frac{\omega_c}{\omega_a} - \omega_c\tau$$

   The phase margin is the distance from −180° at the crossover. With $\omega_c = 6.07$, $\omega_a = 40$ and no delay: $64.1^\circ - 8.6^\circ = 55.5^\circ$.

4. $$\tau_m = \frac{\mathrm{PM}}{\omega_c} = \frac{0.969\ \mathrm{rad}}{6.07\ \mathrm{rad/s}} = 0.160\ \mathrm{s}$$

   A delay adds a phase lag of $\omega_c\tau$ at the crossover without changing the magnitude; the delay that uses up the whole margin is the delay margin.

The lower gain margin comes from the Routh test of the characteristic polynomial. For a second-order polynomial, stability needs all coefficients positive:

1. $$s^2 + (kbK_D)\,s + (kbK_P - \lambda^2)$$

   With the loop gain multiplied by $k$, the closed loop of the plant $\ddot\theta = \lambda^2\theta + b\delta$ with $\delta = -k(K_P\theta + K_D\dot\theta)$ has this polynomial.

2. $$kbK_D > 0 \quad\text{and}\quad kbK_P > \lambda^2$$

   Both coefficients must be positive: the first holds for any positive gain, the second sets the lower limit.

3. $$k_{\min} = \frac{\lambda^2}{bK_P} = \frac{0.49}{19.62\times0.8405} = 0.0297 \;\Rightarrow\; 20\log_{10}\frac{1}{k_{\min}} = 30.5\ \mathrm{dB}$$

   The gain may be reduced to 3% of its design value. At liftoff mass; at burnout ($b = 11.77$) the margin is 21.7 dB at $\lambda = 0.9$.

On the imaginary axis the plant is a fixed −180°, so the phase margin of the pitch loop is the PD lead minus the servo lag minus the delay. The unstable plant also gives a lower gain margin, set by the stiffness needed to hold it up.

A worked example

The reference pitch loop ($K_P = 0.8405$, $K_D = 0.2854$, designed for $b = 19.62$) and the policy: nominal 45°, 6 dB, disk 0.5; dispersed 30°, 4.5 dB, 0.35.

1. The reference loop. With an ideal actuator: crossover 6.13 rad/s, phase margin 64.3°, disk margin 1.12. With the 40 rad/s servo: 6.07 rad/s, 55.5°, disk margin 1.00, lower gain margin 30.5 dB, no upper limit (the loop has no −180° crossing without delay). Delay margin 0.160 s.
2. Nominal check. With the 10 ms delay of the flight computer: 52.0°, 30.5 dB down, 24.8 dB up, disk margin 0.95: the nominal requirement is met with room.
3. The envelope. 81 combinations of authority (11.77, 15.0, 19.62), servo (25, 40, 60 rad/s), delay (10, 30, 60 ms) and divergence rate (0.6, 0.7, 0.9). Worst phase margin 29.7° and worst upper gain margin 9.1 dB, both at full authority with the 25 rad/s servo and 60 ms; worst lower gain margin 21.7 dB at burnout authority and λ = 0.9; worst disk margin 0.52.
4. Violations. Three combinations fail the dispersed requirement, all with the slowest servo and the longest delay at full authority, with divergence rates 0.6, 0.7 and 0.9: phase margin 29.7° to 29.8°. Nothing else in the envelope fails.
5. Time domain. With the rate limit of 3 rad/s and the position limit of 0.15 rad, a 0.05 rad tilt recovers to within 5% in 0.68 s at the nominal corner and 0.97 s at the worst corner. From 1.0 rad the times are 1.06 s and 1.83 s. Raising the gain at the worst corner (no limits) by 8.2 dB leaves the loop stable and by 9.9 dB makes it diverge, against a computed upper margin of 9.1 dB.
6. The derived requirement. The one failing corner is removed by a requirement on a component, not on the controller: the gimbal servo must be faster than 25 rad/s, or the total loop delay must be below about 59 ms. A servo of 40 rad/s with up to 60 ms gives a margin of 34.6°, and 60 rad/s gives 37.5°.

import math, cmath

LAM2 = 0.7 ** 2                       # divergence rate squared of the unstable pitch axis
B0 = 19.62                            # pitch authority at liftoff, m g l / J
KP, KD = 0.8405, 0.2854               # the reference PD of module 4, designed for B0
WS = [10 ** (-1.5 + 4.0 * i / 3999) for i in range(4000)]

def L(w, b=B0, servo=40.0, delay=0.0, lam2=LAM2):
    s = 1j * w
    return (KP + KD * s) * b / (s * s - lam2) * servo / (s + servo) * cmath.exp(-s * delay)

def margins(**kw):
    """Gain crossover, phase margin, lower and upper gain margin (dB) and disk margin of the loop broken at the gimbal."""
    ph = []; un = 0.0; last = None
    for w in WS:
        a = math.degrees(cmath.phase(L(w, **kw)))
        if last is None: un = a
        else:
            d = a - last
            while d > 180: d -= 360
            while d < -180: d += 360
            un += d
        last = a; ph.append(un)
    wc = pm = None; low = up = None
    for i in range(1, len(WS)):
        m0, m1 = abs(L(WS[i - 1], **kw)), abs(L(WS[i], **kw))
        if m0 > 1 >= m1 and wc is None:
            lo, hi = WS[i - 1], WS[i]
            for _ in range(60):
                mid = 0.5 * (lo + hi)
                if abs(L(mid, **kw)) > 1: lo = mid
                else: hi = mid
            wc = 0.5 * (lo + hi)
            f = (math.log(wc) - math.log(WS[i - 1])) / (math.log(WS[i]) - math.log(WS[i - 1]))
            pm = 180 + ph[i - 1] + f * (ph[i] - ph[i - 1])
        # phase crossings of -180 (the phase of an unstable plant starts at -180 and rises)
        if (ph[i - 1] + 180) * (ph[i] + 180) < 0:
            mag = abs(L(WS[i], **kw))
            if mag < 1: up = 20 * math.log10(1 / mag)       # more gain would push the loop to instability
            else: low = 20 * math.log10(1 / mag)            # less gain would
    # the plant is unstable: at zero frequency the loop is real and negative, and the gain may fall until |L(0)| = 1
    low = 20 * math.log10(abs(L(1e-9, **kw)))
    if up is None: up = float('inf')
    peak = max(abs(1 - 2 / (1 + L(w, **kw))) for w in WS)
    return wc, pm, low, up, 2 / peak

wc, pm, low, up, dm = margins(servo=1e6)
print(f"reference loop, ideal actuator and no delay: crossover {wc:.2f} rad/s, phase margin {pm:.1f} deg, gain margin {low:.1f} dB down and {up:.1f} dB up, disk margin {dm:.2f}")
wc, pm, low, up, dm = margins()
print(f"with the 40 rad/s servo: crossover {wc:.2f} rad/s, phase margin {pm:.1f} deg, gain margin {low:.1f} dB down and {up:.1f} dB up, disk margin {dm:.2f}")
print(f"the delay margin of this loop: {math.radians(pm) / wc:.3f} s")

# the policy: requirements for the nominal and for the dispersed loop (illustrative values chosen for this course)
REQ_NOM = dict(pm=45.0, gm=6.0, dm=0.50)
REQ_DISP = dict(pm=30.0, gm=4.5, dm=0.35)
print(f"policy: nominal phase margin >= {REQ_NOM['pm']:.0f} deg, gain margin >= {REQ_NOM['gm']} dB each way, disk margin >= {REQ_NOM['dm']}; dispersed: {REQ_DISP['pm']:.0f} deg, {REQ_DISP['gm']} dB, {REQ_DISP['dm']}")

# the envelope: mass (authority b between 11.77 and 19.62), servo speed, delay, and the divergence rate
grid = []
for b in (11.77, 15.0, 19.62):
    for servo in (25.0, 40.0, 60.0):
        for delay in (0.01, 0.03, 0.06):
            for lam in (0.6, 0.7, 0.9):
                grid.append((b, servo, delay, lam))
res = []
for (b, servo, delay, lam) in grid:
    wc, pm, low, up, dm = margins(b=b, servo=servo, delay=delay, lam2=lam * lam)
    res.append((pm, low, up, dm, wc, b, servo, delay, lam))
print(f"envelope sweep: {len(grid)} combinations of authority, servo speed, delay and divergence rate")
for name, idx, f in (('phase margin', 0, min), ('lower gain margin', 1, min), ('upper gain margin', 2, min), ('disk margin', 3, min)):
    r = f(res, key=lambda t: t[idx])
    print(f"  worst {name:18s}: {r[idx]:6.2f}  at b={r[5]}, servo {r[6]:.0f} rad/s, delay {r[7] * 1000:.0f} ms, lambda {r[8]}")
bad = [r for r in res if r[0] < REQ_DISP['pm'] or r[1] < REQ_DISP['gm'] or r[2] < REQ_DISP['gm'] or r[3] < REQ_DISP['dm']]
print(f"  combinations that violate the dispersed requirement: {len(bad)} of {len(res)}")
for r in bad:
    print(f"    violation: b={r[5]}, servo {r[6]:.0f} rad/s, delay {r[7] * 1000:.0f} ms, lambda {r[8]}: phase margin {r[0]:.1f} deg, gain margin up {r[2]:.1f} dB, disk margin {r[3]:.2f}")
worst = min(res, key=lambda t: t[0])
nom = margins(b=B0, servo=40.0, delay=0.01, lam2=LAM2)
print(f"  nominal with a 10 ms delay: phase margin {nom[1]:.1f} deg, gain margin {nom[2]:.1f} dB down and {nom[3]:.1f} dB up, disk margin {nom[4]:.2f}")

# time domain: the worst-phase-margin corner, with the gimbal limits of the vehicle (rate 3 rad/s, position 0.15 rad)
def step(theta0, b, servo, delay, lam, T=8.0, dt=0.0005, rate_lim=3.0, pos_lim=0.15, gain=1.0):
    th, q, d = theta0, 0.0, 0.0
    nd = int(round(delay / dt)); buf = [0.0] * (nd + 1)
    peak = abs(theta0); tail = None
    for k in range(int(T / dt)):
        u = -gain * (KP * th + KD * q)
        buf.append(u); ucmd = buf.pop(0)
        ucmd = max(-pos_lim, min(pos_lim, ucmd))
        dd = servo * (ucmd - d)
        dd = max(-rate_lim, min(rate_lim, dd))
        d += dd * dt
        q += (lam * lam * th + b * d) * dt
        th += q * dt
        peak = max(peak, abs(th))
        if abs(th) > 20: return None, peak
    return th, peak
def settle(theta0, b, servo, delay, lam, T=8.0, dt=0.0005):
    """Time after which the tilt stays within 5% of its initial value, with the gimbal limits."""
    th, q, d = theta0, 0.0, 0.0
    nd = int(round(delay / dt)); buf = [0.0] * (nd + 1)
    last = 0.0; peak = abs(theta0)
    for k in range(int(T / dt)):
        buf.append(-(KP * th + KD * q)); ucmd = max(-0.15, min(0.15, buf.pop(0)))
        d += max(-3.0, min(3.0, servo * (ucmd - d))) * dt
        q += (lam * lam * th + b * d) * dt
        th += q * dt
        peak = max(peak, abs(th))
        if abs(th) > 0.05 * abs(theta0): last = (k + 1) * dt
    return last, peak
print(f"time domain with the gimbal limits (position 0.15 rad, rate 3 rad/s), tilt recovery from an initial tilt of 0.05 rad:")
for label, c in (('nominal (10 ms delay)', (B0, 40.0, 0.01, 0.7)), (f'worst phase margin corner (b={worst[5]}, servo {worst[6]:.0f}, delay {worst[7] * 1000:.0f} ms, lambda {worst[8]})', (worst[5], worst[6], worst[7], worst[8]))):
    t5, pk = settle(0.05, *c)
    print(f"  {label}: within 5% after {t5:.2f} s, peak tilt {pk:.3f} rad")

# does the time domain agree with the margin? raise the gain of the worst corner with no limits and watch the ring
print("gain raised at the worst-phase-margin corner (no gimbal limits), tilt after 8 s from 0.02 rad:")
up_gain = 10 ** (worst[2] / 20)
for k in (0.9 * up_gain, 1.1 * up_gain):
    th, peak = step(0.02, worst[5], worst[6], worst[7], worst[8], rate_lim=1e9, pos_lim=1e9, gain=k)
    print(f"  gain x{k:.2f} ({20 * math.log10(k):.1f} dB; the upper gain margin is {worst[2]:.1f} dB, x{up_gain:.2f}): " + ("diverges" if th is None or abs(th) > 0.02 else f"decays to {th:.5f} rad"))
print("a large initial tilt of 1.0 rad, at the gimbal limits:")
for label, c in (('nominal', (B0, 40.0, 0.01, 0.7)), ('worst', (worst[5], worst[6], worst[7], worst[8]))):
    t5, pk = settle(1.0, *c)
    print(f"  {label:7s}: within 5% after {t5:.2f} s, peak tilt {pk:.2f} rad")

reference loop, ideal actuator and no delay: crossover 6.13 rad/s, phase margin 64.3 deg, gain margin 30.5 dB down and inf dB up, disk margin 1.12
with the 40 rad/s servo: crossover 6.07 rad/s, phase margin 55.5 deg, gain margin 30.5 dB down and inf dB up, disk margin 1.00
the delay margin of this loop: 0.160 s
policy: nominal phase margin >= 45 deg, gain margin >= 6.0 dB each way, disk margin >= 0.5; dispersed: 30 deg, 4.5 dB, 0.35
envelope sweep: 81 combinations of authority, servo speed, delay and divergence rate
  worst phase margin      :  29.73  at b=19.62, servo 25 rad/s, delay 60 ms, lambda 0.6
  worst lower gain margin :  21.74  at b=11.77, servo 25 rad/s, delay 10 ms, lambda 0.9
  worst upper gain margin :   9.11  at b=19.62, servo 25 rad/s, delay 60 ms, lambda 0.6
  worst disk margin       :   0.52  at b=19.62, servo 25 rad/s, delay 60 ms, lambda 0.6
  combinations that violate the dispersed requirement: 3 of 81
    violation: b=19.62, servo 25 rad/s, delay 60 ms, lambda 0.6: phase margin 29.7 deg, gain margin up 9.1 dB, disk margin 0.52
    violation: b=19.62, servo 25 rad/s, delay 60 ms, lambda 0.7: phase margin 29.8 deg, gain margin up 9.1 dB, disk margin 0.52
    violation: b=19.62, servo 25 rad/s, delay 60 ms, lambda 0.9: phase margin 29.8 deg, gain margin up 9.1 dB, disk margin 0.52
  nominal with a 10 ms delay: phase margin 52.0 deg, gain margin 30.5 dB down and 24.8 dB up, disk margin 0.95
time domain with the gimbal limits (position 0.15 rad, rate 3 rad/s), tilt recovery from an initial tilt of 0.05 rad:
  nominal (10 ms delay): within 5% after 0.68 s, peak tilt 0.050 rad
  worst phase margin corner (b=19.62, servo 25, delay 60 ms, lambda 0.6): within 5% after 0.97 s, peak tilt 0.050 rad
gain raised at the worst-phase-margin corner (no gimbal limits), tilt after 8 s from 0.02 rad:
  gain x2.57 (8.2 dB; the upper gain margin is 9.1 dB, x2.86): decays to 0.00004 rad
  gain x3.14 (9.9 dB; the upper gain margin is 9.1 dB, x2.86): diverges
a large initial tilt of 1.0 rad, at the gimbal limits:
  nominal: within 5% after 1.06 s, peak tilt 1.00 rad
  worst  : within 5% after 1.83 s, peak tilt 1.00 rad

In robotics

A policy is what a design review and a flight-clearance board check.

• Requirements on components. The sweep turned a controller question into a hardware requirement: a servo faster than 25 rad/s or a loop delay under about 59 ms. This is the usual route from control margins to the specification of the actuator and of the flight computer.
• Gain scheduling. Fixed gains cost margin as the mass falls: the authority at burnout is 60% of that at liftoff, and the crossover falls with it. A scheduled gain keeps the margins in the band; the scheduling lesson M5-T2: Gain Scheduling designs it and this policy checks it over the envelope.
• Two breaking points. The same loop cut at the sensor gives a different disk margin, because the sensitivity at the input and the output differ for a multivariable loop. A policy requires both. For Hopper's lateral channel, with the LQR of module 4, the margins at the gimbal and at the position measurement are both stated.
• Margins and the analysis. The delay and servo ranges are inputs to the dispersions of M12-T2: Monte Carlo Dispersion Analysis: a campaign draws the servo speed and the delay from their ranges and the policy says which values must pass.
• When the policy is not enough. Margins are linear and local. The time-domain check addresses saturation, and the robustness tests of module 4 address structured uncertainty; neither replaces testing on the hardware.

The testing that follows is in UAV M8-T3: Verification and validation and UAV M8-T4: SIL, HIL and ground testing.

Common misconceptions

Key takeaways

1. A margin policy fixes the required margins at the nominal point and at every point of an envelope, and the loop-breaking points at which they are measured.
2. For the pitch loop on an unstable axis, the phase margin is the PD lead minus the servo lag minus the delay, and the lower gain margin is $20\log_{10}(bK_P/\lambda^2)$ = 30.5 dB.
3. The reference loop has 55.5° and a delay margin of 0.160 s; with a 10 ms delay, 52.0°.
4. A sweep of 81 combinations found one governing corner: 29.7° with a 25 rad/s servo and 60 ms of delay at full authority.
5. The time domain confirms the margin: stable at 8.2 dB of added gain, unstable at 9.9 dB, against 9.1 dB computed.
6. A violated margin often becomes a requirement on a component: servo speed or total delay.

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.

• K. J. Åström and R. M. Murray, Feedback Systems: An Introduction for Scientists and Engineers, 2nd ed., Princeton University Press, 2021, the chapters on the Nyquist criterion and robust performance. Free PDF from the authors' web site; the book is also sold in print. Free.
• S. Skogestad and I. Postlethwaite, Multivariable Feedback Control: Analysis and Design, 2nd ed., Wiley, 2005, chapter 2 on loop shaping and margins. Paid.
• P. Seiler, A. Packard and P. Gahinet, "An introduction to disk margins", IEEE Control Systems Magazine 40(5), 2020, on the disk margin as a specification. Paid (IEEE); a preprint is on arXiv.
• MIT 16.30 Feedback Control Systems (OpenCourseWare) , for gain and phase margins and the Nyquist criterion. 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.

• Before
  M12-T2: Monte Carlo Dispersion Analysis

  Choose what to disperse and how, compute the number of runs for a stated probability and confidence, and read the tails of the result.

• Lab
  M12-L2: Capstone: Land Hopper

  Write the navigation filter, the guidance law and the controller of Hopper, and pass an acceptance campaign of dispersed flights with wind, sensor errors, mass uncertainty and one beacon fault.

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