The Routh–Hurwitz stability criterion is a cornerstone of classical control theory, providing a direct algebraic method to determine whether a linear time-invariant system is stable without solving for its poles. For integer-order systems—where the characteristic polynomial has only integer powers of the Laplace variable s—the test is both efficient and elegant. However, a growing number of real-world applications, from viscoelastic damping to electrochemical impedance spectroscopy, are described by fractional-order differential equations. These systems introduce derivatives of non-integer order, and their characteristic equations contain fractional powers of s, rendering the traditional Routh array impossible to construct directly. Engineers and researchers must therefore adapt the underlying philosophy of the Routh–Hurwitz analysis to handle fractional-order dynamics. This article explains how to extend those ideas, presenting practical methods that preserve the simplicity of the original criterion while accommodating the richer behavior of fractional systems.

Understanding the Classical Routh–Hurwitz Criterion

Before exploring extensions, it is essential to recap the classical method. For a system with a characteristic polynomial

ansn + an−1sn−1 + … + a0 = 0

the Routh–Hurwitz criterion constructs a table (the Routh array) from the coefficients. The stability condition is that every element in the first column of the array must have the same sign (usually positive). If any sign changes occur, the number of sign changes equals the number of right-half-plane poles. The beauty of the test lies in its purely algebraic nature: no root‑finding, no complex analysis. For integer-order systems, this works perfectly because polynomials have a finite number of roots, and the left-half-plane condition is unambiguous.

The criterion also extends to systems with time delays (via the Padé approximation) and to discrete-time systems (using the bilinear transform), but all those extensions rely on the polynomial framework. Fractional-order systems break that framework, requiring a fresh perspective.

Fractional-Order Systems: A Brief Overview

A fractional-order system is one whose dynamics are described by a fractional differential equation of the form

an Dαn y(t) + … + a0 Dα0 y(t) = bm Dβm u(t) + … + b0 Dβ0 u(t)

where the differentiation orders αi and βi are not necessarily integers. The Laplace transform of such an equation (under zero initial conditions) yields a transfer function containing terms like sγ with γ fractional. For example, a simple fractional-order integrator has the transfer function 1/sα with 0 < α < 2. These systems exhibit properties not seen in integer-order systems: infinite‑dimensional state‑space representation, non‑exponential relaxation, and long‑range memory.

Fractional calculus has found applications in control of flexible structures, modeling of biological tissues, battery impedance characterization, and even financial time series. The need to check stability for such systems is acute, yet the classical Routh–Hurwitz table cannot be built because the characteristic equation is not a polynomial.

Why the Classical Routh–Hurwitz Test Fails for Fractional Systems

The obstacle is fundamental: a fractional-order characteristic equation, such as

s2.5 + 3 s1.2 + 2 = 0

is not a polynomial. The powers of s are not integer multiples of each other. The Routh array requires integer powers so that a systematic elimination can be performed. Moreover, fractional powers lead to branch cuts in the complex plane; the concept of “right-half-plane pole” becomes more subtle because a fractional power can have multiple Riemann sheets. A fractional system may have an infinite number of poles (if the order is irrational) or a finite number if the orders are commensurate.

Because the Routh–Hurwitz criterion is essentially an algebraic test for the Hurwitz property of a polynomial, it cannot be applied directly to a non‑polynomial expression. Instead, modifications are needed that preserve the spirit of the test—determining stability by examining some transformed algebraic condition—while accounting for the fractional nature.

The Matignon Stability Criterion: A Direct Algebraic Replacement

The most well‑known extension is the Matignon stability criterion (also called the Matignon theorem). It provides a necessary and sufficient condition for a large class of fractional‑order systems with commensurate orders. A fractional system is said to have commensurate orders if all fractional powers in the characteristic equation are integer multiples of a base order α (e.g., α, 2α, 3α, …). Many practical systems, especially those derived from fractional PID controllers or from fractional viscoelastic models, are commensurate.

Matignon criterion: For a fractional‑order system with characteristic equation

P(s) = an snα + an−1 s(n−1)α + … + a0 = 0

where α is a positive real number, the system is stable if and only if

|arg(si)| > απ/2 for every root si of P.

In other words, the roots must lie in a region of the complex plane outside a sector of angle απ/2 centered on the negative real axis. This condition reduces to the classical left‑half‑plane condition when α = 1 (since then απ/2 = π/2), but for α < 1 the stability region is a wedge.

The Matignon criterion is powerful because it transforms the stability problem back into a root‑location problem for the transformed polynomial P in the variable w = sα. One can apply the Routh–Hurwitz test to the polynomial P(w) after substituting w = sα, but the condition on the argument of s is not captured by the classical Routh array alone. However, for commensurate systems, the stability condition can be checked by first constructing the Routh array for the polynomial P(w) and then verifying the argument condition.

