HP Forums

Full Version: Control Systems Stability (Routh-Hurwitz Criterion)
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
CAS code: first time: compile before use

eg1. how([1,4,3,2,1,4,4])--> x^6+4*x^5+3*x^4+2*x^3+x^2+4*x+4 is stored in nx after running eg1.

eg2. how([1,4,k,2,1,4,4])--> x^6+k*x^4+4*x^5+2*x^3+x^2+4*x+4

The characteristic equation s^4 + 2s^3 + 3s^2 + 4s + K+2 enter as how([1,2,3,4,k+2]).
This code will give you the Routh-Hurwitz Matrix = M25;

This code will not tell you how many pole(s) lies in the left/right planes:
To find how many Sign changes just use the Hp Prime standard function; sturmab(nx,x,-20,0) as in eg1. but not in eg2.(because of the variable k)

PHP Code:
#pragma mode( separator(.,;) integer(h32) )
exc();
dett();
auxx();
#cas
how(den):=
BEGIN
local c
,a,j,b,n,nl,M,rw;
n:=length(den);
M:=length(den);
nx:=simplify(poly2symb(den));
IF 
even(n-1THEN rw:=((n-1)/2)+1; ELSE rw:=(n)/2 END;
M25:=makemat(n,rw);
c:=0;
WHILE 
n>DO
a:=1;L1:={};c:=c+1;
IF 
n==0 THEN       
L1
(a):=coeff(nx,n);      
ELSE
nl:=n-1;
FOR 
j FROM c TO length(denSTEP 2 DO
L1(a):=coeff(nx,x,nl);
a:=a+1;
nl:=nl-2;
END;         
END;
L2(n):=length(L1);
FOR 
j FROM 1 TO length(L1) DO
M25[c,j]:=L1(j);
END;
n:=n-1;
END;
exc();
END;
#end


#cas  
exc()
BEGIN
local j
,u,C,M81,B;
C:=M-1;
FOR 
j FROM 1 TO M-DO
u:=j+1;
B:=j;
M81:=subMat(M25,j,1..rw,u,1..rw);
dett();     
END;
END;
#end

#cas 
dett()
BEGIN
local j
,A,M99,M91;
M99:=matrix(1,rw);
M91:=matrix(2,2);
B:=B+2;
A:=1;
C:=C-1;
FOR 
j FROM 2 TO rw DO
M91[1,1]:=M81[1,1];
M91[2,1]:=M81[2,1];
M91[1,2]:=M81[1,j];
M91[2,2]:=M81[2,j];           
M25[B,A]:=(-DET(M91)/M81[2,1]);  
A:=A+1;
END;
IF 
M25({B})==M99({1}) THEN auxx(); END
IF 
M25[B,1]==0 THEN M25[B,1]:=ε END
END;
#end


#cas 
auxx()
BEGIN
local j
,D;
D:=C;
M99:=M25(B-1);
FOR 
j FROM 1 TO rw DO
IF 
j==1 THEN L3(j):=M99[j]*x^(D); ELSE L3(j):=M99[j]*x^(DEND;
D:=D-2;
END;
M100:=coeff(diff(sum(L3)));
FOR 
j FROM 1 TO length(M100) DO
M25[B,j]:=M100(j);
END;
END;
#end 
Oh this looks awesome! I will try it out. Thanks for adding more Control System programs to the HP Prime. I can't get enough of them. Thanks so much.
Reference URL's