Introduction to the Routh-Hurwitz Criterion

The Routh-Hurwitz criterion is a classical algebraic method used in control engineering to determine the absolute stability of a linear time-invariant (LTI) system. Given the characteristic polynomial of a system's closed-loop transfer function, this criterion provides a necessary and sufficient condition for stability without requiring explicit calculation of the poles. The core idea is simple: if all roots of the characteristic equation lie in the left half of the complex plane, the system is stable. The Routh-Hurwitz test checks this by constructing a tabular array (the Routh array) from the polynomial coefficients and analyzing sign changes in its first column. The number of sign changes equals the number of poles with positive real parts, directly indicating instability.

For engineering students, mastering this technique is essential for designing controllers and evaluating system robustness. However, manual construction of the Routh array becomes tedious and error-prone for high-order polynomials. MATLAB offers a powerful environment to automate the process, enabling students to focus on interpretation rather than arithmetic. This article provides a comprehensive guide to implementing the Routh-Hurwitz criterion in MATLAB, complete with code examples, special-case handling, and practical insights. MATLAB's built-in stability functions can also be used for comparison, but writing your own implementation deepens understanding.

Why Implement in MATLAB?

Manual Routh-Hurwitz analysis is feasible for second- and third-order systems, but as the order increases, the number of entries in the array grows quadratically, and arithmetic mistakes become common. MATLAB eliminates such human error and allows rapid iteration when analyzing multiple design candidates. Moreover, MATLAB can handle floating-point arithmetic precisely, detect near-zero entries that cause numerical issues, and integrate with other analysis tools like root locus or Bode plots. By implementing the algorithm yourself, you gain a flexible tool that can be extended to symbolic calculations, parameterized studies, or educational demonstrations. Symbolic Math Toolbox can also be employed to keep parameters symbolic during stability analysis.

Step-by-Step Implementation

Defining the Characteristic Polynomial

The first step is to define the characteristic polynomial of the system. In control engineering, the characteristic equation is derived from the denominator of the closed-loop transfer function, typically expressed as aₙ sⁿ + aₙ₋₁ sⁿ⁻¹ + ... + a₀ = 0. The coefficients are stored in a vector coeffs in descending order of powers. For example, the polynomial s³ + 3s² + 10s + 20 is represented as [1 3 10 20]. Ensure that the polynomial is complete (no missing powers) by inserting zeros for missing terms.

Constructing the Routh Array Algorithm

The Routh array is a matrix whose rows correspond to decreasing powers of s (starting from the highest power) and columns correspond to alternating coefficients. The first two rows are filled directly from the polynomial coefficients: the first row contains coefficients of even powers (starting with the highest power), and the second row contains coefficients of odd powers. Subsequent rows are computed using the formula:

routh(i,j) = - (routh(i-2,1)*routh(i-1,j+1) - routh(i-1,1)*routh(i-2,j+1)) / routh(i-1,1)

This formula is applied for each column index j until the row entries become zero or the row length is exhausted. In MATLAB, we initialize the array with zeros and fill iteratively. Special care must be taken when the first element of a row is zero, which would cause division by zero. A common workaround is to replace the zero with a small epsilon (e.g., eps) and continue, though this may introduce minor numerical inaccuracies. More robust approaches include using the auxiliary polynomial method, which we discuss later.

Analyzing the First Column for Sign Changes

Once the array is complete, extract the first column (all rows). Count how many times the sign of the entries changes when moving down the column. A sign change occurs when two consecutive elements have opposite signs (positive to negative or negative to positive). If there are no sign changes, the system is stable; otherwise, the number of sign changes equals the number of roots in the right half-plane, indicating instability.

% Example: Stable system (polynomial s^3 + 3s^2 + 10s + 20)
coeffs = [1 3 10 20]; 
n = length(coeffs);
routh = zeros(n, ceil(n/2));

% First row (even powers)
routh(1, 1:2:end) = coeffs(1:2:end);
% Second row (odd powers)
if length(coeffs) >= 2
    routh(2, 1:2:end) = coeffs(2:2:end);
end

% Fill remaining rows
for i = 3:n
    for j = 1:size(routh,2)-1
        a = routh(i-2, 1);
        b = routh(i-2, j+1);
        c = routh(i-1, 1);
        d = routh(i-1, j+1);
        if c == 0
            c = eps; % replace zero with small value to avoid division by zero
        end
        routh(i, j) = ((c * b) - (a * d)) / c;
    end
end

disp('Routh Array:')
disp(routh)
firstCol = routh(:,1);
signChanges = sum(diff(sign(firstCol)) ~= 0);
if signChanges == 0
    disp('System is stable.');
