Post Reply 
Control System Programs
06-28-2020, 01:19 AM (This post was last modified: 08-10-2020 01:55 PM by Major.)
Post: #1
Control System Programs
Here are programs made for a basic Control Systems class.

Summary
1. Controllability (ctrb) and Is controllable? (isctrb)
2. Observability (obsv) and Is observable? (isobsv)
3. Ackermann's Formula (ACKER)
4. Statespace to Transfer Function (ss2tf)
5. Transfer Function to Statespace (tf2ss)
6. Routh-Hurwitz Criterion (rhs)
7. Static Error Constants, Steady State Error, and System Type (sterr)
8. Second Order System Fundamentals (solveFund)
9. Step Function (stp)
10. Step Info (stepInfo)

1. The first program will take as input the system matrix (a) and the input matrix (b) to find the controllability matrix. With it there is another function that will check if the controllability matrix found is controllable or not.

Code:
EXPORT CTRB(a,b)
BEGIN
LOCAL k,Cm,Cmsolve;
IF (rowDim(a) == colDim(a) AND rowDim(a) == rowDim(b) AND SIZE(SIZE(b)) >1) THEN
Cm := MAKEMAT(0,rowDim(a),rowDim(a)*colDim(b));
 FOR k FROM 1 TO rowDim(a) DO
  Cmsolve := a^(k-1)*b;
  Cm := REPLACE(Cm,{1,k*colDim(b)-(colDim(b)-1)},Cmsolve);
 END;
ELSE
 RETURN "Matrix A should be square, and the rows of Matrices A and B should have the same amount of elements.";
END;
RETURN Cm;
END;

EXPORT isctrb(Cm)
BEGIN
 IF (abs(DET(Cm)) >= 1 AND RANK(Cm) == rowDim(Cm)) THEN
  RETURN "Controllable";
 ELSE
  RETURN "Uncontrollable";
 END;
END;

2. The next program will take the system matrix (a) and the output matrix (c) as inputs to find the observability matrix. Along with it is another function to check if the observability matrix is observable.

Code:
EXPORT OBSV(a,c)
BEGIN
LOCAL k,Om,Omsolve;
IF (SIZE(SIZE(c)) == 1) THEN
 c := TRN(list2mat(c,1));
END;
Om := MAKEMAT(0,rowDim(a)*rowDim(c),rowDim(a));
IF (rowDim(a) == colDim(a) AND rowDim(a) == colDim(c)) THEN
 FOR k FROM 1 TO rowDim(a) DO
  Omsolve := c*a^(k-1);
  Om := REPLACE(Om,{k*rowDim(c)-(rowDim(c)-1),1},Omsolve);
 END;
ELSE
 RETURN "Matrix A should be square, and the columns of Matrices A and C should have the same amount of elements.";
END;
RETURN Om;
END;

EXPORT isobsv(Om)
BEGIN
 IF (abs(DET(Om)) >= 1 AND RANK(Om) == rowDim(Om)) THEN
  RETURN "Observable";
 ELSE
  RETURN "Unobservable";
 END;
END;

3. This function below uses Ackermann's Formula to find the feedback constants. This function requires the use of the CTRB(a,b) function. It takes as input the system matrix (a), the input matrix (b), and the desired characteristic equation coefficients.

Changed determinate to determinant

Code:
EXPORT ACKER(a,b,ce)
BEGIN
LOCAL c,cm,sc,sP,cemsolve,p,icm,k;
cm := CTRB(a,b);
 IF (abs(DET(cm)) >= 1 AND RANK(cm) == rowDim(a)) THEN
  sP := SIZE(ce);
  sP := sP(1);
  cemsolve := 0;
  p := revlist(ce);
  sc := SIZE(cm);
  sc := sc(1);
  c := [0];
  FOR W FROM 1 TO sc-2 DO
   c := append(c,0);
  END;
  c := append(c,1);
  FOR Q FROM sP DOWNTO 1 DO
    cemsolve := cemsolve+p(Q)*a^(Q-1);
  END;
  icm := inv(cm);
  k := c*icm*cemsolve;
  RETURN k;
 ELSE
  RETURN "System is uncontrollable, the determinant is zero and the rank does not match the system matrix";
 END;
END;

4. This program converts the system from a statespace model to a transfer function. It takes as arguments the system matrix (a), the input matrix (b), the output matrix (d) and the feedthrough matrix (d).