This leads to a practical procedure:

  1. Identify the base fractional order α so that the characteristic equation becomes a polynomial in sα.
  2. Define w = sα and rewrite the characteristic equation as P(w) = 0, a conventional polynomial.
  3. Apply the classical Routh–Hurwitz test to P(w). If any sign changes occur in the first column, the system is unstable under the classical polynomial definition in w.
  4. Check the argument condition: For each root wi, compute the corresponding si = wi1/α and verify that |arg(si)| > απ/2. This step may require numerical computation because the roots wi are generally complex.

The Matignon criterion is often used as the primary tool for stability analysis of commensurate fractional‑order systems. For non‑commensurate systems (where the fractional orders are not multiples of a common base), alternative methods such as the frequency‑domain approach or numerical inversion must be used.

Alternative Analytical Approaches

The Mittag–Leffler Function and State‑Space Representation

Fractional‑order systems can be represented in a pseudo‑state‑space form using the fractional derivative of order α. The stability condition becomes a constraint on the eigenvalues of the system matrix. Specifically, for a linear fractional‑order system described by

CDα x(t) = A x(t)

the system is stable if and only if

|arg(λ(A))| > απ/2

where λ(A) are the eigenvalues of the matrix A. This is the direct multi‑variable analogue of the Matignon criterion. For commensurate systems, this condition reduces to the same wedge‑shaped stability region. Engineers can compute the eigenvalues and check their arguments using standard numerical linear algebra.

Transformation to an Integer‑Order System via Variable Substitution

Another analytical path is to introduce a new variable z = sα and then apply the classical Routh–Hurwitz test to the transformed polynomial. However, this approach only checks the “polynomial part” of the stability. A necessary condition is that the polynomial in z has all roots in the left half‑plane, but the sufficient condition requires the additional argument test. In practice, engineers often perform both: first run the Routh array on P(z) to detect obvious polynomial instability; then, for the remaining stable candidates, apply the argument condition.

Frequency‑Domain Techniques

The stability can also be assessed by the Nyquist or Bode criteria adapted to fractional‑order transfer functions. The open‑loop transfer function G(s) of a fractional‑order system is a non‑rational function, but its frequency response can be evaluated numerically. The Nyquist stability criterion still holds: the number of encirclements of the −1 point determines closed‑loop stability. However, the direct application of the Routh–Hurwitz method is often preferred because it avoids plotting and provides an algebraic condition.

Numerical and Computational Tools

Given the complexity of fractional‑order systems, engineering practice relies heavily on software tools. The FOMCON toolbox (Fractional‑Order Modeling and Control) for MATLAB is a popular choice. It provides functions to compute stability checks using the Matignon criterion, to construct approximations (e.g., Oustaloup recursive approximation) of fractional operators, and to perform time‑domain simulations.

Other tools include:

  • CRONE toolbox (Commande Robuste d’Ordre Non Entier): developed at the University of Bordeaux, it supports fractional‑order controller design and stability analysis.
  • Python libraries: fractional (for fractional calculus) and control (for general control system analysis) can be used together to implement the Matignon criterion.
  • Symbolic computation in Mathematica or Maple can handle algebraic manipulations of fractional powers and root‑finding in the complex plane.

For a quick check, many engineers use the following numerical algorithm:

  1. Define the characteristic equation as a symbolic expression in s with fractional powers.
  2. Convert to a polynomial in w = sα if the system is commensurate.
  3. Compute the roots of the polynomial in w using a numerical solver.
  4. For each root, compute the corresponding s and its argument.
  5. Check if |arg(s)| > απ/2 for all roots.

This algorithm can be implemented in a few lines of code and is robust for most commensurate systems. For non‑commensurate systems, one must rely on frequency‑domain methods or on grid‑based numerical evaluation of the characteristic equation in the complex plane.

Practical Step‑by‑Step Procedure for Engineers

