From previous lessons you know that for a continuous-time LTI closed-loop system, asymptotic stability is determined by the locations of the poles of the closed-loop transfer function \( G_{\text{cl}}(s) \). For a unity-feedback loop with plant \( G(s) = \dfrac{N(s)}{D(s)} \), we have

\[ G_{\text{cl}}(s) = \frac{G(s)}{1+G(s)} = \frac{N(s)}{D(s)+N(s)}. \]

The characteristic equation is \( D(s)+N(s)=0 \), and the corresponding characteristic polynomial \( p(s) \) can be written as

\[ p(s) = a_n s^n + a_{n-1} s^{n-1} + \dots + a_1 s + a_0, \quad a_n > 0. \]

Direct numerical computation of all roots of \( p(s) \) can be expensive and is not convenient for hand analysis or parameter design. The Routh–Hurwitz criterion offers a purely algebraic test: from the coefficients \( a_i \) alone, we can determine how many roots lie in the right half-plane (RHP) and thus whether the closed-loop system is stable.

flowchart TD
  A["Closed-loop model G_cl(s)"] --> B["Characteristic polynomial p(s)"]
  B --> C["Arrange coefficients into Routh array"]
  C --> D["Inspect signs in first column"]
  D --> E["Count sign changes"]
  E --> F["Number of RHP roots = sign changes"]
  F --> G["Conclude: stable / unstable / marginal"]
        

Later design-oriented chapters (root locus, PID tuning, loop-shaping) will frequently rely on the Routh–Hurwitz test to guarantee that proposed controller parameters preserve closed-loop stability.

Consider a real polynomial \( p(s) \) of degree \( n \):

\[ p(s) = a_n s^n + a_{n-1} s^{n-1} + \dots + a_1 s + a_0, \quad a_n > 0. \]

A necessary condition for all roots of \( p(s) \) to lie in the open left half-plane (LHP) is that all coefficients be positive: \( a_i > 0 \) for all \( i \). If any coefficient is zero or negative, \( p(s) \) must have at least one root in the closed right half-plane (RHP) or on the imaginary axis.

The Hurwitz approach introduces an associated matrix built from the coefficients. For a 4th-order polynomial \( p(s) = a_4 s^4 + a_3 s^3 + a_2 s^2 + a_1 s + a_0 \), the Hurwitz matrix is

\[ H = \begin{bmatrix} a_3 & a_1 & 0 & 0 \\ a_4 & a_2 & a_0 & 0 \\ 0 & a_3 & a_1 & 0 \\ 0 & a_4 & a_2 & a_0 \end{bmatrix}. \]

Let \( \Delta_k \) denote the leading principal minor of order \( k \) of \( H \):

\[ \Delta_k = \det H_k, \quad H_k = \text{upper-left } k \times k \text{ block of } H. \]

For the 4th-order case:

\[ \begin{aligned} \Delta_1 &= a_3, \\[4pt] \Delta_2 &= \begin{vmatrix} a_3 & a_1 \\ a_4 & a_2 \end{vmatrix} = a_3 a_2 - a_4 a_1, \\[4pt] \Delta_3 &= \det \begin{bmatrix} a_3 & a_1 & 0 \\ a_4 & a_2 & a_0 \\ 0 & a_3 & a_1 \end{bmatrix}, \quad \Delta_4 = \det H. \end{aligned} \]

The Hurwitz stability theorem (for real-coefficient polynomials with \( a_n > 0 \)) states:

\[ \text{All roots of } p(s) \text{ lie in the LHP} \iff a_i > 0 \text{ for all } i \text{ and } \Delta_k > 0 \text{ for } k = 1,\dots,n. \]

In practice, directly computing all determinants \( \Delta_k \) is algebraically heavy for higher-order systems. The Routh algorithm is an efficient way to encode these determinant conditions into a simple tabular procedure.

As an illustration, consider the 3rd-order polynomial \( p(s) = a_3 s^3 + a_2 s^2 + a_1 s + a_0 \). The Hurwitz determinants are