Code:
#cas
ss2tf(a,b,c,d) :=
BEGIN
LOCAL tfnum,tfden,tfmonic,tftemp,tf,k,q;
IF TYPE(SIZE(c)) == 0 THEN
 c := TRN(list2mat(c,1));
END;
CASE
  IF (TYPE(d) == 0 OR TYPE(d) == 3) THEN d := [[d]] END;
  IF (TYPE(SIZE(d)) == 0) THEN d := TRN(d) END;
  DEFAULT BREAK END;
 IF (rowDim(a) == colDim(a) AND rowDim(b) == colDim(c) AND rowDim(a) == rowDim(b) AND rowDim(a) == colDim(c) AND (SIZE(d) == [1, 1] OR (colDim(d) == colDim(b) AND rowDim(c) == rowDim(d)))) THEN
  tf := (c*(((s*IDENMAT(rowDim(a)))-a)^-1)*b)+d;
  FOR k FROM 1 TO rowDim(d) DO
   FOR q FROM 1 TO colDim(d) DO
    tftemp := tf(k,q);
    tfnum := numer(tftemp);
    tfden := denom(tftemp);
    tfmonic := coeff(tfden,s)(1);
    tfnum := approx((tfnum)/tfmonic);
    tfden := approx((tfden)/tfmonic);
    tfnum := expand(tfnum);
    tfden := expand(tfden);
    tftemp := [tfnum/tfden];
    tf := REPLACE(tf,{k,q},tftemp);
   END;
  END;
  tf := approx(exact(tf));
 ELSE
  RETURN "The number of rows and columns in Matrix A need to match those in B and C respectively. The rows and columns of the D Matrix need to match C and B respectively."
 END;
 IF (SIZE(d) == [1, 1]) THEN
  tf := tf(1,1);
 ELSE
  RETURN tf;
 END;
END;
#end

5. This program converts the system from a transfer function to statespace form. It takes as inputs the numerator and denominator either as a function or as a vector containing the coefficients. The third input (form) is either 0 for Phase-Variable form and 1 for observable canonical form. It only outputs the system matrix (a). To output the other matrices just insert b, c, and d into the CAS console. If you want to check the statespace with MATLAB's take the transpose of the output of form 1.

Code:
#cas
tf2ss(num,den,form) :=
BEGIN
LOCAL A1,A2,Ai2,sA1,Snum,Aden;
CASE
 IF (TYPE(SIZE(num)) == 4 OR TYPE(SIZE(den)) == 4) THEN RETURN "Numerator and Denominator must be vectors"; END;
 IF TYPE(den) == 8 OR TYPE(num) == 8 THEN
  Aden := expand(den);
  Aden := symb2poly(den,s);
  Snum := expand(num);
  Snum := symb2poly(num,s);
  Snum := EXPR(STRING(Snum));
  Aden := EXPR(STRING(Aden));
  
 END;
 IF TYPE(den) == 0 OR 1 OR 4 OR 6 OR TYPE(num) == 0 OR 1 OR 4 OR 6 THEN
  Aden := den;
  Snum := num;
 END;
 DEFAULT RETURN "Args should be poly/(s-p1)*(s-pn) or list/vector";
 END;

 IF (Snum == Aden) THEN RETURN "Numerator is the same as the denominator"; END;

WHILE SIZE(Snum) < SIZE(Aden) DO
  Snum := INSERT(Snum,1,0);
  IF SIZE(Snum) == SIZE(Aden) THEN BREAK; END;
