ryoder@ecn.purdue.edu (Robert W Yoder) (07/02/90)
The following is the response I received from Paul Dujmich after I sent him a
Routh-Hurwitz program I had written. I had never posted it because I had not
gotten around to adding the necessary code to cover the special cases. He
expanded it to cover those cases, so here it is.
Robert Yoder
ryoder@ecn.purdue.edu
--------------------------------------------------------------------------------
Bob
I took the liberty of modifying your basic Routh program for the purpose of
adding special case checking. Most of the special case code is in a separate
subroutine I call "ARYCHK". It was just easier to do it this way.
The new code checks for left-most column zero, as well as zeroes all the
way across a row (premature termination). I believe these two are the
only special cases to worry about. I have tested the new version on my
homework problems, which I had worked out by hand. They included all the
special case conditions mentioned above. The program produced accurate
results for every polynomial that was tried.
One note is in order.
For left-most zero, we use the epsilon method, replacing the 0 with epsilon,
a very small number approaching zero.Then we take the limit of the term,
as epsilon approaches 0. In the program, I use +/- 1E-50 to represent
epsilon. This works fine, but produces slightly different numbers in the
finished array (compared to the limit method). But, we really don't care
about the actual values anyway, all we care about is the sign changes, and
the program has been correct in every case tested. For zeroes all the way
across a row, I do the accepted method of forming a new Aux polynomial
from the coefficients in the row just above the row of zeroes. Then I take
the derivative of the Aux polynomial, and use the new coefficients to
replace the row of zeroes. This has also worked flawlessly for every thing
that was tried.
Please feel free to post this version on net.handhelds if you want. I'm sure
there are other EE/BSET people who would like to have it. Thanks again for
your reply to my posting, and for your original program. It sure was a lot
easier not having to start from scratch. Well anyway, here's the updated
version.
Thanks Again
Paul
pauld@fs1.ece.cmu.edu
============================================================================
%%HP: T(3)A(R)F(.);
@ PROGRAM NAME: ROUTH
@ DATE: 6-30-90
@ PURPOSE: Does Routh-Hurwitz Stability Testing
@ ARGUMENTS: 1. An array in level 2 containing coefficients of the
@ characteristic polynomial from the transfer function.
@ 2. The ORDER of the characteristic polynomial in level 1
@ NOTE: This program requires "ARYCHK", a subroutine that checks for special
@ cases, to operate.
@Routh-Hurwitz Program
@This program tests the characteristic polynomial of a transfer function
@for stability. The program expects an array as input, which is made up of
@the coefficients of the characteristic polynomial. The array must be in stack
@level 2 when the program is executed. The program also expects the ORDER of
@the characteristic polynomial in level 1.
@The program outputs a Routh-Hurwitz array as it's output.
@The stability of the system being tested can be obtained
@by examining the left-most column of the output array. If all elements in the
@left-most column of the array are of the same sign, the system is STABLE, and
@all of the polynomial roots are in the left-half plane of the complex plane.
@If there are any sign changes within these elements, then the system is
@UNSTABLE, and the polynomial has roots in the RIGHT HALF PLANE of the complex
@plane, or
@it has roots on the IMAGINARY AXIS of the complex plane. The NUMBER of sign
@changes within the left-most array elements give the number of roots in the
@right half plane. The program tests each input array for the following
@SPECIAL CASES, and applies a fix, so that the array may be completed:
@ Special Case 1. Zero in left-most array element
@ Special Case 2. Zeroes all the way across a row
\-> ord @load in polynomial ORDER
\<< DUP @make copy of input array on stack
SIZE @return dimensions of input array {2 2} on stack
LIST\-> @dump the array dimensions from list to stack
@ 2 (rows)
@ 2 (columns)
@ 2 (there were 2 elements in the list)
DROP @throw away the #elements in list, don't need it
'cols' STO @save # columns of input array
DROP @drop off # columns from stack
cols 3 * @new # of rows = 3* original # columns
cols @orig number of columns
2 \->LIST @put last 2 elements in a list
RDM @redimension input array to size specified by list
'a' STO @save the redimensioned array
3 'r' STO @start calculating at row 3
1 CF @clear flag 1 in case it is set
@**************** calculate # of coefficients in orig array *************
ord 1 + @increment polynomial order to get # coefficients
'coef' STO @save it
@************** Check row 2 of input array for special cases ***************
CLEAR @clear stack
3 'r' STO
ARYCHK @check row 2 of the input array for special cases
@if they exist-fix them.
@********************** main DO LOOP ****************************************
3 'r' STO @start calculating at row 3
DO
1 cols 1 - @set up FOR loop for 1 less than number of columns
FOR c
a @put array on stack
r @put current row on stack
c @put current column within current row on stack
2 \->LIST @make a list of last 2
@************* array element calculation ************************
'(a(r-1,1)*a(r-2,c+1)-a(r-2,1)*a(r-1,c+1))/a(r-1,1)'
EVAL
PUT @put the calculated element back into the array
'a' STO @save the updated array
NEXT
1 'r' STO+ @increment to next row
ARYCHK @check for special cases
UNTIL 1 FS? @run DO loop until flag 1 is set
END
a @ put array on stack
r 1 - cols 2 \->LIST
RDM @ redimension the array
{ a cols r s cnt coef hc pwr} @list of variables to nuke
PURGE @nuke them
1 CF @clear flag 1
\>>
=============================================================================
%%HP: T(3)A(R)F(.);
@ PROGRAM NAME: ARYCHK (subroutine for Routh-Hurwitz Program)
@ DATE: 6-30-90
@ PURPOSE: checks a row of a Routh-Hurwitz array for special cases
@ARGUMENTS: Variable 'r' must exist and contain the row of the Routh-
@Hurwitz array (+1) to be checked.
\<<
IF 'a(r-1,1)==0' @if leftmost element of last row == 0
THEN
IF 'r-1 > coef' @normal valid array termination
THEN 1 SF @set the flag
@special case checking
ELSE 0 'cnt' STO @clear cnt
1 cols @FOR loop goes from 1 to
@ # columns
FOR c
'a(r-1,c)' EVAL @count # valid array
@elements in row where 0
@occurred
0 \=/ @element !=0?
'cnt' STO+ @if yes, then valid
NEXT
IF 'cnt==0' @the whole row had 0's
@premature termination of array (all 0's in row)
THEN
coef r 2 - - @calculate first power of
@S for new polynomial
'pwr' STO @save it
1 cols @FOR loop goes from 1 to
@width of a column
FOR c
'a(r-2,c)'
EVAL @get term coefficient
"*S^"
+ @add it to string
pwr @get the exponent
+ @add it to string
"'"
SWAP
+ @make algebraic object
@ inside a string
OBJ\-> @remove string quotes
@object is now algebraic
1 'S' STO
'S' @variable of differentiation
\.d @take the derivative
EVAL @evaluate at S = 1
a SWAP @put array on stack
r 1 - c 2 \->LIST @elem to change
SWAP
PUT @correct the element
'a' STO @save updated array
pwr 2 - @decrement exponent by 2
'pwr' STO @save new power
IF pwr 0 < @if exponent < 0 then
@terminate the FOR loop
THEN cols 1 +
'c' STO
END
NEXT
'S' PURGE @ done with S for now
@zero in left-most column only
ELSE 0 'cnt' STO @clear cnt
1 coef @FOR loop goes from 1 to
@ # of coefficients
FOR c
'a(c,1)' @count the number of left
EVAL @most elements that are
0 < @negative
'cnt' STO+ @cnt has total
NEXT
coef 2 / @find half #coefficients
IP @only the integer part
'hc' STO @save it in hc
IF 'cnt > hc' @more than 1/2 the
@coefficients are negative
THEN a @get array
r 1 - 1 2
\->LIST @array element to change
-1.E-50 @make it neg
PUT @put in array
'a' STO @save changed array
ELSE a @get array
r 1 - 1 2 @select element
\->LIST
1.E-50 @make it pos
PUT @put in array
'a' STO @save changed array
END
END
END
END
\>>
==============================================================================
--
--------------------------------------------------------------------------------
Robert Yoder
Internet: ryoder@ecn.purdue.edu