Control Systems Stability (Routh-Hurwitz Criterion)

[toshk]

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)

Code:


#pragma mode( separator(.,;) integer(h32) )
exc();
dett();
auxx();
#cas
how(den):=
BEGIN
local c,a,j,b,n,nl;
n:=length(den);
M:=length(den);
nx:=simplify(poly2symb(den));
IF even(n-1) THEN rw:=((n-1)/2)+1; ELSE rw:=(n)/2 END;
M25:=makemat(n,rw);
c:=0;
WHILE n>0 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(den) STEP 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();
return "Routh-Hurwitz Matrix"=M25
END;
#end


#cas  
exc()
BEGIN
local j,u;
C:=M-1;
FOR j FROM 1 TO M-2 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:=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^(D) END;
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