END;

 Snum := approx(Snum./Aden(1));
 Aden := approx(Aden./Aden(1));

 CASE
  IF (form == 0) THEN
   A1 := -(REVERSE(Aden));
   sA1 := size(A1)-1;
   A2 := A1(1..(sA1));
   Ai2 := IDENMAT(sA1);
   a := ADDROW(Ai2,A2,sA1+1);
   a := DELROW(a,1);

   b := MAKEMAT(0,sA1,1);
   b := REPLACE(b,sA1,[1]);

   qrinfo := quorem(Snum,Aden);
   c := qrinfo(2);
   d := TRN(TRN(qrinfo(1)));
   IF (SIZE(qrinfo(2)) != rowDim(a)) THEN
    WHILE SIZE(qrinfo(2)) < rowDim(a) DO
     c := INSERT(c,1,0);
     IF SIZE(c) == rowDim(a) THEN BREAK; END;
    END;
   END;
  c := TRN(TRN(revlist(c)));
  RETURN a;
 END;

 IF (form == 1) THEN
   A1 := -(Aden);
   sA1 := size(A1)-1;
   A2 := A1(2..(sA1+1));
   Ai2 := IDENMAT(sA1);
   a := ADDCOL(Ai2,A2,1);
   a := DELCOL(a,sA1+1);

   qrinfo := quorem(Snum,Aden);
    b := qrinfo(2);
    d := TRN(TRN(qrinfo(1)));
    IF (SIZE(qrinfo(2)) != rowDim(A3)) THEN
     WHILE SIZE(qrinfo(2)) < rowDim(a) DO
      b := INSERT(b,1,0);
      IF SIZE(b) == rowDim(a) THEN BREAK; END;
     END;
    END;
   b := TRN(b);

   c := MAKEMAT(0,1,sA1);
   c := REPLACE(c,1,[1]);

   RETURN a;
  END;
   DEFAULT RETURN "Invalid form, choose 0 for Phase Var. Canonical or 1 for Obsv. Canon. Form"
  END;
END;
#end

6. This program creates the Routh table from the characteristic equation. Insert the characteristic equation as a function of s or as a list or vector of coefficients. You can also enter the system matrix (a) instead of the characteristic equation. This program uses the CAS so if you want to find k just enter k, but you may have to simplify the answer within the table.

Tries to take care of the Epsilon Method (where zero is in first column) and Polynomial Method (Row of Zeros). Does not take into account the reverse coefficient method - this would have to be entered by the user.

There also may be differences as someone may simplify in different ways compared to what is in this table. Also take into consideration if working with numbers if you don't get the correct results in exact mode you may be able to receive the correct results in approximate mode due to simplification.

Updated to account for a characteristic equation with the second element being zero resulting in undefined.

Fixed a two by two routh matrix whose determinant is equal to zero from using row of zeros function instead of epsilon


Code:
#cas
rhs(ce) :=
BEGIN
LOCAL rh,sumsrhb,srhb,ssprh,ct,Srh,nPFlag,scol;
LOCAL csprh,csprhlist,order,q,w,rhdet,sprh,sprhcheck;
LOCAL qw;
CASE
IF TYPE(ce) == 8 AND TYPE(SIZE(ce)) == 0 THEN ce := coeff(ce,s); END;
IF TYPE(ce) == 4 AND TYPE(SIZE(ce)) == 4 THEN ce:= coeff(det(s*IDENMAT(rowDim(ce))-ce),s); END;
IF TYPE(ce) == 4 OR TYPE(ce) == 6 AND TYPE(SIZE(ce)) == 0 THEN ce:=ce END;
DEFAULT RETURN "The characteristic equation needs to either be an equation, list/vector, or a system matrix (a)"; END;
nPFlag := 0;
IF ce(2) == 0 THEN
 ce(2):= ɛ;
END;
IF odd(SIZE(ce)) THEN
 rh := TRN(list2mat(ce,2));
ELSE
 rh := TRN(list2mat(ce,2));
END;
ct := MAKEMAT(0,SIZE(ce)-rowDim(rh),colDim(rh));
FOR q FROM 1 TO rowDim(ct) DO
 rh := ADDROW(rh,ct(q),(1+rowDim(rh)));
END;
Srh := colDim(rh);
FOR q FROM 1 to rowDim(rh) DO
 IF (MIN(rh(q)) != 1 AND MIN(rh(q)) != 0) THEN
  srhb := (rh(q)/content(rh(q)));
  sumsrhb := ΣLIST(srhb);
  IF FP(sumsrhb) == 0 AND FP(srhb(q)) == 0 THEN
   REPLACE(rh,{q,1},srhb);
  END;
 END;
