Control System Programs
06-28-2020, 01:19 AM (This post was last modified: 08-10-2020 01:55 PM by Major.)
Post: #1
 Major Junior Member Posts: 21 Joined: Jan 2020
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
08-20-2020, 08:31 PM
Post: #2
 victorvbc Member Posts: 70 Joined: Jun 2019
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!
08-21-2020, 12:30 AM
Post: #3
 Major Junior Member Posts: 21 Joined: Jan 2020
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
08-26-2020, 12:13 AM
Post: #4
 dah145 Junior Member Posts: 19 Joined: Aug 2020
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.
11-12-2020, 05:49 PM
Post: #5
 victorvbc Member Posts: 70 Joined: Jun 2019
RE: Control System Programs
loui_chakra Wrote:Hey Victor,

I saw your message that you can share the control system programs. Can you kindly share them?

Cheers,
Louis

Hi!
Apparently you have private messaging disabled, so I'll just send it over here.

I'll send you the backup file from my calculator, as it already has all the CAS functions compiled.

Here's a walkthrough video that I recorded to show you how to use some of the functions:
https://youtu.be/0nizY3YbHDE

Cheers!

Attached File(s)