\[ \Delta_1 = a_2, \qquad \Delta_2 = \begin{vmatrix} a_2 & a_0 \\ a_3 & a_1 \end{vmatrix} = a_2 a_1 - a_3 a_0, \qquad \Delta_3 = a_3 \Delta_2. \]

Hence, for 3rd-order asymptotic stability we must have:

\[ a_3 > 0,\quad a_2 > 0,\quad a_1 > 0,\quad a_0 > 0,\quad a_2 a_1 - a_3 a_0 > 0. \]

We will see shortly that these same inequalities appear naturally from the Routh array.

Given the polynomial \( p(s) = a_n s^n + a_{n-1} s^{n-1} + \dots + a_0 \) with \( a_n > 0 \), the Routh array is an \( (n+1) \times m \) table, where \( m = \left\lceil \dfrac{n+1}{2} \right\rceil \). The first two rows are filled using the coefficients of \( p(s) \):

\[ \begin{array}{c|cccc} s^n & a_n & a_{n-2} & a_{n-4} & \dots \\ s^{n-1} & a_{n-1} & a_{n-3} & a_{n-5} & \dots \\ s^{n-2} & r_{3,1} & r_{3,2} & r_{3,3} & \dots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array} \]

Each new row is computed from the two rows above via a simple recurrence. Let \( r_{i,j} \) denote the element in row \( i \), column \( j \) (with \( i = 1 \) for the top row \( s^n \)). For \( i \ge 3 \) and \( j \ge 1 \):

\[ r_{i,j} = \frac{ r_{i-1,1} \, r_{i-2,j+1} - r_{i-2,1} \, r_{i-1,j+1} }{ r_{i-1,1} }. \]

If an element \( r_{i-2,j+1} \) or \( r_{i-1,j+1} \) does not exist (beyond the available columns), it is treated as zero.

For a 4th-order polynomial \( p(s) = a_4 s^4 + a_3 s^3 + a_2 s^2 + a_1 s + a_0 \), the Routh array has the form:

\[ \begin{array}{c|ccc} s^4 & a_4 & a_2 & a_0 \\ s^3 & a_3 & a_1 & 0 \\ s^2 & b_1 & b_2 & 0 \\ s^1 & c_1 & c_2 & 0 \\ s^0 & d_1 & 0 & 0 \end{array} \]

where, by the recurrence,

\[ \begin{aligned} b_1 &= \frac{a_3 a_2 - a_4 a_1}{a_3}, \quad b_2 = \frac{a_3 a_0 - a_4 \cdot 0}{a_3} = a_0, \\[4pt] c_1 &= \frac{b_1 a_1 - a_3 b_2}{b_1}, \quad c_2 = 0, \\[4pt] d_1 &= \frac{c_1 b_2 - b_1 c_2}{c_1} = b_2 = a_0. \end{aligned} \]

For a 3rd-order polynomial \( p(s) = a_3 s^3 + a_2 s^2 + a_1 s + a_0 \), the Routh array is:

\[ \begin{array}{c|cc} s^3 & a_3 & a_1 \\ s^2 & a_2 & a_0 \\ s^1 & \dfrac{a_2 a_1 - a_3 a_0}{a_2} & 0 \\ s^0 & a_0 & 0 \end{array} \]

Note the term \( \dfrac{a_2 a_1 - a_3 a_0}{a_2} \) is exactly \( \Delta_2 / \Delta_1 \) from the Hurwitz analysis above. In general, there is a deep connection between the first column of the Routh array and the Hurwitz determinants; Routh’s method is a numerically convenient encoding of these positivity conditions.

flowchart TD
  A["Polynomial coefficients {a_n,...,a_0}"] --> B["Fill first row: a_n, a_(n-2), ..."]
  B --> C["Fill second row: a_(n-1), a_(n-3), ..."]
  C --> D["Apply recurrence r_{i,j} = (r_{i-1,1} r_{i-2,j+1} - r_{i-2,1} r_{i-1,j+1}) / r_{i-1,1}"]
  D --> E["Obtain full Routh array"]
  E --> F["Examine signs of first column"]
        