END;
rh := ADDCOL(rh,mat2list(MAKEMAT(0,rowDim(rh),1)),colDim(rh)+1);
FOR q FROM 1 TO SIZE(ce)-2 DO
 ssprh := rh(q+1);
 IF  MIN(ssprh) != 1 AND FP(MIN(ssprh)) == 0 AND SIZE(remove(0,ssprh)) != 1 THEN
  ssprh := ssprh/content(ssprh);
  sprh := ssprh;
  REPLACE(rh,{q+1,1},sprh);
 ELSE
  sprh := ssprh;
 END;
 sprhcheck := rh(q);
 sprhcheck := sprhcheck/content(sprhcheck);
 FOR w FROM 1 TO Srh DO
  rhdet := -det([[rh(q,1),rh(q,w+1)],[rh(q+1,1),rh(q+1,w+1)]])/rh(q+1,1);
  IF nPFlag == 0 THEN
     IF (SIZE(ce) != 4 AND w == 1 AND rhdet == 0 AND (sprh == sprhcheck OR ΣLIST(abs(sprh) .- abs(sprhcheck)) == 0) OR DIFFERENCE(abs(sprh) .- abs(sprhcheck)) == 0) THEN
   REPLACE(rh,{q+1,1},sprh);
   order := rowDim(rh)-2-(q-1);
   csprhlist := {};
   FOR qw FROM 1 TO SIZE(sprh) DO
    csprh := sprh(qw)*s^order;
    csprhlist := CONCAT(csprhlist,csprh);
    order := order-2;
   END;
   csprhlist := diff(csprhlist,s);
   csprhlist := ΣLIST(csprhlist);
   csprhlist := remove(0,coeff(csprhlist,s));
   IF SIZE(csprhlist) != 1 THEN
    csprhlist := csprhlist/content(csprhlist);
   END;
   REPLACE(rh,{q+2,1},csprhlist);
   q := q+1;
   w:= 1;
   nPFlag := 1;
   rhdet := -det([[rh(q,1),rh(q,w+1)],[rh(q+1,1),rh(q+1,w+1)]])/rh(q+1,1);
  END;
  END;
  IF (w == 1 AND rhdet == 0 AND (sprh != sprhcheck OR SIZE(ce) == 4)) THEN
   rhdet := ε;
  END;
  REPLACE(rh,{q+2,w},[rhdet]);
  END;
 END;
 scol := MAKELIST(s^q,q,rowDim(rh)-1,0);
 ADDCOL(rh,scol,1);
 DELCOL(rh,colDim(rh));
RETURN simplify(rh);
END;
#end


7. The program sterr takes the transfer function (xfcn) and input signal (inp) as inputs. Normally the input signal constant will be 1 but you can change this to any number. This program will find the Static Error Constants (Kp - position, Kv - velocity, Ka - acceleration) and also finds their respective steady state errors (essp - step, essv - ramp, essa - parabola). The program also tells you the system type (0 - (s^0 = 1) step, 1 - (s^1 = s) ramp, 2 - (s^2) parabola). Outputs a table containing all the above information.

Code:
#cas
sterr(xfcn,inp) :=
BEGIN
LOCAL kp,essp,kv,essv,ka,essa,ty;
LOCAL sterrlabel,sterrval;
kp := limit(xfcn,s,0);
essp := (inp*1/(1+kp));
kv := limit(s*xfcn,s,0);
essv := (inp*1/kv);
ka := limit(s^2*xfcn,s,0);
essa := (inp*1/ka);
IF ((kv == 0 AND essv == ±∞ OR ∞) AND (ka == 0 AND essa == ±∞ OR ∞)) THEN ty := 0; kp := approx(kp); essp := approx(essp); END;
IF ((kp == ±∞ OR ∞ AND essp == 0) AND (ka == 0 AND essa == ±∞ OR ∞)) THEN ty := 1; kv := approx(kv); essv := approx(essv); END;
IF ((kp == ±∞ OR ∞ AND essp == 0) AND (kv == ±∞ OR ∞ AND essv == 0)) THEN ty := 2; ka := approx(ka); essa := approx(essa); END;
sterrlabel := list2mat({"Static Error Constant","Kp","Kv","Ka","Type"},1);
sterrval := {"Static Error Constant Value",kp,kv,ka,ty};
sterrlabel := ADDCOL(sterrlabel,sterrval,2);
sterrval := {"Static Error", "essp", "essv", "essa","----------"};
sterrlabel := ADDCOL(sterrlabel,sterrval,3);
sterrval := {"----------",essp,essv,essa,"----------"};
sterrlabel := ADDCOL(sterrlabel,sterrval,4);
RETURN sterrlabel;
END;
#end

8. This program tries to solve for all possible input combinations to solve for a second order system. Inputs are s, zeta, wn, perOS, Tr, Tp, and Ts with the leftmost variable taking precedence. It will output s, zeta, wn, perOS, Tr, Tp, Ts, K, and the transfer function.

