Routh-Hurwitz Stability Criterion

Required imports

from sympy import *
from sympy.abc import s, epsilon, k

from mathprint import *

Motivation

Given the following transfer function:

\[ G(S) = \frac{ N(s) }{D(s)}\]

where \(N(s)\) is the numerator and \(D(s)\) is the denominator. Both are polynomials. The roots of the numerator defines the system stabilty.

Routh-Hurwitz method allows us to conclude the system’s stability without actually calculating the polynomial roots of the denominator.

The Routh table

Given the following system:

\(\boxed{G(s) = \frac{ \text{num}(s)}{a_ns^n + a_{n-1}s^{n-1} + \dots + a_1s + a_0}}\)

As an example, let us take \(n=4\). The RH-table is as follows:

\(s^4\)

\(a_4\)

\(a_2\)

\(a_0\)

\(s^3\)

\(a_3\)

\(a_1\)

0

\(s^2\)

\(\frac{-\left|\begin{array}{ll}a_4 & a_2 \\ a_3 & a_1\end{array}\right|}{a_3}=b_1\)

\(\frac{-\left|\begin{array}{cc}a_4 & a_0 \\ a_3 & 0\end{array}\right|}{a_3}=b_2\)

\(\frac{-\left|\begin{array}{ll}a_4 & 0 \\ a_3 & 0\end{array}\right|}{a_3}=0\)

\(s^1\)

\(\frac{-\left|\begin{array}{ll}a_3 & a_1 \\b_1 & b_2\end{array}\right|}{b_1}=c_1\)

\(\frac{-\left|\begin{array}{ll}a_3 & 0 \\ b_1 & 0\end{array}\right|}{b_1}=0\)

\(\frac{-\left|\begin{array}{ll}a_3 & 0 \\ b_1 & 0\end{array}\right|}{b_1}=0\)

\(s^0\)

\(\frac{-\left|\begin{array}{cc}b_1 & b_2 \\c_1 & 0\end{array}\right|}{c_1}=d_1\)

\(\frac{-\left|\begin{array}{ll}b_1 & 0 \\ c_1 & 0\end{array}\right|}{c_1}=0\)

\(\frac{-\left|\begin{array}{ll}b_1 & 0 \\ c_1 & 0\end{array}\right|}{c_1}=0\)
\[ \boxed{ \text{ The denominator polynomial is stable iff all elements of the first column of the Routh table are nonzero and have the same sign. } } \]

As we can see, the procedures to build a Routh table is repetitive. Thus, we will use Python SymPy to automate the the procedure.

Implementations

The implementation involves four functions:

  • simplify_line

  • make_routh_table

  • evaluate_table

  • print_table

'''
This function is to simplify one line/row by dividing all its elements 
with the greatest common divisor. However, before that, we must exclude
any element that is zero.
'''
def simplify_line(line):
    l = []
    for k in range(len(line)):
        if line[k] != 0.0:
            l.append(line[k])

    return gcd(l)

'''
This is the core function that generates the RH-table.
'''
def make_routh_table(sys):
    table = []

    # convert into fraction, 
    # take the coefficients of the denominator
    f = fraction(sys)
    den = Poly(f[1], s)
    a = den.all_coeffs()
    
    norder = len(a) - 1 # system order number
    nrow = norder + 1
    ncol = round(norder / 2)
    if ncol % 2 == 0:
        ncol = ncol + 1

    table = zeros(nrow, ncol)

    # Fill the first two rows
    table[0, 0:len(a[0::2])] = [a[0::2]]
    table[1, 0:len(a[1::2])] = [a[1::2]]
    
    for j in range(norder-1):
        r = simplify_line(table[j,:])
        table[j,:] = table[j,:] / r

        if table[j+1,:].is_zero_matrix == False:  # NOT a row of all zeros
            if (table[j+1,0]== 0):                # case 1: first column is zero, avoid division by zero
                table[j+1,0] = epsilon
        elif table[j+1,:].is_zero_matrix == True: # case2: a row of all zeros  
            expr = 0.0
            for k in range(0,ncol):
                expr =  expr + table[j,k]*s**((norder-2)-2*k)

            a = Poly(diff(expr,s), s).all_coeffs()[0::2]
            table[j+1,0:len(a)] = [a]

        for k in range(ncol-1):
            A = Matrix([[table[j,0],    table[j,k+1]],
                       [table[j+1,0],  table[j+1,k+1]]])
            
            table[j+2,k] = simplify(-A.det()/A[1,0])
        
    return table