Below is a pragmatic workflow that combines the classical Routh–Hurwitz spirit with the needed fractional modifications.

  1. Obtain the characteristic equation from the system’s transfer function. For a fractional‑order system, this will involve terms like sγ with γ non‑integer.
  2. Check for commensurability. Determine if all fractional orders are integer multiples of a common base α. For example, orders 0.5, 1.5, 2.0 are multiples of 0.5; orders 0.7 and 1.2 are not multiples of a single base (unless the base is 0.1, but that is rarely rational). If the system is commensurate, proceed to the polynomial transformation. If not, skip to step 6 (frequency‑domain method).
  3. Define the polynomial in w = sα. Write the characteristic equation as P(w) = 0. For instance, s2.5 + 3 s1.2 + 2 is not commensurate (2.5 and 1.2 are not multiples of a common base), so this step would fail. A better example: s1.5 + 3 s0.5 + 2 = 0 with α=0.5 gives w3 + 3w + 2 = 0.
  4. Apply the classical Routh–Hurwitz test to P(w). Build the Routh array. If any sign change occurs in the first column, the system is unstable (under the polynomial interpretation). If all first‑column entries have the same sign, proceed to the next step.
  5. Verify the argument condition. Compute the roots of P(w). For each root wi, compute si = wi1/α. The complex power is multivalued; use the principal value. Then check |arg(si)| > απ/2. If this holds for all roots, the fractional system is stable. If any root violates the condition, the system is unstable.
  6. For non‑commensurate systems, use a frequency‑domain approach. Evaluate the characteristic equation along the boundary of the stability region, which for fractional systems is a sector. The Nyquist plot of the open‑loop transfer function can indicate closed‑loop stability. Alternatively, use numerical root‑finding in the s‑plane by discretizing the fractional operator, for example using the Oustaloup approximation. Then apply the classical Routh–Hurwitz test to the approximated integer‑order system. This method gives an approximate result and should be validated with simulations.

Many engineers find that the commensurate case covers a large portion of practical fractional‑control designs, especially fractional PID controllers (where the derivative and integral orders are multiples of a common α, often 0.5 or 0.1). The procedure above is therefore widely applicable.

Illustrative Case Study: Fractional‑Order PID Controller

Consider a fractional‑order PID controller (PIαDβ) controlling a first‑order plant. The overall closed‑loop characteristic equation may be

s2.5 + 3 s1.2 + 0.5 s0.1 + 1 = 0

This system is non‑commensurate because the fractional orders 2.5, 1.2, 0.1 are not integer multiples of a common base. We cannot directly apply the Matignon criterion. Instead, we approximate the fractional operators using an Oustaloup filter of order, say, 5 over a frequency range [0.01, 100] rad/s. This yields a high‑order integer‑order transfer function. The approximated characteristic polynomial will be of order 5×5 = 25 for each fractional term, but we can reduce the order using balanced truncation. Then the classical Routh–Hurwitz test can be applied to the reduced polynomial. A stability check via simulation confirms the result. This approximate method is the de facto standard in industry when analytic tools are unavailable.

In contrast, if the controller is designed with commensurate orders (e.g., α = 0.5, so orders are 0.5, 1.0, 1.5, 2.0), the system becomes commensurate. Then we can use the Matignon criterion directly, obtaining an exact necessary and sufficient condition (modulo the branch‑cut interpretation). This illustrates the advantage of designing fractional controllers with commensurate orders when possible.

Limitations and Ongoing Research

The extension of Routh–Hurwitz analysis to fractional‑order systems is not a single formula but a collection of techniques, each with its own domain of applicability. The Matignon criterion is exact only for commensurate systems. For non‑commensurate systems, the approximated integer‑order approach may introduce errors, especially if the system has lightly damped high‑frequency modes. Researchers are exploring direct algebraic criteria for non‑commensurate systems using the concept of fractional‑order Routh arrays (e.g., using the fractional‑order variable transformation s = (z)1/α and applying a modified array). Another promising direction is the use of linear matrix inequalities (LMIs) adapted to fractional‑order systems.

Furthermore, the stability region for fractional systems with irrational orders is not a simple wedge; it can become a more complicated region in the Riemann surface. Practical engineers usually rely on numerical evaluation of the characteristic equation’s magnitude and phase along the stability boundary, which is computationally intensive but yields reliable results for specific parameter values.

Conclusion

Extending the Routh–Hurwitz analysis to fractional‑order systems requires moving beyond the classical polynomial framework. The Matignon stability criterion for commensurate systems re‑casts the problem into a root‑argument condition that can be checked with the help of the classical Routh array and numerical root‑finding. For non‑commensurate systems, frequency‑domain methods and numerical approximations remain the primary tools. By understanding both the theoretical basis and the practical numerical workflow, engineers can confidently perform stability analysis for a wide range of fractional‑order systems—bridging the gap between classical control and emerging fractional‑order applications.

For further reading, readers may consult the Wikipedia entry on the Routh–Hurwitz criterion and the fractional calculus overview. The original work by Matignon is detailed in this paper (DOI link). Additionally, the FOMCON toolbox for MATLAB provides practical implementation of the methods described here.