else
    disp(['System is unstable with ', num2str(signChanges), ' sign change(s).']);
end

Example: Stable System

Consider the polynomial [1 3 10 20] used in the code above. Running the script yields a Routh array with the following first column values: 1, 3, 3.333, 20. All are positive, so zero sign changes. The system is stable. You can verify using MATLAB's isstable command on a transfer function with that denominator (assuming a proper numerator).

Example: Unstable System

Now test an unstable polynomial: [1 2 3 4 5] (fourth-order). The first column might show sign changes. For instance, if the first column becomes 1, 2, -3, 5, 2, there are two sign changes (from 2 to -3, then from -3 to 5), indicating two right-half-plane poles. Such a system would have increasing oscillations or divergence. Interpreting these results helps students connect algebraic conditions with dynamic behavior.

Handling Special Cases

Zero in the First Column

When the first element of a row is zero but the row is not entirely zero, the standard algorithm breaks down. The simplest fix is to replace that zero with a very small positive number (like eps). However, this can sometimes produce incorrect results if the zero arises from a symmetrical root placement. A more rigorous method is to use the epsilon method with symbolic small parameter or to form an auxiliary polynomial from the previous row. For robust code, you can check if c == 0 and then replace it with eps. This works for most academic examples but may fail in pathological cases. Advanced implementations can detect this situation and apply the auxiliary polynomial method (see this control theory resource).

Row of Zeros

A row entirely composed of zeros indicates that the polynomial has roots that are symmetric about the origin (e.g., imaginary axis poles). In such cases, the system is marginally stable or unstable depending on the sign of the auxiliary polynomial's derivative. The standard approach is to form an auxiliary polynomial from the row above the zero row, differentiate it, and replace the zero row with the coefficients of the derivative. This procedure identifies roots on the imaginary axis or symmetrically placed real roots. Implementing this in MATLAB requires dynamic row insertion and polynomial differentiation. For brevity, many textbooks treat this as an special exercise. A practical student implementation might print a warning and stop, asking the user to apply the auxiliary polynomial manually.

Advanced Implementation with Symbolic Variables

Using the Symbolic Math Toolbox, you can create a Routh-Hurwitz function that accepts symbolic coefficients. This is useful when a parameter (e.g., a controller gain K) appears in the characteristic equation. The symbolic array will contain expressions rather than numbers, and sign analysis becomes a problem of determining for what ranges of the parameter the first-column entries are all positive. You can then solve inequalities to find the stability margin. A simple symbolic version might replace numeric entries with symbolic variables and then evaluate sign(expr) after substituting candidate values. MATLAB's symbolic documentation explains how to create and manipulate symbolic expressions. However, note that symbolic reasoning about sign changes is more complex; often you'll substitute numeric ranges.

Comparison with Other Stability Analysis Methods

The Routh-Hurwitz criterion is algebraic and does not require plotting or frequency response, making it ideal for pencil-and-paper analysis. However, MATLAB also provides other stability tools:

  • Pole-zero map (pzmap): Directly shows pole locations; stable if all poles are in the left half-plane.
  • Root locus (rlocus): Shows how poles move with a parameter, giving insight into relative stability.
  • Nyquist plot (nyquist): Uses frequency response to determine stability via the Nyquist criterion, advantageous for systems with delays.
  • Bode plot (bode): Gain and phase margins indicate stability margins.

While Routh-Hurwitz is straightforward, it cannot reveal damping ratios or resonance peaks. For comprehensive analysis, combine the Routh criterion with other MATLAB tools. For example, after verifying stability with Routh, use step to visualize transient response.

Practical Tips for Engineering Students

  • Always start by checking if all coefficients of the characteristic polynomial are positive; if not, the system is certainly unstable (necessary condition).
  • Use routh = routh(1:end, 1:end); to trim trailing zeros after construction for cleaner display.
  • For high-degree polynomials, consider using poly and roots to verify the Routh result.
  • Incorporate the code into a reusable function routhHurwitz(coeffs) that returns stability status and the array for debugging.
  • Understand the limitations: the algorithm assumes real coefficients and continuous-time systems. For discrete-time, use the bilinear transform or direct Jury stability test.

Conclusion

Implementing the Routh-Hurwitz criterion in MATLAB transforms a manual algebraic chore into an automated, reliable analysis tool. Engineering students gain deepened understanding of stability conditions and learn to handle special cases like zero rows and division errors. The code presented here serves as a solid foundation that can be extended to symbolic analysis and integrated into larger control design workflows. By combining this implementation with MATLAB's extensive visualization and analysis functions, students can rapidly evaluate system stability and focus on higher-level design decisions. As you progress, explore the MATLAB File Exchange contributions for more polished versions and consider contributing your own refinements.