'''
This function does not conclude the stability based on the generated Routh table. 
It simply generates and print the first-column signs which we can use to draw the conclusion.
'''
def evaluate_table(table):
    signs = []
    for k in range(len(table[:,0])):
        table[k,0] = table[k,0].subs(epsilon, 0.001) 
        if table[k,0].evalf() < 0:
            signs.append('-')
        elif table[k,0].evalf() > 0:
            signs.append('+')

    sp = ''
    cnt = 0
    for k, s in enumerate(signs):
        if (k > 0) and (s != sp):
            cnt = cnt + 1
        sp = s

    print(cnt, 'sign changes.')
    
    return signs

'''
This prints the Routh table
'''
def print_table(table):
    mprint(latex(Matrix(table)))

Examples

Example 1

\[ \frac{1}{a_0s^3 + a_1s^2 + a_2s + a_3}\]
a0, a1, a2, a3 = symbols("a0 a1 a2 a3")
table = make_routh_table(1 / (a0*s**3 + a1*s**2 + a2*s + a3))
print_table(table)
\[\begin{split}\displaystyle \left[\begin{matrix}a_{0} & a_{2} & 0\\a_{1} & a_{3} & 0\\- \frac{a_{0} a_{3}}{a_{1}} + a_{2} & 0 & 0\\a_{3} & 0 & 0\end{matrix}\right]\end{split}\]

Therefore, to guarantee stabilty:

table[2,0]>0
\[\displaystyle - \frac{a_{0} a_{3}}{a_{1}} + a_{2} > 0\]

Example 2

\[ \frac{1000}{s^3 + 10s^2 + 31s + 1030}\]
table = make_routh_table(1000 / (s**3 + 10*s**2 + 31*s + 1030))
print_table(table)
evaluate_table(table)
\[\begin{split}\displaystyle \left[\begin{matrix}1 & 31 & 0\\1 & 103 & 0\\-72 & 0 & 0\\103 & 0 & 0\end{matrix}\right]\end{split}\]
2 sign changes.

['+', '+', '-', '+']

Example 3

\[G(s)= \frac{10}{s^5+2s^4+3s^3+6s^2+5s+3}\]
table = make_routh_table(10 / (s**5 + 2*s**4 + 3*s**3 + 6*s**2 + 5*s + 3))
print_table(table)
evaluate_table(table)
\[\begin{split}\displaystyle \left[\begin{matrix}1 & 3 & 5\\2 & 6 & 3\\\epsilon & \frac{7}{2} & 0\\6 - \frac{7}{\epsilon} & 3 & 0\\\frac{- 6 \epsilon \left(\epsilon - 7\right) - 49}{2 \left(6 \epsilon - 7\right)} & 0 & 0\\3 & 0 & 0\end{matrix}\right]\end{split}\]
2 sign changes.

['+', '+', '+', '-', '+', '+']

Example 4

\[ \frac{128}{s^8 + 3s^7 + 10s^6 + 24s^5 + 48s^4 + 96s^3 + 128s^2 + 192s + 128}\]
table = make_routh_table(128 / (s**8 + 3*s**7 + 10*s**6 + 24*s**5 + 48*s**4 + 96*s**3 + 128*s**2 + 192*s + 128))
print_table(table)
evaluate_table(table)
\[\begin{split}\displaystyle \left[\begin{matrix}1 & 10 & 48 & 128 & 128\\1 & 8 & 32 & 64 & 0\\1 & 8 & 32 & 64 & 0\\3 & 16 & 32 & 0 & 0\\1 & 8 & 24 & 0 & 0\\-1 & -5 & 0 & 0 & 0\\1 & 8 & 0 & 0 & 0\\3 & 0 & 0 & 0 & 0\\8 & 0 & 0 & 0 & 0\end{matrix}\right]\end{split}\]
2 sign changes.

['+', '+', '+', '+', '+', '-', '+', '+', '+']

Example 5

\[ \frac{k}{s^3 + 10s^2 + 31s + 30+k}\]

Find the largest \(k\) that does not cause the closed-loop system unstable.

table = make_routh_table(k / (s**3 + 10*s**2 + 31*s + 30+k))
print_table(table)
\[\begin{split}\displaystyle \left[\begin{matrix}1 & 31 & 0\\10 & k + 30 & 0\\28 - \frac{k}{10} & 0 & 0\\k + 30 & 0 & 0\end{matrix}\right]\end{split}\]

Our main interest lies in row #3, column #1: it must be a positive number.

table[2,0]
\[\displaystyle 28 - \frac{k}{10}\]

Hence, we create the inequality and solve for \(k\):

reduce_inequalities(table[2,0] > 0, k)
\[\displaystyle -\infty < k \wedge k < 280\]

Good readings