Control Systems Stability (Routh-Hurwitz Criterion)
01-16-2016, 07:31 AM (This post was last modified: 01-16-2016 07:40 AM by toshk.)
Post: #1
 toshk Member Posts: 195 Joined: Feb 2015
Control Systems Stability (Routh-Hurwitz Criterion)
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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)