Stability and the Routh–Hurwitz criterion

Lesson 1 warned that feedback can make a system unstable, and lesson 4 assumed that the closed-loop poles were safely in the left half of the s-plane. This lesson makes stability precise and gives the first of three tests for it. The Routh–Hurwitz criterion decides, from the coefficients of the characteristic polynomial alone, how many closed-loop poles lie in the right half plane, without finding the roots. Its real value in design is that it works with symbols: it gives the range of a gain or a time constant for which a loop is stable, which is exactly the question an engineer asks when choosing a controller.

What stability means

A linear time-invariant system is bounded-input, bounded-output stable if every bounded input produces a bounded output. For systems described by rational transfer functions this is equivalent to a simple condition on the poles: every pole must have a strictly negative real part. A pole in the right half plane contributes a growing exponential. A simple pole on the imaginary axis contributes a mode that neither grows nor decays, such as a constant or a sustained sinusoid; such a system is called marginally stable, and it is not BIBO stable, because an input at the pole's frequency produces an output that grows without limit. Repeated poles on the imaginary axis grow like t or t².

The PT-1's open-loop plant has a pole at the origin, so it is marginally stable: a constant voltage produces a ramping angle. That is not a problem, because the loop is closed. What matters is the location of the closed-loop poles, the roots of the characteristic equation 1 + C(s)G(s) = 0, and those must all lie in the open left half plane. The question for the designer is which controllers achieve that.

A necessary condition

Write the characteristic polynomial as a_n sⁿ + a_{n−1}sⁿ⁻¹ + … + a₀. If all its roots have negative real parts, the polynomial is a product of factors (s + a) with a > 0 and (s² + bs + c) with b, c > 0, and multiplying such factors gives coefficients that are all positive. So a necessary condition for stability is that every coefficient is present and has the same sign. A missing term or a sign change means at least one root is in the right half plane or on the imaginary axis. For first- and second-order polynomials the condition is also sufficient. For third order and above it is not: s³ + s² + s + 10 has all positive coefficients but two roots in the right half plane.

The Routh array

Edward Routh in 1877 and, independently, Adolf Hurwitz in 1895 found a complete test. The Routh array is built from the coefficients. Its first row holds a_n, a_{n−2}, a_{n−4}, …; its second row holds a_{n−1}, a_{n−3}, …. Each following row is computed from the two above it: each entry is a 2 × 2 determinant of the first column and the next column of those rows, divided by the first entry of the row immediately above, with the sign arranged as below. The array has n + 1 rows, the last containing a single number.

The criterion is then a counting rule: the number of roots in the right half plane equals the number of sign changes in the first column. The system is stable if and only if every entry in the first column is positive. For a third-order polynomial a₃s³ + a₂s² + a₁s + a₀ with positive coefficients, the array reduces to a single condition that is worth remembering: a₂a₁ > a₃a₀. For s³ + s² + s + 10 that is 1 > 10, false, and the first column is 1, 1, −9, 10, with two sign changes, confirming two unstable roots.

Special cases

Two things can go wrong while building the array. If a first-column entry is zero but the rest of its row is not, the division by zero can be avoided by replacing the zero with a small positive number ε and continuing; the signs as ε → 0 give the count. A zero in the first column always indicates a root in the right half plane or on the imaginary axis. If an entire row becomes zero, the polynomial has roots placed symmetrically about the origin, such as a pair on the imaginary axis. The row above the zero row then defines an auxiliary polynomial whose roots are those symmetric roots; its derivative replaces the zero row and the array continues. In design this second case is the useful one, because a pair of imaginary roots is exactly what appears at the edge of stability, and the auxiliary polynomial gives the frequency of the oscillation.

The pan servo's gain limit

Lesson 3's proportional pan servo, analysed with the first-order motor model, has the characteristic polynomial s² + 21s + 8.33K, which is stable for every positive gain K. That cannot be right physically, and it is not: the model omitted the electrical pole. With the full plant 16 667/(s(s + 21.2)(s + 1979)), the characteristic polynomial is s³ + 2000.2s² + 41 975s + 16 667K. The third-order condition gives the gain limit, and the auxiliary polynomial 2000.2s² + 16 667K at that limit gives the frequency at which the loop oscillates.

Amara confirms the result numerically, building the Routh array for several gains and checking it against the roots found by the Durand–Kerner iteration. She also checks two other loops that come up in the design: the potentiometer sensor with a 10 ms filter that lesson 3 considered, and proportional-plus-integral control of position, which lesson 6 will want for its disturbance rejection.

