/* Maxima Control Engineering Package COMA ("COntrol engineering with MAxima") */ /* Wilhelm Haager, 2019-06-14 (wilhelm.haager@htlstp.ac.at) */ load(coma); /* Plot-Routines*/ /* The list "coma_defaults" holds default options (mainly graphic options) for COMA functions in the style of an associative array: */ coma_defaults; /* Plotting FOUR functions with ONE variable each: */ plot([sin(5*x)**2,0.8*sin(5*y)**2,x,0.8*y],xrange=[-0.5,1.5], color=[red,green,brown,blue],line_type=[solid,solid,dots,dots])$ /* Plotting TWO functions with TWO variables each: */ plot([sin(5*x)**2+0.8*sin(5*y)**2,x+0.8*y], xrange=[-0.5,1.5],surface_hide=true,color=[red,green,brown,blue])$ /* Plotting the damping of a transfer function with respect to the parameter "a": As the function "damping" requires a transfer function with numeric coefficients, it has to be quoted. */ f:(s+a)/(s^3+a*s^2+2*s+a); plot('damping(f),xrange=[0,5]); plot('damping(f),logx=true,logy=true,xrange=[0.05,50]); /* Isolines of the damping of a transfer function with respect to two parameters "a" and "b": */ f:1/(s^5+s^4+6*s^3+a*s^2+b*s+1); contourplot('damping(f),a,b, xrange=[0,10], yrange=[0,10], contours=[0.10,0.05,0,-0.05,-0.1,-0.15], color=[green,green,black,red,red,red])$ /* Electrical Engineering Operators and Functions*/ /* Parallel connection of impedances: */ R1//R2; /* Complex quantities in polar coordinates with the cis-Operator "/_" for input and the postfix operator "//_" as well as the function "cisform" for output: */ U1:230; U2:230 /_ 240; U3:230 /_ 120; cisform(U1-U3); U1-U3 //_ ; /* Data Processing*/ /* Moving Average */ /* List with values of a noisy sine: */ sinli:makelist(float(sin(t)-0.1+random(0.2)),t,-%pi/10,21/10*%pi,%pi/60)$ /* Smoothed list: */ sinlis:moving_average(sinli,15)$ plot([points(sinli),points(sinlis)],point_type=0,points_joined=true, yrange=[-1.1,1.1]); /* Regression Analysis */ /* List of points: */ pli:[[0,0.2],[2,0.5],[4,1.4],[5,2],[6,3],[7,4.5],[8,6.5]]; /* Polynomial Regression: */ p:regan(pli,[1,x,x**2,x**3]); plot([points(pli),p],point_type=7,points_joined=false, xrange=[0,10],yrange=[-1,10]); /* Laplace-Transform, Step Response*/ /* Laplace transform and inverse Laplace transform are part of Maxima: */ laplace(t^2*sin(a*t),t,s); ilt(1/(s^3+2*s^2+2*s+1),s,t); ilt(1/((s+a)^2*(s+b)),s,t); /* Inverse Laplace transform fails at higher order transfer functions: */ ilt(1/(s^3+2*s^2+3*s+1),s,t); /* nilt calculates inverse Laplace transform with numerically calculated poles. As the result is valid only for t>=0. If unit_step_included is set to "true" (default), it is multiplied with the unit step function: */ nilt(1/(s^3+2*s^2+3*s+1),s,t); nilt(1/(s^3+2*s^2+3*s+1),s,t),unit_step_included=false; /* nilt tries to perform the back-transformation analytically, unless the variable try_symbolic_ilt is set to false; the additional parameters for the Laplace variable and the time variable can be omitted (if assumed to be "s" and "t", respectively): */ nilt(1/(s^3+2*s^2+2*s+1)); nilt(1/(s^3+2*s^2+2*s+1)),try_symbolic_ilt=false; /* List of second order systems: */ pt2li:create_list(1/(s^2+2*d*s+1), d, [0.0001,0.1,0.2,0.3,0.4,0.5]); /* Unit step responses of the second order systems: */ step_response(pt2li)$ /* Response on a piecewise defined input function: */ u:5*unit_step(t)-5*unit_step(t-1)+5*unit_step(t-1.5)-5*unit_step(t-2); plot(u,xrange=[-1,5],yrange=[-2,6]); U:laplace(u,t,s); X:U*1/(1+s); x:nilt(X); plot([u,x],xrange=[-1,5],yrange=[-2,6]); response([explicit(u,t,-1,5),X],xrange=[-1,5],yrange=[-2,6]); /* Frequency Response*/ /* Magintude plots of second order systems: */ magnitude_plot(pt2li, yrange=[0.1,10])$ /* Magnitude plots with linear scale: */ magnitude_plot(pt2li,logy=false,yrange=[0,5])$ /* Magnitude plot scaled in dB: */ magnitude_plot(pt2li,logy=dB,yrange=[-20,20]); /* Phase plots of second order systems: */ phase_plot(pt2li)$ /* Bode plot (magnitude plot and phase plot together: */ bode_plot(pt2li,yrange=[[0.1,10],[-180,0]])$ /* Asymptotic characteristics */ f:ratsimp(2*(1+1/(5*s)+s/(1+0.1*s))); bode_plot([f,asymptotic(f)],yrange=[[0.1,100],[-120,120]]); /* Nyquist plots (frequency response locus) of second order systems: */ nyquist_plot(pt2li, xrange=[-3,3], yrange=[-5,1], dimensions=[500,500], nticks=2000)$ /* Nichols plot: */ nichols_plot(pt2li,xrange=[-200,20],yrange=[0.01,10]); /* Calculation of gain and phase in degrees of a frequency response: */ F:5*(1+5*s)/((1+3*s)*(1+s)); absval(F); phase(F),numer; /* Investigations in the s-Plane*/ fli:makelist(stable_rantranf(5),k,1,2); /* Calculating the zeros and poles of a list of transfer functions: */ zeros(fli); poles(fli); /* Zeros and poles of a list of transfer functions in the complex plane: */ poles_and_zeros(fli)$ /* Calculation of closed loop transfer functions: */ fli:closed_loop(makelist(k*((s-a)*(s+1)) /(s*(s-2)*(s+7)),a,-11,-8)); /* Root locus plots with respect to the open loop gain of a list of transfer functions in the complex plane: */ root_locus(fli,xrange=[-17,3], yrange=[-6,6],nticks=5000)$ /* Root locus plots with respect to the damping ratio of a second order of transfer function: */ root_locus(1/(s**2+a*s+1),xrange=[-4,1], trange=[1e-4,3],nticks=5000); /* Transfer Functions:*/ /* Random Transfer Functions */ /* List of transfer functions with random coefficients: */ fli:makelist(tf(random,3),k,1,4); /* Checking the stability of the transfer functions: */ stablep(fli); /* List of stable transfer functions with random coefficients: */ fli:makelist(tf(random,stable,3),k,1,4); /* Checking the stability of the transfer functions: */ stablep(fli); /* Random normalized transfer function: */ tf(random,stable,normalized,7); /* Filters */ /* Butterworth filter: */ poles_and_zeros(tf(butterworth,5)); step_response(makelist(tf(butterworth,k,2),k,3,10),xrange=[0,10],yrange=[0,1.3]); /* Bessel filter */ poles_and_zeros(tf(bessel,5)); step_response(makelist(tf(bessel,k,2),k,3,10),xrange=[0,10],yrange=[0,1.3]); /* Cebychev filter */ poles_and_zeros(tf(chebychev,5)); step_response(makelist(tf(chebychev,k,2),k,3,10),xrange=[0,10],yrange=[0,1.3]); magnitude_plot(makelist(tf(chebychev,k,2),k,3,10),yrange=[0.01,10]); magnitude_plot(makelist(tf(chebychev,6,2,eps),eps,0.1,0.5,0.1),yrange=[0.5,5]); /* PTn with poles on a semicircle (similar to Bessel filter), but poles confined to an angle (e.g.30 degrees): */ poles_and_zeros([tf(ptn,10,1,30), parametric(r*cos(150*%pi/180),r*sin(150*%pi/180),r,0,1.2), parametric(r*cos(-150*%pi/180),r*sin(-150*%pi/180),r,0,1.2)], yrange=[-1,1],xrange=[-2,1],color=[red,gray,gray]); /* Basic control system elements */ [tf(pt1),tf(dt1),tf(pi),tf(pd)]; [tf(pt2),tf(pt2,2,1,5),tf(pt2,prod),tf(pt2,prod,2,3,4)]; [tf(pid),tf(pid,prod)]; [tf(pidt1),tf(pidt1,prod)]; /* List with open loop transfer functions: */ fo:[k/s,5/(s*(s+3)),1-b/s]; /* Calculation of the closed loop transfer functions: */ fw:closed_loop(fo); /* Determining the types of the transfer functions: */ tftype(fo); tftype(fw); open_loop(fw); /* Various transfer functions */ [tf(generic),tf(generic,3,5,a,b)]; /* Transfer function of an impedance chain with series resistances and traverse capacitors: */ impedance_chain(R,1/(s*C),4); /* Fifth order Pade approximation of a time delay: */ time_delay(T,5); /* Transfer function with coefficients optimized according to the ISE-criterion: */ [tf(ise,5),tf(ise,5,2)]; step_response([tf(ise,5),tf(ise,5,2)]); /* Canonical forms of transfer functions */ f:rantranf(8); sum_form(f,1); sum_form(f,2); sum_form(f,3); sum_form(f,4); product_form(f,1); product_form(f,2); /* Stability Behaviour*/ /* Fifth order transfer function and closed loop transfer function (with additional open loop gain k): */ f:(4*s^2+9*s+6)/(2*s^5+6*s^4+9*s^3+10*s^2+5*s+1); fw:closed_loop(k*f); /* Stability limit with respect to the open loop gain k: */ lim:stability_limit(fw,k); /* Calculation of the stability using the Hurwitz criterion: */ ratsimp(hurwitz(denom(fw))); /* List of three transfer functions: below, at and above the stability limit: */ kli:ev([1.1*k,k,0.9*k],lim,eval,numer); foli:create_list(k*f,k,kli); fwli:closed_loop(foli)$ /* Checking the stability using the predicate function "stablep": */ stablep(fwli); step_response(fwli,xrange=[0,50])$ foli; /* Calculation of the gain crossover frequencies: */ float(gain_crossover(foli)); /* Calculation of the phase margin (unstable closed loop systems have a negative value): */ phase_margin(foli); /* Calculation of the phase crossover frequencies: */ phase_crossover(foli); /* Calculation of the gain margin (unstable closed loop systems have a value less than one): */ gain_margin(foli); /* Calculation of the (absolute) damping (unstable systems have a negative value): */ damping(fwli); /* Calculation of the (relative) damping (i.e. damping ratio), unstable systems have a negative value: */ damping_ratio(fwli); /* List of transfer functions with two symbolic coefficients a and b and another varying coefficient: */ fli:makelist(1/(s^5+3*s^4+k*s^3+a*s^2+b*s+1),k,2,6); /* Graphical representation of the stability limit with respect to the two coefficients a and b: */ stable_area_outlined(fli,a,b,xrange=[0,20],yrange=[0,10]); stable_area([fli[1],fli[2]],a,b,xrange=[2,10],yrange=[0,3]); unstable_area([fli[1],fli[2]],a,b,xrange=[2,10],yrange=[0,3]); /* Controller Design*/ /* Transfer function of a plant to be controlled: */ fs:2/((1+5*s)*(1+s)**2*(1+0.3*s)); /* List with three controllers: I-, PI- and PID-controller: */ [fri,frpi,frpid]:[1/(s*Ti), kr*(1+1/(s*Tn)),(1+s*Ta)*(1+s*Tb)/(s*Tc)]; /* Calculation of the I-controller according to the gain optimum: */ g1:gain_optimum(fs,fri); /* Calculation of the PI-controller according to the gain optimum: */ g2:gain_optimum(fs,frpi); /* Calculation of the PID-controller according to the gain optimum: */ g3:gain_optimum(fs,frpid); /* List with all three dimensioned controllers: */ reli:float(ev([fri,frpi,frpid],[g1,g2,g3])); /* Comparison of the closed loop step responses of the controlled systems: */ step_response(closed_loop(reli*fs)); /* Transfer function of a plant with symbolic coefficients: */ fs:2/((1+a*s)*(1+s**2)*(1+b*s)); /* Gain optimum also handles plants with symbolic coefficients: */ gain_optimum(fs,frpi); /* Optimization*/ /* Transfer function with two symbolic parameters a and b: */ f:1/(s**3+a*s**2+b*s+1); /* Deviation from the stationary value in response to an input step: */ xdev:ratsimp((1-f)/s); /* Integral of squared error (ISE-criterion): */ iise:ise(xdev); /* Optimization of the two symbolic parameters a and b: */ abl:ratsimp(jacobian([iise],[a,b])); realonly:true; res:solve(abl[1],[a,b]); /* Optimum transfer function: */ fopt:ev(f,res); /* State Space Methods*/ /* System matrix A of an RLC-circuit: */ A:matrix([-R/L,0,-1/L,0],[0,-R/L,1/L,-1/L], [1/C1,-1/C1,0,0],[0,1/C1,0,0]); /* Input matrix B of an RLC-circuit: */ B:matrix([1/L],[0],[0],[0]); /* Output matrix C of an RLC-circuit: */ C:matrix([0,0,0,1]); /* A list of the state matrices A, B, and C forms a "system", which can be checked with the predicate functions "systemp" and "nsystemp": */ circuit:[A,B,C]; systemp(circuit); nsystemp(circuit); /* Calculation of the transfer function from the state matrices A, B and C, which can be performed in several ways: */ f:transfer_function(circuit); transfer_function(A,B,C); /* Calculation of the transfer function of the circuit using "impedance_chain" gives the same result, expectedly: */ f:impedance_chain(R+s*L,1/(s*C1),2); /* Transformation of the state matrices into the controller canonical form: */ circ1:controller_canonical_form(f); /* Transformation of the state matrices into the observer canonical form: */ circ2:observer_canonical_form(f); /* Calculation of the controllability matrix: */ h1:ratsimp(controllability_matrix(A,B)); /* Calculation of the observability matrix: */ h2:observability_matrix(circuit); /* The rank of those matrices show, whether the system is controllable and observable: */ [rank(h1),rank(h2)]; kill(A,B,C); /* Time delayed Systems*/ /* (added 2019-01-05) */ /* Open loop system (IT1) with additional dead time: */ Fo:1/(s*(s+3))*exp(-2*s); /* Frequency Response */ absval(Fo); phase(Fo),numer; bode_plot(Fo); nyquist_plot(Fo,xrange=[-1.5,0.5],yrange=[-0.7,0.5]); nichols_plot(Fo,xrange=[-500,0],yrange=[0.1,10]); /* Stability Performance */ gain_crossover(Fo); phase_margin(Fo); phase_crossover(Fo); gain_margin(Fo); /* Closed-Loop System */ /* The closed-loop transfer function can be represented exactly as an (infinite) sum - one addend for each time step. For t