Last week's tutorial for a control systems unit I'm taking at the moment consisted almost entirely of practising the Routh-Hurwitz method for determining polynomial stability. For anyone who's ever had to spend time constructing Routh tables, you'll know just how tedious it is. Which is why during that tute I spent my time much more profitably writing this small python program for calculating Routh tables.
Of course, there are already plenty of scripts/programs on the internets that can do this (eg this one, for Matlab/Octave). So why bother writing another one?
Well, firstly for my own education and entertainment. But also because the script linked above has many deficiencies. It is completely unable to deal with the case of symetrically distributed roots or marginal stability. Its method for handling the case of a zero in the first column is extremely kludgy. The correct way to do this is to replace the zero with a small quantity, epsilon, and then consider the limit as epsilon -> 0. The above script simply replaces the zero with the fixed quantity 0.01 - which is fine if your polynomial coefficients are on the order of 1 or higher, but not so fine if you have coefficients on the order of 0.001 or smaller. This is not just me being pedantic, either - these situations occur quite regularly.
Thus I decided to write my own python script. In order to handle all the arithmetic exactly (ie to avoid problems with floating point round-off) I decided to use the sympy library for python. You will need to have this library installed to run my script (on debian et al, it's just a matter of sudo apt-get install python-sympy). Using sympy has a further advantage - it allows symbolic variables to be entered as polynomial coefficients.
There are several classes of (useful) problem which are best solved through this approach. Eg this one
Q. A control system has an open loop transfer function of k/(s^3 + 3s^2 + s). For what range of gains is the unity feedback system stable?
You can solve this with my script easily enough. First deduce the closed loop transfer function (ie k/(s^2 + 3s^2 + s + k)). Enter the coefficients of the denominator polynomial when prompted (ie 1, 3, 1, k). This produces the output
1 1
3 k
1-k/3
k
For stability, all the numbers in the first column must be positive. Thus from the Routh table, we can conclude k>0 and 1-k/3 > 0 => 0 < k < 3
My script also works in the case where a polynomial has roots distributed symmetrically in the complex plane. This condition is indicated by a row in the Routh table consisting entirely of zeros. In these cases, the auxiliary polynomial is calculated, and the row of zeros is replaced with the derivative of the auxiliary polynomial. The program displays the auxiliary polynomial on-screen, allowing the user to check for the case of marginal stability (ie purely imaginary roots).
Here is the script:
http://dl.dropbox.com/u/1614464/routhHurwitz/routhHurwitz.py