Please be aware that the equations to solve the system may not exactly adhere to what you may be using. For example, Settling time (Ts) is solved within 2% of the final value. The equation to solve for that can either be -ln(0.02)/zeta*wn or approximately 4/zeta*wn.

Added negative in front of sigma variable in case where Tp and Ts are given

Code:
#cas
solveFund(s,zeta,wn,perOS,Tr,Tp,Ts) :=
BEGIN
LOCAL sigma,wd,funlabel,funval,k,xfcn;
CASE
IF s THEN 
 zeta := MAX(solve(-(RE(s)/IM(s))^2*z^2+(RE(s)/IM(s))^2=z^2,z));
 wn := -RE(s)/zeta;
 perOS := exp((-zeta*3.1415926539)/(sqrt(1-zeta^2)))*100;
 Tr := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/wn;
 Ts := -ln(0.02)/(zeta*wn);
 Tp := 3.1415926539/IM(s);
END;

IF zeta AND Tr THEN
 wn := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/Tr;
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
 Tp := 3.1415926539/(wn*sqrt(1-zeta^2));
 Ts := -ln(0.02)/(zeta*wn);
 perOS := exp((-zeta*3.1415926539)/(sqrt(1-zeta^2)))*100;
END;

IF zeta AND Tp THEN
 wn := 3.1415926539/(Tp*sqrt(1-zeta^2));
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
 Tr := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/wn;
 Ts := -ln(0.02)/(zeta*wn);
 perOS := exp((-zeta*3.1415926539)/(sqrt(1-zeta^2)))*100;
END;

IF zeta AND Ts THEN
 wn := (-ln(0.02*sqrt(1-zeta^2))/(zeta*Ts));
 perOS := exp((-zeta*3.1415926539)/(sqrt(1-zeta^2)))*100;
 Tr := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/wn;
 Tp := 3.1415926539/IM(s);
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
END;

IF zeta AND wn THEN
 Tp := 3.1415926539/(wn*sqrt(1-zeta^2));
 Ts := -ln(0.02)/(zeta*wn);
 Tr := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/wn;
 perOS := exp((-zeta*3.1415926539)/(sqrt(1-zeta^2)))*100;
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
END;

IF wn AND perOS THEN
 zeta := -ln(perOS/100)/(sqrt(3.1415926539^2+(ln(perOS/100))^2));
 Tp := 3.1415926539/(wn*sqrt(1-zeta^2));
 Ts := -ln(0.02)/(zeta*wn);
 Tr := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/wn;
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
END;

IF wn AND Tr THEN
 zeta := MAX(RE(solve((1.76*z^3-0.417*z^2+1.039*z+1)/Tr=wn,z)));
 Tp := 3.1415926539/(wn*sqrt(1-zeta^2));
 Ts := -ln(0.02)/(zeta*wn);
 perOS := exp((-zeta*3.1415926539)/(sqrt(1-zeta^2)))*100;
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
END;

IF wn AND Tp THEN
 zeta := MAX(fsolve((3.1415926539/(Tp*sqrt(1-z^2)))=wn,z));
 Tr := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/wn;
 Ts := -ln(0.02)/(zeta*wn);
 perOS := exp((-zeta*3.1415926539)/(sqrt(1-zeta^2)))*100;
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
END;

IF wn AND Ts THEN
 zeta := MAX(fsolve((-ln(0.02*sqrt(1-zeta^2))/(z*Ts))=wn,z));
 Tr := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/wn;
 Tp := 3.1415926539/(wn*sqrt(1-zeta^2));
 perOS := exp((-zeta*3.1415926539)/(sqrt(1-zeta^2)))*100;
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
END;

IF Tr AND perOS THEN
 zeta := -ln(perOS/100)/(sqrt(3.1415926539^2+(ln(perOS/100))^2));
 wn := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/Tr;
 Tp := 3.1415926539/(wn*sqrt(1-zeta^2));
 Ts := -ln(0.02)/(zeta*wn);
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
END;

IF Tr AND Tp THEN
 zeta := MIN(fsolve((1.76*z^3-0.417*z^2+1.039*z+1)/Tr=3.1415926539/(Tp*sqrt(1-z^2)),z));
 wn := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/Tr;
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
 Ts := -ln(0.02)/(zeta*wn);
 perOS := exp((-zeta*3.1415926539)/(sqrt(1-zeta^2)))*100;
