Stability analysis of control loops

A controller that makes the robot faster but unstable is worse than useless; it is dangerous. Stability analysis tells us, before we touch the hardware, whether a loop will settle, how much gain we can add before it oscillates, and how much unmodeled delay or lag it can tolerate. ROB-230 introduced poles, Bode plots and gain and phase margins. This lesson adds the classical tools engineers use to analyze closed loops as a gain changes: the Routh–Hurwitz criterion, which tests stability from polynomial coefficients without finding roots; the root locus, which shows how closed-loop poles move as gain rises; and the Nyquist criterion, which connects open-loop frequency response to closed-loop stability and justifies the margins. Each is applied to the joint loops of lesson 4.

Stability of a closed loop

A linear time-invariant closed loop is stable when all roots of its characteristic polynomial, the closed-loop poles, have negative real parts. For a unity-feedback loop with controller C(s) and plant P(s), the characteristic equation is

where L(s) is the loop transfer function. Writing \( L(s) = KN(s)/D(s) \), the closed-loop poles are the roots of \( D(s) + KN(s) = 0 \). For a first or second-order polynomial with all coefficients positive, stability is automatic. From third order up it is not, which is why PID loops on joints with lags, and any loop with delay, need analysis.

The Routh–Hurwitz criterion

Routh's 1877 method decides how many roots of a polynomial lie in the right half-plane using only arithmetic on the coefficients. For \( a_ns^n + a_{n-1}s^{n-1} + \cdots + a_0 \), a necessary condition is that all coefficients be present and have the same sign. To test sufficiency, build the Routh array: the first two rows hold alternate coefficients, and each later entry is computed from the two rows above by

The number of sign changes in the first column equals the number of roots in the right half-plane. For a cubic \( a_3s^3 + a_2s^2 + a_1s + a_0 \) with positive coefficients, the array reduces to a single condition:

which is the PID stability condition \( (b + K_d)K_p > JK_i \) quoted in lesson 4. The real strength of the method is that it handles symbolic gains, so it gives the range of a gain for which a loop is stable.

The root locus

Routh tells us whether the loop is stable; the root locus shows how the closed-loop poles move as a gain K rises from zero to infinity, so we can choose K for the response we want. Evans' rules, published in 1948, let us sketch the locus by hand for \( L(s) = KN(s)/D(s) \) with n poles and m zeros:

• The locus starts (K = 0) at the open-loop poles and ends (K → ∞) at the open-loop zeros, or at infinity along asymptotes for the n − m poles that have no zero to reach.
• Points on the real axis lie on the locus if an odd number of real poles and zeros lie to their right.
• The asymptotes leave from the centroid \( \sigma_a = (\sum p_i - \sum z_i)/(n - m) \) at angles \( (2k + 1)180°/(n - m) \).
• Branches leave the real axis at breakaway points, where \( dK/ds = 0 \).
• Imaginary-axis crossings, where the loop becomes unstable, are found with Routh or by substituting \( s = j\omega \).

The root locus also shows how compensation helps. Adding a zero, as derivative action does, pulls the locus to the left and allows more gain with good damping. Adding a pole at the origin, as integral action does, pulls it to the right, which is why PID designers place the integral zero well below the crossover frequency.

The Nyquist criterion

Root locus and Routh need a polynomial model. Nyquist's criterion, from 1932, works directly with the open-loop frequency response \( L(j\omega) \), which can be measured on the robot with a sine sweep even when the model is uncertain and even when the loop contains a pure delay \( e^{-sT} \), which has no polynomial form. The Nyquist plot traces \( L(j\omega) \) in the complex plane as ω runs from −∞ to +∞. The criterion states:

where P is the number of open-loop poles in the right half-plane, N is the number of clockwise encirclements of the critical point −1 by the Nyquist plot, and Z is the number of closed-loop poles in the right half-plane. For the usual case of a stable open loop (P = 0, treating poles at the origin with a small detour), the closed loop is stable if and only if the plot does not encircle −1. The margins of ROB-230 measure how far the plot stays from −1: the gain margin is how much the gain can be multiplied before the plot reaches −1 along the negative real axis, and the phase margin is how much extra phase lag, such as delay, can be added before the point at unit magnitude rotates onto −1.

For robot control, the practical outcome is a set of design targets that keep the loop robust to modeling error: a phase margin of at least 45° and a gain margin of at least 6 dB (a factor of 2) are typical, with more margin where the plant varies a lot with payload or configuration.

Delay and the delay margin

Delays are everywhere in digital robot control: sampling and computation, communication buses, sensor processing and, for teleoperation, the network. A pure delay T leaves magnitude unchanged and adds phase lag \( \omega T \) radians, growing without limit with frequency. If a loop has phase margin \( \phi_m \) (in radians) at crossover frequency \( \omega_c \), the largest extra delay it can tolerate before instability is the delay margin

This gives a vivid design rule: faster loops, with higher crossover, tolerate proportionally less delay. A joint position loop with 20 rad/s crossover and 45° margin tolerates about 39 ms of delay; a current loop crossing over at 6,000 rad/s with the same margin tolerates only 0.13 ms, which is why current loops run in dedicated hardware next to the motor.

Analyzing loops numerically

In practice engineers combine these tools with numerical computation. The code below finds closed-loop poles across a range of gains for the example loop, locates the gain at which a pole crosses the imaginary axis, and computes the gain and phase margins from the frequency response for K = 20, confirming the hand results.

import numpy as np
den = np.array([1, 7, 10, 0])             # s(s+2)(s+5)
def poles(K):
    return np.roots(den + np.array([0, 0, 0, K]))
for K in (5, 20, 60, 70, 80):
    print(f"K={K:3d}: max real part {max(poles(K).real):+.3f}")
w = np.logspace(-2, 2, 20000)
K = 20
L = K/(1j*w*(1j*w + 2)*(1j*w + 5))
mag, ph = np.abs(L), np.unwrap(np.angle(L))*180/np.pi
ic = np.argmin(abs(mag - 1)); ip = np.argmin(abs(ph + 180))
print(f"crossover {w[ic]:.2f} rad/s, phase margin {180 + ph[ic]:.1f} deg")
print(f"phase crossover {w[ip]:.2f} rad/s, gain margin {1/mag[ip]:.2f} ({-20*np.log10(mag[ip]):.1f} dB)")
print(f"delay margin {np.radians(180 + ph[ic])/w[ic]*1000:.0f} ms")

K=  5: max real part -0.856
K= 20: max real part -0.561
K= 60: max real part -0.088
K= 70: max real part +0.000
K= 80: max real part +0.082
crossover 1.52 rad/s, phase margin 35.8 deg
phase crossover 3.16 rad/s, gain margin 3.50 (10.9 dB)
delay margin 410 ms

The largest real part crosses zero between K = 60 and 80, at 70, and the frequency-response margins agree with the hand calculation. The delay margin of about 0.4 s is large because this example loop is slow; real joint loops cross over ten to a hundred times faster.