import math, cmath
def routh(p, eps=1e-9):
    """Routh array for coefficients p (highest power first). Returns the first column."""
    n = len(p)
    rows = [p[0::2], p[1::2] + [0.0] * (len(p[0::2]) - len(p[1::2]))]
    for _ in range(n - 2):
        a, b = rows[-2], rows[-1]
        if abs(b[0]) < eps: b = [eps] + b[1:]         # zero in the first column
        new = [(b[0] * a[k + 1] - a[0] * b[k + 1]) / b[0] for k in range(len(a) - 1)] + [0.0]
        rows.append(new)
    return [r[0] for r in rows]
def roots(p, iters=500):
    """All roots of a polynomial by the Durand-Kerner iteration."""
    a = [c / p[0] for c in p]; n = len(a) - 1
    z = [(0.4 + 0.9j) ** k for k in range(n)]
    f = lambda x: sum(c * x ** (n - i) for i, c in enumerate(a))
    for _ in range(iters):
        z = [zi - f(zi) / math.prod(zi - zj for j, zj in enumerate(z) if j != i) for i, zi in enumerate(z)]
    return sorted(z, key=lambda r: r.real, reverse=True)
def sign_changes(col): return sum(1 for x, y in zip(col, col[1:]) if x * y < 0)
# 1. pan servo, P control, full model with the electrical pole: 16667 K / (s (s + 21.2)(s + 1979))
for K in (100.0, 2000.0, 5000.0, 6000.0):
    p = [1.0, 21.2 + 1979.0, 21.2 * 1979.0, 16667.0 * K]
    col = routh(p)
    print(f"K = {K:6.0f} V/rad: first column signs {['+' if c > 0 else '-' for c in col]}, RHP poles {sign_changes(col)}, "
          f"rightmost roots {roots(p)[0].real:.1f} ± j{abs(roots(p)[0].imag):.1f}")
Kmax = (21.2 + 1979.0) * 21.2 * 1979.0 / 16667.0
print(f"limit K < {Kmax:.0f} V/rad, oscillating at {math.sqrt(21.2*1979.0):.0f} rad/s")
# 2. potentiometer with a 10 ms filter in the feedback path
Kmax2 = 1.21 * 21.0 / (0.01 * 8.33)
print(f"with the 10 ms sensor lag: K < {Kmax2:.0f} V/rad, oscillating at {math.sqrt(21/0.01):.1f} rad/s")
# 3. PI position control: s^3 + 21 s^2 + 8.33 Kp s + 8.33 Kp / Ti
for Ti in (0.02, 0.04, 0.06, 0.2):
    col = routh([1.0, 21.0, 8.33 * 27, 8.33 * 27 / Ti])
    print(f"PI with Kp = 27, Ti = {Ti:.4f} s: RHP poles {sign_changes(col)}")
# 4. cabinet with 90 s dead time under proportional control: exact limit versus Pade approximations
K, tau, th = 0.1, 3000.0, 90.0
lo, hi = 1e-4, 1.0
for _ in range(100):                                      # phase crossover by bisection
    w = (lo + hi) / 2
    ph = -math.atan(tau * w) - th * w
    lo, hi = (w, hi) if ph > -math.pi else (lo, w)
Ku = math.sqrt(1 + (tau * w) ** 2) / K
print(f"cabinet exact: ultimate gain {Ku:.0f} W/K at period {2*math.pi/w:.0f} s")
# first-order Pade: (tau s + 1)(th/2 s + 1) + K Kc (1 - th/2 s) = 0
Kc1 = (tau + th / 2) / (K * th / 2)
print(f"first-order Pade predicts {Kc1:.0f} W/K")
# second-order Pade: e^{-th s} ~ (1 - th s/2 + th^2 s^2/12) / (1 + th s/2 + th^2 s^2/12)
def unstable(Kc):
    d = [th**2 / 12, th / 2, 1.0]; n = [th**2 / 12, -th / 2, 1.0]
    den = [tau * d[0], tau * d[1] + d[0], tau * d[2] + d[1], d[2]]
    cl = [x + K * Kc * y for x, y in zip(den, [0.0] + n)]
    return sign_changes(routh(cl)) > 0
lo, hi = 1.0, 2000.0
for _ in range(60):
    mid = (lo + hi) / 2
    lo, hi = (lo, mid) if unstable(mid) else (mid, hi)
print(f"second-order Pade predicts {lo:.0f} W/K")