END;

IF Tr AND Ts THEN
 zeta := MAX(RE(fsolve((1.76*z^3-0.417*z^2+1.039*z+1)/Tr=(-ln(0.02)/(z*Ts)),z)));
 wn := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/Tr;
 Tp := 3.1415926539/(wn*sqrt(1-zeta^2));
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
 perOS := exp((-zeta*3.1415926539)/(sqrt(1-zeta^2)))*100;
END;

IF Tp AND perOS THEN
 zeta := -ln(perOS/100)/(sqrt(3.1415926539^2+(ln(perOS/100))^2));
 wn := 3.1415926539/(Tp*sqrt(1-zeta^2));
 Tr := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/wn;
 Ts := -ln(0.02)/(zeta*wn);
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
END;

IF Tp AND Ts THEN
 sigma := -ln(0.02)/Ts;
 wd := 3.1415926539/Tp;
 s := -sigma+wd*i;
 zeta := MAX(solve(-(RE(s)/IM(s))^2*z^2+(RE(s)/IM(s))^2=z^2,z));
 wn := 3.1415926539/(Tp*sqrt(1-zeta^2));
 Tr := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/wn;
 perOS := exp((-zeta*3.1415926539)/(sqrt(1-zeta^2)))*100;
END;

IF Ts AND perOS THEN
 zeta := -ln(perOS/100)/(sqrt(3.1415926539^2+(ln(perOS/100))^2));
 wn := -ln(0.02)/(zeta*Ts);
 Tr := (1.76*zeta^3-0.417*zeta^2+1.039*zeta+1)/wn;
 Ts := -ln(0.02)/(zeta*wn);
 s := -zeta*wn+(wn*sqrt(1-zeta^2))*i;
END;

DEFAULT RETURN "This combination of inputs does not provide enough information!"; END;

k := wn^2;
xfcn := 0;
funlabel := list2mat({"S (σ+*ωd)", "Zeta (ζ)", "Natural Frequency (ωn)", "Percent Overshoot (%OS)", "Rise Time (Tr)", "Peak Time (Tp)", "Settling Time (Ts)"},1);
funval := {s,zeta,wn + " rad/s",perOS + " %",Tr + " s",Tp + " s",Ts + " s"};
funlabel := ADDCOL(funlabel,funval,2);
purge(s);
xfcn := k/(s^2+2*zeta*wn*s+k);
REPLACE(funlabel,{9,2},[xfcn]);
RETURN funlabel;

END;
#end

9. The stp (step) program graphs the step response of the transfer function. It takes as input the transfer function, gain, and a boolean for if the system is an open loop or closed loop unity feedback system.

*This program is a bit unstable and may cause your calculator to crash* Reasons for this may be due to the complexities of the transfer function or errors within or when running the program.
HP Prime G1 users may have troubles with underdamped systems. Unfortunately, I cannot find a way to always distinguish between infinity and undefined.
*HP Prime G2 users may not always experience the same issue.*


The program outputs the step response plot in F2, the upper and lower rise time F4 and F5, the upper and lower settling times (within 2% of final value) F7 and F8, and the final value F0. These are provided in case stepInfo is innacurate, verify, or prefer to find the values manually.

Code:
#cas
stp(xfcn,k,typesys) :=
BEGIN
LOCAL Gfn,Hfn,usrg,dxfcn,nxfcn,roots,ceroots,minev;
LOCAL Tfs,y,Tfinal;
  F2 := "";
  F4 := "";
  F5 := "";
  F7 := "";
  F8 := "";
  F0 := "";
Gfn := xfcn;
Hfn := 1;
IF (k <= 0) THEN
 RETURN "Gain (K) must be greater than or equal to one (1). If gain is already taken into account in the transfer function then leave as one (1)";
END;
CASE
IF (typesys == 0) THEN
 solveTf := (k*Gfn);
END;
IF (typesys == 1) THEN
 solveTf := ((k*Gfn)/(1+Hfn*k*Gfn));
