3. Why do we need the Routh Array

In a previous notebook we showed that we can calculate the roots of the denominator of a closed loop transfer function to determine stability regions as a function of \(K_c\). However, it became clear that analytic calculation of the roots would only work for lower-order systems.

Using numeric methods seemed to work OK, but involved trial-and-error.

Numeric root finding algorithms are also problematic. Consider finding the roots of \((1 + s)^{10}\). We can see that they should all be -1. Let’s see how well numpy.roots does in finding them.

import numpy

numpy.roots([1, 2, 1])

array([-1., -1.])

polynomial = [1]
term = [1, 1]

for i in range(10):
    polynomial = numpy.convolve(polynomial, term)

polynomial

array([  1,  10,  45, 120, 210, 252, 210, 120,  45,  10,   1])

roots = numpy.roots(polynomial)
roots

array([-1.0486659 +0.01614412j, -1.0486659 -0.01614412j,
       -1.02925286+0.04166079j, -1.02925286-0.04166079j,
       -0.99899397+0.05030124j, -0.99899397-0.05030124j,
       -0.9701264 +0.03974754j, -0.9701264 -0.03974754j,
       -0.95296087+0.01496287j, -0.95296087-0.01496287j])

We’re making up to 5% error and reporting non-negligable imaginary components, when we know the roots are actually real. So it’s not that easy to make a call about the nature of the roots of high order polynomials by calculating them numerically. And it’s not just because the algorithm isn’t good enough. Evaluating one of the roots gives zero to many decimals. The problem is that computers use finite representations of these numbers.

numpy.polyval(polynomial, roots[0])

(-1.0769163338864018e-13+6.760217385881617e-15j)

4. A better way

The Routh-Hurwitz stabilbility criterion provides an efficient check of stability for closed loop systems which avoids calculating the roots of a higher-order polynomial and is therefore less error prone if we have numeric coefficients and actually possible if we have symbolic coefficients (recall we cannot calculate the roots analytically for orders higher than 4).

import sympy

sympy.init_printing()

s = sympy.Symbol('s')

a_0, a_1, a_2, a_3, a_4 = sympy.symbols('a_0:5')
p = a_0 + a_1*s**1 + a_2*s**2 + a_3*s**3 + a_4*s**4

Note that we have to convert the expression above to a Poly object to recover all the coefficients.

p = sympy.Poly(p, s)
p

$$\operatorname{Poly}{\left( a_{4} s^{4} + a_{3} s^{3} + a_{2} s^{2} + a_{1} s + a_{0}, s, domain=\mathbb{Z}\left[a_{0}, a_{1}, a_{2}, a_{3}, a_{4}\right] \right)}$$

This function constructs the Routh array as given in Seborg.

from tbcontrol.symbolic import routh

help(routh)

Help on function routh in module tbcontrol.symbolic:

routh(p)
    Construct the Routh-Hurwitz array given a polynomial in s

    Input: p - a sympy.Poly object
    Output: The Routh-Hurwitz array as a sympy.Matrix object

routh(p)

$$\left[\begin{matrix}a_{4} & a_{2} & a_{0}\\a_{3} & a_{1} & 0\\- \frac{a_{1} a_{4}}{a_{3}} + a_{2} & a_{0} & 0\\\frac{a_{0} a_{3}^{2} + a_{1} \left(a_{1} a_{4} - a_{2} a_{3}\right)}{a_{1} a_{4} - a_{2} a_{3}} & 0 & 0\\a_{0} & 0 & 0\end{matrix}\right]$$

Let’s try this on example 10.1

K_c = sympy.Symbol('K_c')

ce = 10*s**3 + 17*s**2 + 8*s + 1 + K_c

A = routh(sympy.Poly(ce, s))
A

$$\left[\begin{matrix}10 & 8\\17 & K_{c} + 1\\- \frac{10 K_{c}}{17} + \frac{126}{17} & 0\\K_{c} + 1 & 0\end{matrix}\right]$$

For stability, the left hand column must have entries with all the same signs:

sympy.solve([e > 0 for e in A[:, 0]], K_c)

$$-1 < K_{c} \wedge K_{c} < \frac{63}{5}$$