# 63. 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 = 
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)

:

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


# 64. 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}$$