END;
DEFAULT RETURN "Enter either zero (0) for open feedback or one (1) for closed (Unity) feedback system";
END;
usrg := solveTf*(1/s);
usrg := collect(approx(usrg));
usrg := ilaplace(usrg,s,t);
dxfcn := collect(denom(solveTf));
dxfcn := cZeros(dxfcn,s);
nxfcn := collect(numer(solveTf));
nxfcn := cZeros(dxfcn,s);
roots := EVAL(mat2list(dxfcn+k*nxfcn));
ceroots := EVAL(mat2list(dxfcn));
minev := MIN(abs(RE(ceroots)));
 IF (minev == 0) THEN
  Tfinal := 10;
 ELSE
  Tfs := 7/minev;
  y := 10^(CEILING(log10(Tfs))-1);
  Tfinal := y*CEILING(Tfs/y);
 END;
ese := string(usrg);
ese := REPLACE(ese,"t","X");
ese := cat("(", ese, ")", "*Heaviside(X)");
F2 := ese;
Yfinal := approx(RE(F2(999999999999999999999999999999)));
IF (Yfinal > 999999999 OR Yfinal < -999999999) THEN
  BREAK;
ELSE
  F4 := STRING(Yfinal*0.9);
  F5 := STRING(Yfinal*0.1);
  F7 := STRING(Yfinal-Yfinal*0.02);
  F8 := STRING(Yfinal+Yfinal*0.02);
  F0 := STRING(Yfinal);
END;
CHECK(2);
Xmin := 0;
Xmax := approx(Tfinal);
Ymin := 0;
IF (Yfinal > 999999999 OR Yfinal < -999999999) THEN
 Ymax := 999;
ELSE
 Ymax := approx(1.10*Yfinal);
END;
IF (Yfinal == 0) THEN
 Ymax := 10;
END;
RETURN "Press the Plot button to view the graph";
END;
#end

10. The step info (stepInfo) program tries to find the step response characteristics of the step response graph. It has no input, and outputs a table containing the Rise Time (Tr), Peak Time (Tp), Peak Value, Settling Time, Settling Value, Percent Overshoot, and Final Value.

This program may take a while to run, and may be unstable. From my experience it seems to run faster when "Use i to factor polynomials" or the "Use i:" CAS setting is unchecked.

This program requires the Step function (stp) to be run first.


If it returns error or incorrect results you can try repositioning the step response plot and running it again.

Code:
#cas
stepInfo() :=
BEGIN
LOCAL Yfinal,Tfinal,Tr,Tp,Ts,perOS;
LOCAL ess,Trupper,Trlower,Tsu,Tsl,TsTwoMax,TsTwoMin;
LOCAL steplabel,stepval,stepinfo,xsteplist,ysteplist,dt;
Yfinal := approx(RE(F2(999999999999999999999999999999)));
Tfinal := Xmax;
IF (Yfinal > 999999999 OR Yfinal < -999999999) THEN
  Tr := ∞;
  Tp := ∞;
  Tpvalue := ∞;
  Ts := ∞;
  Tsvalue := ∞;
  perOS := ∞;
  Yfinal := ∞;
  ess := ∞;
ELSE
 ess := (1-Yfinal)*100;
 xsteplist := {};
 ysteplist := {};
dt := 0.1;
 FOR q FROM 0 TO Tfinal STEP dt DO
  xsteplist := CONCAT(xsteplist,q);
  ysteplist := CONCAT(ysteplist,F2(q));
 END;
 Tp := MAX(ysteplist);
 Tpvalue := POS(ysteplist,Tp);
 Tpvalue := approx(xsteplist(Tpvalue));
 Tpvalue := EXTREMUM(F2(X),Tpvalue);
 IF (Tp <= Yfinal) THEN
  Tp := Yfinal;
 ELSE
   Tp := subst(F2(x),x = Tpvalue);
 END;
 Trupper := MIN(solve(F2(X)=Yfinal*0.9,X=0..Tpvalue));
 Trlower := MIN(solve(F2(X)=Yfinal*0.1,X=0..Tpvalue));
 Tr := Trupper - Trlower;
 Tsu := (Yfinal+(Yfinal*0.02));
 Tsl := (Yfinal-(Yfinal*0.02));
 TsTwoMax := solve(F2(X)=Tsu,X = 0 .. Tfinal);
 TsTwoMin := solve(F2(X)=Tsl,X = 0 .. Tfinal);
 Ts := MAX(CONCAT(TsTwoMin, TsTwoMax));
 Tsvalue := F2(Ts);
 IF (Tp > Yfinal) THEN
  perOS := ((Tp-Yfinal)/(Yfinal))*100;
 ELSE
  perOS := 0;
 END;