K =    100 V/rad: first column signs ['+', '+', '+', '+'], RHP poles 0, rightmost roots -10.4 ± j27.1
K =   2000 V/rad: first column signs ['+', '+', '+', '+'], RHP poles 0, rightmost roots -6.3 ± j129.4
K =   5000 V/rad: first column signs ['+', '+', '+', '+'], RHP poles 0, rightmost roots -0.1 ± j204.1
K =   6000 V/rad: first column signs ['+', '+', '-', '+'], RHP poles 2, rightmost roots 2.0 ± j223.4
limit K < 5035 V/rad, oscillating at 205 rad/s
with the 10 ms sensor lag: K < 305 V/rad, oscillating at 45.8 rad/s
PI with Kp = 27, Ti = 0.0200 s: RHP poles 2
PI with Kp = 27, Ti = 0.0400 s: RHP poles 2
PI with Kp = 27, Ti = 0.0600 s: RHP poles 0
PI with Kp = 27, Ti = 0.2000 s: RHP poles 0
cabinet exact: ultimate gain 530 W/K at period 356 s
first-order Pade predicts 677 W/K
second-order Pade predicts 534 W/K

The pan servo is stable up to 5035 V/rad, and at 5000 V/rad its dominant roots sit at −0.1 ± j204, a barely damped oscillation at 205 rad/s, about 33 Hz. At 6000 V/rad the first column changes sign twice and the dominant pair has crossed into the right half plane. The limit is far above the 133 V/rad of the cascade design, so the electrical pole is not a practical constraint. The potentiometer's 10 ms filter is a different matter: it cuts the limit to 305 V/rad, with oscillation at 46 rad/s, which would have made the cascade gains marginal. That result, more than resolution, is why Amara chose the encoder.

The PI result is a useful design rule. With C(s) = K_p(1 + 1/(T_is)), the characteristic polynomial is s³ + 21s² + 8.33K_ps + 8.33K_p/T_i, and the third-order condition 21 × 8.33K_p > 8.33K_p/T_i reduces to T_i > 1/21 = 0.0476 s, independent of K_p. The integral time must be longer than the motor's mechanical time constant, or the loop is unstable at any gain. The code confirms this, with two right-half-plane roots for T_i = 0.02 s and 0.04 s and none for 0.06 s and 0.2 s. Lesson 11 shows that a good integral time is several times longer still.

Delay and the cabinet

The Routh criterion applies to polynomials, and a pure time delay, with transfer function e^{−θs}, is not a polynomial. The CC-1 cabinet, with its 90 s sensor delay, needs an approximation. The usual one is the Padé approximation, a ratio of polynomials that matches the delay's frequency response at low frequency. The first-order Padé form is (1 − θs/2)/(1 + θs/2); the second-order form adds the terms ±θs/2 + θ²s²/12. Both have the delay's exact magnitude of one at all frequencies, and their phase matches well up to a frequency of about 1/θ and 2/θ respectively.

The code compares the exact gain limit for the cabinet under proportional control with the two approximations. The exact limit, found by locating the frequency at which the loop's phase reaches −180°, a method that lesson 9 develops, is 530 W/K, with an oscillation period of 356 s. The first-order Padé approximation predicts 677 W/K, a dangerous overestimate of 28 %, because the oscillation frequency is high enough that the approximation's phase is poor there. The second-order form predicts 534 W/K, within 1 %. When delay matters, use a higher-order approximation or the frequency-domain methods of lesson 9, which handle delay exactly.

Stability in practice

A loop that is stable in theory may still be unusable. Poles just inside the left half plane give a lightly damped, ringing response, and a small change in the plant can push them across. The Routh criterion answers a yes-or-no question; it does not say how close the loop is to instability. One way to use it for margin is to shift the imaginary axis: substituting s = z − σ into the characteristic polynomial and applying the test to the result checks whether all poles lie to the left of −σ, which guarantees a minimum decay rate. For the PT-1 the settling specification of lesson 4 requires all poles to the left of −16. Gain and phase margins, introduced in lesson 9, are the more common measures, because they can be read from measured frequency responses.

Digital implementation adds delay that the continuous model does not show. The PT-1 controller samples every millisecond and applies its output through a zero-order hold, which on average delays the signal by half a sample, 0.5 ms, and the computation itself takes some fraction of a sample. At the cascade design's speeds this costs little, but a design that relied on stability margins near 200 rad/s would be in trouble. Lessons 10 and 11 account for it.