In the next section we use this array to count unstable poles.

Let the Routh array of \( p(s) \) have first-column entries \( r_{1,1}, r_{2,1}, \dots, r_{n+1,1} \). Assume for now that none of these entries is zero (special cases with zeros will be discussed briefly later and in a subsequent lesson).

The central result is:

\[ \text{Number of roots of } p(s) \text{ in the RHP} = \text{number of sign changes in } \bigl( r_{1,1}, r_{2,1}, \dots, r_{n+1,1} \bigr). \]

For an asymptotically stable LTI system (all poles in the LHP, no purely imaginary poles), the number of RHP roots must be zero. Therefore, a necessary and sufficient condition for such stability is:

\[ r_{i,1} > 0 \quad \text{for all } i = 1,\dots,n+1. \]

Because the first row is \( a_n \), and the last row is \( a_0 \), the Routh test automatically includes the basic requirement that all coefficients of \( p(s) \) be positive.

A proof outline (details are beyond the scope of this introductory lesson) uses the Cauchy argument principle. Consider the rational function \( F(s) = \dfrac{p'(s)}{p(s)} \). The change in the argument of \( F(s) \) along the imaginary axis can be expressed as a function of the Routh array entries, and is proportional to the difference between the number of roots in the RHP and LHP. Routh transformed the argument principle computation into the algebraic recurrence presented above; the sign changes in the first column count how many roots lie in the RHP.

For 3rd-order systems, the Routh test recovers the Hurwitz inequalities derived earlier. From the cubic array:

\[ \begin{array}{c|cc} s^3 & a_3 & a_1 \\ s^2 & a_2 & a_0 \\ s^1 & \dfrac{a_2 a_1 - a_3 a_0}{a_2} & 0 \\ s^0 & a_0 & 0 \end{array} \]

we see that stability requires:

\[ a_3 > 0,\quad a_2 > 0,\quad \frac{a_2 a_1 - a_3 a_0}{a_2} > 0,\quad a_0 > 0, \]

which is equivalent to the set \( a_i > 0 \) and \( a_2 a_1 - a_3 a_0 > 0 \) from the Hurwitz theorem.

Consider a simple robot joint modeled (after linearization) as a DC motor driving an inertial load. A common approximation for the plant (from earlier modeling chapters) is

\[ G(s) = \frac{K}{s(s+2)(s+4)}, \]

where \( K > 0 \) is a proportional controller gain that we wish to tune. With unity feedback, the closed-loop characteristic equation is:

\[ 1 + G(s) = 0 \quad \Longleftrightarrow \quad s(s+2)(s+4) + K = 0. \]

Expanding:

\[ s(s+2)(s+4) = s(s^2 + 6s + 8) = s^3 + 6s^2 + 8s, \]

so the characteristic polynomial is

\[ p(s) = s^3 + 6 s^2 + 8 s + K. \]

We want the range of gains \( K \) for which all closed-loop poles lie in the LHP.

Construct the Routh array:

\[ \begin{array}{c|cc} s^3 & 1 & 8 \\ s^2 & 6 & K \\ s^1 & \dfrac{6 \cdot 8 - 1 \cdot K}{6} & 0 \\ s^0 & K & 0 \end{array} = \begin{array}{c|cc} s^3 & 1 & 8 \\ s^2 & 6 & K \\ s^1 & \dfrac{48 - K}{6} & 0 \\ s^0 & K & 0 \end{array}. \]

Stability requires all first-column entries to be positive:

\[ \begin{aligned} &1 > 0 &\quad \text{(always true)}, \\ &6 > 0 &\quad \text{(always true)}, \\ &\frac{48 - K}{6} > 0 \quad &\Longrightarrow\quad K < 48, \\ &K > 0. && \end{aligned} \]

Thus the necessary and sufficient condition for asymptotic closed-loop stability is

\[ 0 < K < 48. \]

In robot joint control, this inequality defines a safe range of proportional gains. Larger gains improve stiffness and tracking up to a point, but if \( K \) exceeds \( 48 \), a closed-loop pole moves into the RHP and the joint position response diverges.

Two special cases frequently appear in practice:

1. Zero in the first element of a row: If \( r_{i,1} = 0 \) while other elements in the row are nonzero, the recurrence formula leads to division by zero. A common workaround is to replace \( r_{i,1} \) by a small positive parameter \( \varepsilon \), perform the computations symbolically, and take the limit as \( \varepsilon \to 0^+ \). The sign of each first-column entry in that limit still correctly indicates the RHP root count.
2. Entire row of zeros: This indicates that \( p(s) \) has symmetrical roots with respect to the imaginary axis, often including purely imaginary roots. One then constructs an auxiliary polynomial from the row above the zero row (by interpreting its entries as even/odd coefficients), differentiates it, and replaces the zero row with the coefficients of that derivative. Detailed analysis of such marginally stable cases will be continued in Lesson 5.

These special cases are important in motion-control applications where controller parameters are pushed close to instability to obtain high performance; small modeling errors can then move poles across the imaginary axis. Understanding these numerical behaviors is therefore crucial for robust robotic control.

Python, together with numpy and control-oriented libraries like python-control and roboticstoolbox, is widely used in prototyping robotic controllers. Below is a generic Routh-array implementation operating on a list of coefficients \( [a_n, a_{n-1}, \dots, a_0] \):

import numpy as np

def routh_array(coeffs, eps=1e-9):
    """
    Compute Routh array for a polynomial with real coefficients.

    coeffs: list or 1D array [a_n, ..., a_0], a_n > 0
    returns: R (2D numpy array), n_rhp (number of RHP roots by Routh)
    """
    coeffs = np.array(coeffs, dtype=float)
    n = len(coeffs) - 1  # degree
    m = int(np.ceil((n + 1) / 2.0))
    R = np.zeros((n + 1, m), dtype=float)

    # Fill first two rows: even and odd coefficients
    R[0, :len(coeffs[0::2])] = coeffs[0::2]   # a_n, a_(n-2), ...
    R[1, :len(coeffs[1::2])] = coeffs[1::2]   # a_(n-1), a_(n-3), ...

    # Build subsequent rows
    for i in range(2, n + 1):
        if abs(R[i - 1, 0]) < eps:
            # Special case: zero leading element, use epsilon-perturbation
            R[i - 1, 0] = eps
        for j in range(0, m - 1):
            a = R[i - 1, 0]
            b = R[i - 2, 0]
            c = R[i - 2, j + 1]
            d = R[i - 1, j + 1]
            R[i, j] = (a * c - b * d) / a

    # Count sign changes in first column (ignoring tiny values)
    col = R[:, 0]
    # Replace very small values by eps for sign computation
    col[np.abs(col) < eps] = eps
    signs = np.sign(col)
    sign_changes = 0
    for i in range(len(signs) - 1):
        if signs[i] * signs[i + 1] < 0:
            sign_changes += 1

    return R, sign_changes

# Example: closed-loop for a DC motor joint
import control  # python-control library

K = 10.0
G = control.tf([K], [1, 6, 8, 0])  # G(s) = K / (s(s+2)(s+4))
G_cl = control.feedback(G, 1)      # unity feedback

# Extract characteristic polynomial coefficients
num_cl, den_cl = control.tfdata(G_cl)
den_cl = np.squeeze(den_cl)  # [a_n, ..., a_0]

R, n_rhp = routh_array(den_cl)
print("Routh array:\n", R)
print("Number of RHP roots:", n_rhp)
if n_rhp == 0:
    print("Closed-loop is asymptotically stable.")
else:
    print("Closed-loop is unstable.")
      

In a robotic application using python-control and roboticstoolbox, one typically derives the joint-space plant model, designs a feedback controller, forms the closed-loop transfer function, and then calls the above routine to ensure the chosen gains keep all poles in the LHP.

In real-time robot control (e.g., via ros_control in ROS), C++ is the dominant language. The Routh test can be implemented as a lightweight utility to verify stability of SISO loops embedded in larger robotic controllers. Below is a minimal C++ implementation (using the standard library; linear algebra libraries like Eigen are not strictly required here):

#include <vector>
#include <iostream>
#include <cmath>

struct RouthResult {
    std::vector<std::vector<double>> table;
    int n_rhp;
};

RouthResult routhArray(const std::vector<double>& coeffs, double eps = 1e-9)
{
    int n = static_cast<int>(coeffs.size()) - 1; // degree
    int m = (n + 2) / 2; // ceil((n+1)/2)
    std::vector<std::vector<double>> R(n + 1, std::vector<double>(m, 0.0));

    // Fill first two rows
    int idx = 0;
    for (int j = 0; j < m && idx < static_cast<int>(coeffs.size()); ++j, idx += 2) {
        R[0][j] = coeffs[idx]; // a_n, a_(n-2), ...
    }
    idx = 1;
    for (int j = 0; j < m && idx < static_cast<int>(coeffs.size()); ++j, idx += 2) {
        R[1][j] = coeffs[idx]; // a_(n-1), a_(n-3), ...
    }

    // Build subsequent rows
    for (int i = 2; i <= n; ++i) {
        if (std::fabs(R[i - 1][0]) < eps) {
            R[i - 1][0] = eps;
        }
        for (int j = 0; j < m - 1; ++j) {
            double a = R[i - 1][0];
            double b = R[i - 2][0];
            double c = R[i - 2][j + 1];
            double d = R[i - 1][j + 1];
            R[i][j] = (a * c - b * d) / a;
        }
    }

    // Count sign changes in first column
    int sign_changes = 0;
    auto sign = [eps](double x) {
        if (x > eps) return 1;
        if (x < -eps) return -1;
        return 1; // treat near-zero as small positive
    };
    for (int i = 0; i < n; ++i) {
        int s1 = sign(R[i][0]);
        int s2 = sign(R[i + 1][0]);
        if (s1 * s2 < 0) {
            ++sign_changes;
        }
    }

    return {R, sign_changes};
}

int main()
{
    // Example: p(s) = s^3 + 6 s^2 + 8 s + K with K = 10
    std::vector<double> coeffs = {1.0, 6.0, 8.0, 10.0};
    RouthResult res = routhArray(coeffs);

    std::cout << "Routh array:\n";
    for (const auto& row : res.table) {
        for (double val : row) {
            std::cout << val << "\t";
        }
        std::cout << "\n";
    }
    std::cout << "Number of RHP roots: " << res.n_rhp << "\n";

    if (res.n_rhp == 0) {
        std::cout << "Closed-loop is asymptotically stable.\n";
    } else {
        std::cout << "Closed-loop is unstable.\n";
    }

    return 0;
}
      

In a ROS-based robot controller, such a function can be used offline (e.g., in parameter tuning tools) or in safety monitors that verify updated gains maintain stability before being applied to the live control loop.

Java is popular in educational robotics (e.g., FRC robots via WPILib) and in some industrial mechatronic controllers. The following Java method constructs a Routh array and counts RHP roots for a polynomial with real coefficients:

public class RouthHurwitz {

    public static class Result {
        public double[][] table;
        public int nRHP;
        public Result(double[][] t, int n) {
            this.table = t;
            this.nRHP = n;
        }
    }

    public static Result routhArray(double[] coeffs, double eps) {
        int n = coeffs.length - 1;
        int m = (n + 2) / 2;
        double[][] R = new double[n + 1][m];

        // First row: a_n, a_(n-2), ...
        int idx = 0;
        for (int j = 0; j < m && idx < coeffs.length; ++j, idx += 2) {
            R[0][j] = coeffs[idx];
        }

        // Second row: a_(n-1), a_(n-3), ...
        idx = 1;
        for (int j = 0; j < m && idx < coeffs.length; ++j, idx += 2) {
            R[1][j] = coeffs[idx];
        }

        // Build remaining rows
        for (int i = 2; i <= n; ++i) {
            if (Math.abs(R[i - 1][0]) < eps) {
                R[i - 1][0] = eps;
            }
            for (int j = 0; j < m - 1; ++j) {
                double a = R[i - 1][0];
                double b = R[i - 2][0];
                double c = R[i - 2][j + 1];
                double d = R[i - 1][j + 1];
                R[i][j] = (a * c - b * d) / a;
            }
        }

        // Count sign changes in first column
        int signChanges = 0;
        for (int i = 0; i < n; ++i) {
            double x1 = R[i][0];
            double x2 = R[i + 1][0];
            int s1 = (x1 > eps) ? 1 : ((x1 < -eps) ? -1 : 1);
            int s2 = (x2 > eps) ? 1 : ((x2 < -eps) ? -1 : 1);
            if (s1 * s2 < 0) {
                signChanges++;
            }
        }

        return new Result(R, signChanges);
    }

    public static void main(String[] args) {
        // Example polynomial: s^3 + 6 s^2 + 8 s + 10
        double[] coeffs = {1.0, 6.0, 8.0, 10.0};
        Result res = routhArray(coeffs, 1e-9);

        System.out.println("Routh array:");
        for (double[] row : res.table) {
            for (double val : row) {
                System.out.print(val + "\t");
            }
            System.out.println();
        }
        System.out.println("Number of RHP roots: " + res.nRHP);
    }
}
      

In Java-based robotics stacks, low-level system identification (e.g., via WPILib) can be combined with Routh-based gain checks to ensure stability of position/velocity loops designed for motors and actuators.

MATLAB is a standard tool in control engineering, with extensive support for robotics via the Robotics System Toolbox. Below is a simple script implementing the Routh array and applying it to the robot joint example:

function [R, nRHP] = routh_array(coeffs)
% ROUTH_ARRAY Construct Routh array for polynomial with real coefficients.
%   coeffs: [a_n, ..., a_0], a_n > 0

n = numel(coeffs) - 1;
m = ceil((n + 1) / 2);
R = zeros(n + 1, m);

% First two rows
R(1,1:2:end) = coeffs(1:2:end);
R(2,1:2:end) = coeffs(2:2:end);

eps_val = 1e-9;

for i = 3:n+1
    if abs(R(i-1,1)) < eps_val
        R(i-1,1) = eps_val;
    end
    for j = 1:m-1
        a = R(i-1,1);
        b = R(i-2,1);
        c = R(i-2,j+1);
        d = R(i-1,j+1);
        R(i,j) = (a*c - b*d)/a;
    end
end

% Count sign changes
col = R(:,1);
col(abs(col) < eps_val) = eps_val;
s = sign(col);
nRHP = 0;
for i = 1:n
    if s(i)*s(i+1) < 0
        nRHP = nRHP + 1;
    end
end
end

%% Example: DC motor robot joint with proportional control
K = 10;
num = K;
den = [1 6 8 0];
G = tf(num, den);       % plant
Gcl = feedback(G, 1);   % unity feedback

[~, den_cl] = tfdata(Gcl, 'v');
[R, nRHP] = routh_array(den_cl);

disp('Routh array:');
disp(R);
fprintf('Number of RHP roots: %d\n', nRHP);
if nRHP == 0
    fprintf('Closed-loop is asymptotically stable.\n');
else
    fprintf('Closed-loop is unstable.\n');
end
      

In Simulink, the same plant and controller can be modeled with standard blocks. One can then extract the linearized closed-loop system using linmod or linearize, obtain its transfer function, and call routh_array to verify stability as controller parameters are varied in a robotic arm or mobile robot model.

Wolfram Mathematica provides symbolic computation capabilities that are useful for deriving parameter conditions for stability (e.g., ranges of gain \( K \)). A simple Routh-array constructor can be programmed as follows:

Clear[RouthArray]
RouthArray[coeffs_List, eps_:10^-9] := Module[
  {n, m, R, i, j, a, b, c, d, col, signs, signChanges},
  n = Length[coeffs] - 1;
  m = Ceiling[(n + 1)/2];
  R = ConstantArray[0.0, {n + 1, m}];

  (* First two rows *)
  R[[1, 1 ;; ;; 2]] = coeffs[[1 ;; ;; 2]];
  If[Length[coeffs] > 1,
    R[[2, 1 ;; ;; 2]] = coeffs[[2 ;; ;; 2]];
  ];

  For[i = 3, i <= n + 1, i++,
    If[Abs[R[[i - 1, 1]]] < eps, R[[i - 1, 1]] = eps];
    For[j = 1, j <= m - 1, j++,
      a = R[[i - 1, 1]];
      b = R[[i - 2, 1]];
      c = R[[i - 2, j + 1]];
      d = R[[i - 1, j + 1]];
      R[[i, j]] = (a*c - b*d)/a;
    ];
  ];

  col = R[[All, 1]];
  col = col /. x_ /; Abs[x] < eps :> eps;
  signs = Sign[col];
  signChanges = Count[Partition[signs, 2, 1], {x_, y_} /; x*y < 0];

  {R, signChanges}
]

(* Example: p(s) = s^3 + 6 s^2 + 8 s + K *)
Clear[K]
coeffs = {1, 6, 8, K};
{R, nRHP} = RouthArray[coeffs];
R // MatrixForm
nRHP
      

Because Mathematica is symbolic, one can often compute inequalities on parameters (e.g., find all \( K \) such that every entry in the first column of R is positive), recovering conditions like \( 0 < K < 48 \) symbolically rather than numerically.

Problem 1 (Stability of a Third-Order System): Determine whether the system with characteristic polynomial \( p(s) = s^3 + 4 s^2 + 5 s + 2 \) is asymptotically stable.

Solution:

Construct the Routh array:

\[ \begin{array}{c|cc} s^3 & 1 & 5 \\ s^2 & 4 & 2 \\ s^1 & \dfrac{4 \cdot 5 - 1 \cdot 2}{4} & 0 \\ & \dfrac{18}{4} & 0 \\ s^0 & 2 & 0 \end{array} \]

All first-column entries are positive: \( 1, 4, \dfrac{18}{4}, 2 > 0 \). Therefore there are zero sign changes and zero RHP roots. The system is asymptotically stable.

Problem 2 (Gain Range for Stability): A unity-feedback servo has characteristic polynomial \( p(s) = s^3 + 3 s^2 + (2+K) s + K \), where \( K \) is a proportional gain. Find the range of \( K \) for which the system is asymptotically stable.

Solution:

Routh array:

\[ \begin{array}{c|cc} s^3 & 1 & 2+K \\ s^2 & 3 & K \\ s^1 & \dfrac{3(2+K) - 1 \cdot K}{3} & 0 \\ & \dfrac{6 + 2K}{3} & 0 \\ s^0 & K & 0 \end{array} \]

Stability requires:

\[ 1 > 0,\quad 3 > 0,\quad \frac{6 + 2K}{3} > 0,\quad K > 0. \]

The inequality \( \dfrac{6 + 2K}{3} > 0 \) gives \( 6 + 2K > 0 \), i.e. \( K > -3 \). Together with \( K > 0 \), the final condition is simply:

\[ K > 0. \]

Thus any positive gain yields an asymptotically stable closed-loop system in this case.

Problem 3 (Fourth-Order Polynomial): Test the stability of the system with characteristic polynomial \( p(s) = s^4 + 2 s^3 + 3 s^2 + 4 s + 5 \).

Solution:

The Routh array is:

\[ \begin{array}{c|ccc} s^4 & 1 & 3 & 5 \\ s^3 & 2 & 4 & 0 \\ s^2 & b_1 & b_2 & 0 \\ s^1 & c_1 & c_2 & 0 \\ s^0 & d_1 & 0 & 0 \end{array} \]

Compute:

\[ \begin{aligned} b_1 &= \frac{2 \cdot 3 - 1 \cdot 4}{2} = \frac{6 - 4}{2} = 1, \\ b_2 &= \frac{2 \cdot 5 - 1 \cdot 0}{2} = 5, \\ c_1 &= \frac{1 \cdot 4 - 2 \cdot 5}{1} = -6, \\ c_2 &= 0, \\ d_1 &= \frac{(-6) \cdot 5 - 1 \cdot 0}{-6} = 5. \end{aligned} \]

So the completed Routh array is

\[ \begin{array}{c|ccc} s^4 & 1 & 3 & 5 \\ s^3 & 2 & 4 & 0 \\ s^2 & 1 & 5 & 0 \\ s^1 & -6 & 0 & 0 \\ s^0 & 5 & 0 & 0 \end{array} \]

The first column is \( [1, 2, 1, -6, 5] \), with sign sequence \( +, +, +, -, + \). There are two sign changes: \( + \to - \) and \( - \to + \). Hence there are two RHP roots and the system is unstable.

Problem 4 (Zero First-Column Entry): Consider \( p(s) = s^3 + 3 s + 2 \). Use the Routh array with the \( \varepsilon \)-perturbation to determine whether the system is stable.

Solution:

Coefficients are \( [1, 0, 3, 2] \). Routh array:

\[ \begin{array}{c|cc} s^3 & 1 & 3 \\ s^2 & 0 & 2 \\ s^1 & \cdot & 0 \\ s^0 & 2 & 0 \end{array} \]

The \( s^2 \) row has zero in the first column. Replace it by \( \varepsilon \), \( \varepsilon > 0 \):

\[ \begin{array}{c|cc} s^3 & 1 & 3 \\ s^2 & \varepsilon & 2 \\ s^1 & \dfrac{\varepsilon \cdot 3 - 1 \cdot 2}{\varepsilon} & 0 \\ & 3 - \dfrac{2}{\varepsilon} & 0 \\ s^0 & 2 & 0 \end{array} \]

For very small \( \varepsilon \), the term \( 3 - \dfrac{2}{\varepsilon} \) is negative (large negative magnitude). Thus the first column is approximately \( [1, \varepsilon, \text{large negative}, 2] \), with signs \( +, +, -, + \) and two sign changes. Therefore there are two RHP roots and the system is unstable.

Problem 5 (Linking Routh and Hurwitz for a Cubic): For \( p(s) = a_3 s^3 + a_2 s^2 + a_1 s + a_0 \) with all \( a_i > 0 \), show explicitly that the condition \( a_2 a_1 - a_3 a_0 > 0 \) arises both from the Hurwitz determinants and from the positivity of the Routh array first column.

Solution:

From the Hurwitz construction, we have \( \Delta_1 = a_2 \), \( \Delta_2 = a_2 a_1 - a_3 a_0 \). Hurwitz stability requires \( \Delta_1 > 0 \) and \( \Delta_2 > 0 \), giving \( a_2 a_1 - a_3 a_0 > 0 \).

From the Routh array,

\[ \begin{array}{c|cc} s^3 & a_3 & a_1 \\ s^2 & a_2 & a_0 \\ s^1 & \dfrac{a_2 a_1 - a_3 a_0}{a_2} & 0 \\ s^0 & a_0 & 0 \end{array}, \]

positivity of the first-column entries requires \( \dfrac{a_2 a_1 - a_3 a_0}{a_2} > 0 \), and since \( a_2 > 0 \), this is equivalent to \( a_2 a_1 - a_3 a_0 > 0 \). Thus, for cubic systems the Hurwitz determinants and the Routh first-column positivity condition encode the same inequality.

In this lesson we introduced the Routh–Hurwitz stability criterion for continuous-time LTI systems. Starting from the characteristic polynomial of the closed-loop transfer function, we defined Hurwitz determinants and then developed the Routh array as a computationally efficient tool that avoids explicit root calculation.

We established that the number of sign changes in the first column of the Routh array equals the number of right-half-plane poles, and that asymptotic stability is equivalent to positivity of all first-column entries (when there are no imaginary-axis roots). We applied the method to obtain gain ranges for stable proportional control of a robot joint, discussed numerical special cases (zeros in the Routh array), and implemented the algorithm in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, with connections to control software used in robotics. Subsequent lessons on relative stability and special cases will build on these foundations.