END;
steplabel := list2mat(["Rise Time (Tr)","Peak Time (Tp)", "Peak Value","Settling Time (Ts)", "Settling Value","Percent Overshoot","Final Value"],1);
stepval := [Tr,Tpvalue,Tp,Ts,Tsvalue,perOS,Yfinal];
IF (Tp <= Yfinal) THEN
 REPLACE(stepval,2,[∞]);
END;
stepinfo := ADDCOL(steplabel,stepval,2);
RETURN stepinfo;
END;
#end

Control System Programs
Find all posts by this user
Quote this message in a reply
08-20-2020, 08:31 PM
Post: #2
RE: Control System Programs
Pretty cool! I think the Prime is a very good environment to solve Ctrl Systems problems outside of Matlab.

Ever since getting the Prime I've built my own library myself, but never thought to share it because commenting and documenting it would be a bit of a pain. A lot of those programs were written directly on the calculator as well, which doesn't make for the prettiest presentation.

They include: ctrb, obsv, tf2ss, ss2tf, place (acker), LQR/LQE solver, H-infinity norm (bisection), Disk Margin, ctrb/obsv discrete/continuous Gramians, Lyapunov LME solver, minimal SS representation of transfer matrices, Kalman Decomposition, Similarity Transformations to Canonical Forms, ZOH discretization of SS model, discrete/continuous SS simulation and plot, ZOH/Tustin/Prewarp/Euler TF discretization, RLocus/Bode/Nyquist/Nichols/Step plots for both s and z domains, Routh criterion and Modified Routh criterion, and many others...

Let me know if you want me to share some of those.

Cheers!
Find all posts by this user
Quote this message in a reply
08-21-2020, 12:30 AM
Post: #3
RE: Control System Programs
(08-20-2020 08:31 PM)victorvbc Wrote:  Pretty cool! I think the Prime is a very good environment to solve Ctrl Systems problems outside of Matlab.

Ever since getting the Prime I've built my own library myself, but never thought to share it because commenting and documenting it would be a bit of a pain. A lot of those programs were written directly on the calculator as well, which doesn't make for the prettiest presentation.

They include: ctrb, obsv, tf2ss, ss2tf, place (acker), LQR/LQE solver, H-infinity norm (bisection), Disk Margin, ctrb/obsv discrete/continuous Gramians, Lyapunov LME solver, minimal SS representation of transfer matrices, Kalman Decomposition, Similarity Transformations to Canonical Forms, ZOH discretization of SS model, discrete/continuous SS simulation and plot, ZOH/Tustin/Prewarp/Euler TF discretization, RLocus/Bode/Nyquist/Nichols/Step plots for both s and z domains, Routh criterion and Modified Routh criterion, and many others...

Let me know if you want me to share some of those.

Cheers!

Holy moly you are awesome. Wanted to do more but I don't have the knowledge and still have some bugs with some of mine. Really just did this since I liked my Controls Class so much. A lot of the stuff you have I haven't heard of; we only just started some of the discrete (z domain) stuff at the end of the semester.

I mean it's up to you if you want to share them, but I'd definitely be interested in them all if you do! Mainly just to play around.

Control System Programs
Find all posts by this user
Quote this message in a reply
08-26-2020, 12:13 AM
Post: #4
RE: Control System Programs
(08-20-2020 08:31 PM)victorvbc Wrote:  Pretty cool! I think the Prime is a very good environment to solve Ctrl Systems problems outside of Matlab.

Ever since getting the Prime I've built my own library myself, but never thought to share it because commenting and documenting it would be a bit of a pain. A lot of those programs were written directly on the calculator as well, which doesn't make for the prettiest presentation.

They include: ctrb, obsv, tf2ss, ss2tf, place (acker), LQR/LQE solver, H-infinity norm (bisection), Disk Margin, ctrb/obsv discrete/continuous Gramians, Lyapunov LME solver, minimal SS representation of transfer matrices, Kalman Decomposition, Similarity Transformations to Canonical Forms, ZOH discretization of SS model, discrete/continuous SS simulation and plot, ZOH/Tustin/Prewarp/Euler TF discretization, RLocus/Bode/Nyquist/Nichols/Step plots for both s and z domains, Routh criterion and Modified Routh criterion, and many others...

Let me know if you want me to share some of those.

Cheers!


Would really apreciate if you could share all of them too.
Find all posts by this user
Quote this message in a reply
Post Reply 




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