Control engineer: MQL recipes

Recipes

Gain and phase margins of a type 1 loop
PI control of a first-order plant: exact and simulated step response
State feedback by pole placement (Ackermann)
Nyquist plot of an open loop
Lyapunov stability by solving At P + P A = -Q
Time-domain specs from the closed-loop poles
Continuous poles of an identified discrete model
State transition matrix by modal decomposition
Relay autotuning of a temperature loop
Relative gain array of a distillation column

executed on MentDB

Gain and phase margins of a type 1 loop

For the open loop L(s) = 8 / (s (s+1) (s+5)), find the gain crossover frequency, evaluate L(jw) with complex arithmetic to get the phase margin and the delay margin (also the margin left with a 0.2 s actuator delay e^(-j w tau)), compute the gain margin at the phase crossover w = sqrt(5), and the closed-loop poles with the damping of the dominant pair.

Pipeline

01 · crossover|L(jw)| = 1num solve→02 · L(jw)Product and quotient of complex numberscomplex mul · complex add · complex div→03 · marginsArgument and modulus, effect of a delaycomplex arg · complex abs · complex to_polar · complex exp→04 · polesRoots of s^3 + 6 s^2 + 5 s + 8num poly_roots · complex re

Script

admin

#1. Open loop L of s = K / s s+1 s+5 with K = 8;
-> "[K]" 8;
#2. Gain crossover: modulus of L at j w equal to 1;
-> "[wc]" (num solve (concat [K] "/(x*sqrt(x^2+1)*sqrt(x^2+25))-1") "x" 0.1 10);
-> "[jw]" (concat [wc] "i");
-> "[L]" (complex div [K] (complex mul [jw] (complex mul (complex add [jw] 1) (complex add [jw] 5))));
-> "[pm]" (math calc (concat "180+" (math rad_to_deg (complex arg [L]))));
#3. Delay margin PM / wc, and phase margin with a 0.2 s actuator delay exp of -j w tau;
-> "[dm]" (math calc (concat (math deg_to_rad [pm]) "/[wc]"));
-> "[pmd]" (math calc (concat "180+" (math rad_to_deg (complex arg (complex mul [L] (complex exp (concat "-" (math calc "0.2*[wc]") "i")))))));
#4. Phase crossover at w = sqrt 5 and gain margin;
-> "[w180]" (math sqrt 5);
-> "[L180]" (complex div [K] (complex mul (concat [w180] "i") (complex mul (complex add (concat [w180] "i") 1) (complex add (concat [w180] "i") 5))));
-> "[gm]" (/ 1 (complex abs [L180]));
#5. Closed-loop poles: roots of s^3 + 6 s^2 + 5 s + K, damping of the dominant pair;
-> "[roots]" (num poly_roots (concat "[" [K] ", 5, 6, 1]"));
-> "[p]" (json load select "rx_cl" [roots] "/[1]");
concat "wc = " (math round [wc] 3) " rad/s, |L(jwc)| = " (math round (json load select "rx_lp" (complex to_polar [L]) "/r") 6) ", phase margin " (math round [pm] 1) " deg, delay margin " (math round [dm] 3) " s, phase margin with 0.2 s delay " (math round [pmd] 1) " deg, gain margin " (math round [gm] 3) " = " (math round (math calc "20*log10([gm])") 1) " dB; closed-loop poles " [roots] ", dominant pair zeta = " (math round (math calc (concat "-" (complex re [p]) "/" (complex abs [p]))) 3) ", wn = " (math round (complex abs [p]) 3) " rad/s"

Result

mentdb

wc = 1.069 rad/s, |L(jwc)| = 1, phase margin 31 deg, delay margin 0.507 s, phase margin with 0.2 s delay 18.8 deg, gain margin 3.75 = 11.5 dB; closed-loop poles ["-5.344537368239695","-0.3277313158801529-1.178748342804254i","-0.32773131588015286+1.178748342804254i"], dominant pair zeta = 0.268, wn = 1.223 rad/s

Go further

Draw the Nyquist curve (see the recipe of this profile)
Find the gain K that gives a 45 deg phase margin with num solve
Check the Routh limit K = 30 with num poly_roots

executed on MentDB

PI control of a first-order plant: exact and simulated step response

A plant 1/(2s+1) under PI control (Kp = 3, Ki = 4) gives the closed loop T(s) = (3s+4)/(2s^2+4s+4). The CAS gives the exact poles, the final value (final value theorem) and the step response by inverse Laplace transform; a numerical simulation of the loop gives the overshoot and the peak time.

Pipeline

01 · polesRoots of 2 s^2 + 4 s + 4cas solve→02 · final valuelim s->0 of s T(s) / scas limit→03 · responseInverse Laplace transform, simplifiedcas eval · cas simplify→04 · simulatePlant state and error integralnum ode→05 · overshootMax of y and first peakstat max · dsp peaks

Script

admin

#1. Plant 1 / 2s+1 with a PI controller Kp = 3, Ki = 4: closed loop T = 3s+4 / 2s^2+4s+4;
-> "[T]" "(3*s+4)/(2*s^2+4*s+4)";
#2. Exact poles, final value and step response with the CAS;
-> "[poles]" (cas solve "2*s^2+4*s+4 = 0" "s");
-> "[final]" (cas limit (concat "s*" [T] "/s") "s" 0);
-> "[y]" (cas eval (concat "InverseLaplaceTransform(" [T] "/s, s, t)"));
-> "[y]" (cas simplify [y]);
#3. Time simulation: plant state y and integral of the error z;
-> "[sim]" (num ode "{\"y\": \"(-y+3*(1-y)+4*z)/2\", \"z\": \"1-y\"}" "{\"y\": 0, \"z\": 0}" 0 8 "{\"points\": 161}");
-> "[os]" (math calc (concat "100*(" (stat max (json load select "rx_sim" [sim] "/y")) "-1)"));
-> "[tp]" (math calc (concat (json load select "rx_pk" (dsp peaks (json select "rx_sim" "/y")) "/indices[0]") "*0.05"));
concat "poles " [poles] ", final value " [final] ", y(t) = " [y] ", overshoot " (math round [os] 1) " % at t = " (math round [tp] 2) " s"

Result

mentdb

poles [{"s":"-1-I"},{"s":"-1+I"}], final value 1, y(t) = 1+(-2*cos(t)+sin(t))/(2*e^t), overshoot 11.9 % at t = 1.9 s

Go further

Draw the simulation with plot ode
Add the derivative term and tune Kd to remove the overshoot
Check the anti-windup by saturating the command in the ODE with min and max

executed on MentDB

State feedback by pole placement (Ackermann)

A DC motor model (position, speed) must get closed-loop poles at -20 +/- 20j. The script checks the controllability, expands the desired characteristic polynomial with the CAS, applies Ackermann's formula K = [0 1] C^-1 phi(A), and verifies the eigenvalues, the trace and the determinant of A - BK.

Pipeline

01 · controllabilityC = [B AB] and its rankmatrix mul · matrix get · matrix rank→02 · polynomial(s+20-20j)(s+20+20j)cas expand→03 · phi(A)A^2 + 40 A + 800 Imatrix pow · matrix scalar_mul · matrix add · matrix identity→04 · gainLast row of C^-1 times phi(A)matrix inverse · matrix row→05 · checkEigenvalues of A - BKmatrix sub · matrix eigen · matrix trace · matrix det

Script

admin

#1. DC motor state x = position, speed: A and B;
-> "[A]" "[[0, 1], [0, -10]]"; -> "[B]" "[[0], [100]]";
#2. Controllability matrix B, AB and its rank;
-> "[AB]" (matrix mul [A] [B]);
-> "[C]" (concat "[[0, " (matrix get [AB] 1 1) "], [100, " (matrix get [AB] 2 1) "]]");
-> "[rk]" (matrix rank [C]);
#3. Ackermann: desired poles -20 +/- 20j, characteristic polynomial, phi of A = A^2 + 40 A + 800 I;
-> "[cp]" (cas expand "(s+20-20*I)*(s+20+20*I)");
-> "[phi]" (matrix add (matrix add (matrix pow [A] 2) (matrix scalar_mul [A] 40)) (matrix scalar_mul (matrix identity 2) 800));
-> "[K]" (matrix mul (concat "[" (matrix row (matrix inverse [C]) 2) "]") [phi]);
#4. Check: eigenvalues, trace and determinant of A - B K;
-> "[Acl]" (matrix sub [A] (matrix mul [B] [K]));
concat "rank " [rk] ", desired polynomial " [cp] ", K = " [K] ", A - BK = " [Acl] ", eigenvalues " (json load select "rx_eg" (matrix eigen [Acl]) "/values[0]") " and " (json select "rx_eg" "/values[1]") ", trace " (matrix trace [Acl]) ", det " (matrix det [Acl])

Result

mentdb

rank 2, desired polynomial 800+40*s+s^2, K = [[8,0.3]], A - BK = [[0,1],[-800,-40]], eigenvalues -20+20i and -20-20i, trace -40, det 800

Go further

Simulate the closed loop with num ode and draw it with plot ode
Design an observer with the same formula on the transposed system (matrix transpose)
Compute the LQR gain for comparison

executed on MentDB

Nyquist plot of an open loop

Draw the Nyquist curve of L(jw) = 8 / (jw (jw+1) (jw+5)) for w from 0.3 to 20 rad/s on a logarithmic grid, computing each point with complex arithmetic: the curve crosses the real axis at -0.267 (gain margin 3.75), on the right of -1, so the closed loop is stable.

Pipeline

01 · gridw = 0.3 x 10^(k/33)math calc→02 · L(jw)Complex product and quotientcomplex mul · complex add · complex div→03 · partsReal and imaginary partscomplex re · complex im→04 · plotIm against Replot xy

Script

admin

#1. Open loop L of s = 8 / s s+1 s+5: Nyquist curve for w from 0.3 to 20 rad/s;
-> "[re]" "["; -> "[im]" "["; -> "[sep]" "";
for (-> "[k]" 0) (<= [k] 60) (++ "[k]") {
  -> "[jw]" (concat (math calc "0.3*10^([k]/33)") "i");
  -> "[L]" (complex div 8 (complex mul [jw] (complex mul (complex add [jw] 1) (complex add [jw] 5))));
  -> "[re]" (concat [re] [sep] (complex re [L])); -> "[im]" (concat [im] [sep] (complex im [L])); -> "[sep]" ",";
};
#2. The curve crosses the real axis at -1/GM = -0.2667: it does not encircle -1;
plot xy (concat [re] "]") (concat [im] "]")

Result

mentdb

The chart is shown in the Plot tab of the editor.

Go further

Add the mirror image for negative frequencies
Draw the Bode magnitude with plot function and 20*log10
Find the real axis crossing exactly with num solve on the imaginary part

executed on MentDB

Lyapunov stability by solving At P + P A = -Q

For the linear system x' = A x with A = [[0, 1], [-2, -3]], solve the Lyapunov equation At P + P A = -I for the symmetric matrix P, check the residual, and prove that P is positive definite with its Cholesky factor, so that V = xt P x is a Lyapunov function.

Pipeline

01 · solveEquations of the entries of Pnum solve_system→02 · residualAt P + P A + Imatrix transpose→03 · normResidual normmatrix norm_frobenius→04 · definiteP = L Lt existsmatrix cholesky

Script

admin

#1. Linear plant x' = A x and weighting Q = I;
-> "[A]" "[[0, 1], [-2, -3]]";
#2. Lyapunov equation At P + P A = -Q, P symmetric with unknowns p1 p2 p3;
json load "rx_lyap" (num solve_system "[\"-4*p2 = -1\", \"p1-3*p2-2*p3 = 0\", \"2*p2-6*p3 = -1\"]" "{\"p1\": 1, \"p2\": 0, \"p3\": 1}");
-> "[p1]" (json select "rx_lyap" "/solution/p1");
-> "[p2]" (json select "rx_lyap" "/solution/p2");
-> "[p3]" (json select "rx_lyap" "/solution/p3");
-> "[P]" (concat "[[" [p1] ", " [p2] "], [" [p2] ", " [p3] "]]");
#3. Check the residual At P + P A + I with matrix products;
-> "[R]" (matrix add (matrix add (matrix mul (matrix transpose [A]) [P]) (matrix mul [P] [A])) (matrix identity 2));
#4. P positive definite if its Cholesky factor exists: V = xt P x is a Lyapunov function;
concat "P = " [P] " | residual norm = " (matrix norm_frobenius [R]) " | Cholesky L = " (matrix cholesky [P]) " | det P = " (matrix det [P]) " > 0: the origin is asymptotically stable"

Result

mentdb

P = [[1.25, 0.25], [0.25, 0.25]] | residual norm = 0 | Cholesky L = {"L": [[1.118033988749895,0],[0.22360679774997896,0.4472135954999579]]} | det P = 0.25 > 0: the origin is asymptotically stable

Go further

Check the result with matrix eigen on A: both eigenvalues have a negative real part.
Use the same equations for the observability Gramian of a state model.

executed on MentDB

Time-domain specs from the closed-loop poles

The closed loop has the characteristic polynomial s^2 + 4s + 16.25. The script finds the dominant poles, reads the natural frequency and damping from their position in the s plane and predicts overshoot, peak time and 2 % settling time.

Pipeline

01 · polesRoots of the characteristic polynomialnum poly_roots→02 · positionsigma and omega_dcomplex im→03 · geometrywn = hypot, zeta = cos(atan2)math atan2→04 · specsMp, tp, tsmath exp

Script

admin

#1. Closed-loop characteristic polynomial s^2 + 4 s + 16.25;
json load "rx_p4f_poles" (num poly_roots "[16.25, 4, 1]");
-> "[p]" (json select "rx_p4f_poles" "/[1]");
json unload "rx_p4f_poles";
-> "[sigma]" (- 0 (complex re [p])); -> "[wd]" (complex im [p]);
#2. Natural frequency = distance to the origin, damping = cosine of the pole angle;
-> "[wn]" (math hypot [sigma] [wd]);
-> "[zeta]" (math cos (math atan2 [wd] [sigma]));
#3. Overshoot, peak time and settling time at 2 pct;
-> "[Mp]" (math exp (/ (* -1 (* (math pi) [zeta])) (math sqrt (- 1 (* [zeta] [zeta])))));
-> "[tp]" (/ (math pi) [wd]);
-> "[ts]" (/ 4 [sigma]);
concat "poles -" [sigma] " +/- j" [wd] " | wn = " (math round [wn] 4) " rad/s, zeta = " (math round [zeta] 4) " | overshoot = " (math round (* 100 [Mp]) 2) " pct, tp = " (math round [tp] 3) " s, ts = " (math round [ts] 2) " s"

Result

mentdb

poles -2 +/- j3.5 | wn = 4.0311 rad/s, zeta = 0.4961 | overshoot = 16.61 pct, tp = 0.898 s, ts = 2 s

Go further

Check the overshoot on the simulated step response with num ode and stat max
Find the gain that gives zeta = 0.7 with num solve

executed on MentDB

Continuous poles of an identified discrete model

An ARX model identified on a test bench at T = 50 ms has the denominator z^2 - 1.6z + 0.68. The script maps its poles back to the s plane with the complex logarithm s = ln(z)/T and reads the natural frequency, damping and settling time of the plant.

Pipeline

01 · rootsPoles in the z planenum poly_roots→02 · maps = ln(z) / Tcomplex log→03 · scaleDivide by the sample timecomplex div→04 · readwn = |s|, zeta = -Re(s)/|s|complex abs

Script

admin

#1. Denominator of the identified model, sample time 50 ms;
-> "[T]" 0.05;
json load "rx_p4f_z" (num poly_roots "[0.68, -1.6, 1]");
-> "[z]" (json select "rx_p4f_z" "/[1]");
json unload "rx_p4f_z";
#2. Back to the s plane: s = ln z / T, principal value;
-> "[s]" (complex div (complex log [z]) [T]);
#3. Natural frequency, damping and settling time of the plant;
-> "[wn]" (complex abs [s]);
-> "[zeta]" (/ (- 0 (complex re [s])) [wn]);
concat "z = " (math round (complex re [z]) 4) " +/- j" (math round (complex im [z]) 4) " (|z| = " (math round (complex abs [z]) 4) ") -> s = " (math round (complex re [s]) 3) " +/- j" (math round (complex im [s]) 3) " | wn = " (math round [wn] 3) " rad/s = " (math round (/ [wn] (* 2 (math pi))) 3) " Hz, zeta = " (math round [zeta] 3) ", ts = " (math round (/ 4 (* [zeta] [wn])) 2) " s"

Result

mentdb

z = 0.8 +/- j0.2 (|z| = 0.8246) -> s = -3.857 +/- j4.9 | wn = 6.235 rad/s = 0.992 Hz, zeta = 0.619, ts = 1.04 s

Go further

Check that the sample rate is 10 to 30 times the bandwidth
Go the other way with complex exp to discretize a continuous design

executed on MentDB

State transition matrix by modal decomposition

For the state model x' = Ax with A = [[0, 1], [-2, -3]], the script diagonalizes A, builds e^(At) = V diag(e^(lambda t)) V^-1 at t = 0.5 s and computes the free response from x0 = [1, 0], checked against the exact solution.

Pipeline

01 · modesEigenvalues and eigenvectorsmatrix eigen→02 · diagonaldiag(e^(lambda t))matrix diag→03 · transitionV D V^-1matrix inverse→04 · responsex(t) = Phi x0matrix mul

Script

admin

#1. Modes of A, eigenvectors given as rows;
-> "[A]" "[[0, 1], [-2, -3]]"; -> "[t]" 0.5;
json load "rx_p4f_eig" (matrix eigen [A]);
-> "[l1]" (json select "rx_p4f_eig" "/values/[0]"); -> "[l2]" (json select "rx_p4f_eig" "/values/[1]");
-> "[V]" (matrix transpose (json select "rx_p4f_eig" "/vectors"));
json unload "rx_p4f_eig";
#2. Phi = V . diag of e^lambda t . V^-1;
-> "[D]" (matrix diag (concat "[" (math exp (* [l1] [t])) ", " (math exp (* [l2] [t])) "]"));
-> "[Phi]" (matrix mul (matrix mul [V] [D]) (matrix inverse [V]));
#3. Free response from x0 = 1, 0 and exact position 2 e^-t - e^-2t;
-> "[x]" (matrix mul [Phi] "[1, 0]");
-> "[exact]" (math calc "2*exp(-[t])-exp(-2*[t])");
json load "rx_p4f_x" [x];
concat "modes " [l1] ", " [l2] " | trace Phi(0.5) = " (math round (matrix trace [Phi]) 6) " = e^-0.5 + e^-1 | x(0.5) = [" (math round (json select "rx_p4f_x" "/[0]") 6) ", " (math round (json select "rx_p4f_x" "/[1]") 6) "], exact position " (math round [exact] 6)

Result

mentdb

modes -1, -2 | trace Phi(0.5) = 0.97441 = e^-0.5 + e^-1 | x(0.5) = [0.845182, -0.477302], exact position 0.845182

Go further

Repeat with matrix pow of the discrete transition matrix to step in time
Check the result with num ode on the same system

executed on MentDB

Relay autotuning of a temperature loop

A relay feedback test (amplitude 10 %) puts a temperature loop into a limit cycle. From the measured oscillation amplitudes and half periods, the script computes the ultimate gain with the describing function of the relay, the ultimate period, and the Ziegler-Nichols and Tyreus-Luyben PID settings.

Pipeline

01 · amplitudeMean amplitude of the limit cyclemath avg→02 · periodPu and its spreadstat pvariance→03 · gainKu = 4d / (pi a)math pi→04 · tuningZN and Tyreus-Luyben PIDmath round

Script

admin

#1. Relay amplitude d = 10 pct, peak to peak amplitudes of the limit cycle in degC;
-> "[d]" 10;
-> "[a]" (/ (math avg "[7.9, 8.1, 8.0, 8.2, 7.8]") 2);
#2. Half periods between zero crossings in s, and their spread;
-> "[half]" "[41.5, 42.3, 41.8, 42.6, 41.9, 42.1]";
-> "[Pu]" (* 2 (math avg [half]));
-> "[cv]" (/ (math sqrt (stat pvariance [half])) (math avg [half]));
#3. Ultimate gain from the describing function of the relay;
-> "[Ku]" (/ (* 4 [d]) (* (math pi) [a]));
#4. Ziegler-Nichols and Tyreus-Luyben PID settings;
concat "Ku = " (math round [Ku] 3) " pct/degC, Pu = " (math round [Pu] 1) " s, spread " (math round (* 100 [cv]) 2) " pct | ZN: Kp = " (math round (* 0.6 [Ku]) 3) ", Ti = " (math round (/ [Pu] 2) 1) " s, Td = " (math round (/ [Pu] 8) 1) " s | TL: Kp = " (math round (/ [Ku] 2.2) 3) ", Ti = " (math round (* 2.2 [Pu]) 1) " s, Td = " (math round (/ [Pu] 6.3) 1) " s"

Result

mentdb

Ku = 3.183 pct/degC, Pu = 84.1 s, spread 0.84 pct | ZN: Kp = 1.91, Ti = 42 s, Td = 10.5 s | TL: Kp = 1.447, Ti = 184.9 s, Td = 13.3 s

Go further

Add a hysteresis to the relay and correct Ku with the arcsine term
Simulate the tuned loop with num ode to compare the overshoot of both tunings

executed on MentDB

Relative gain array of a distillation column

For the steady-state gains of the Wood-Berry column (top and bottom compositions vs reflux and steam), the script computes the determinant from the LU pivots, the relative gain lambda11 and chooses the input-output pairing.

Pipeline

01 · decomposeP.G = L.Umatrix lu→02 · determinantProduct of the pivots, sign of Pjson select→03 · rgalambda11 = g11 g22 / detmatrix get→04 · checkg11 times inverse 11matrix inverse

Script

admin

#1. Steady-state gains, outputs xD and xB, inputs reflux and steam;
-> "[G]" "[[12.8, -18.9], [6.6, -19.4]]";
#2. LU with partial pivoting: det = product of the pivots times det P;
json load "rx_p4f_lu" (matrix lu [G]);
-> "[u11]" (json select "rx_p4f_lu" "/U/[0]/[0]"); -> "[u22]" (json select "rx_p4f_lu" "/U/[1]/[1]");
-> "[sp]" (if (== (json select "rx_p4f_lu" "/P/[0]/[0]") 1) {1} {-1});
json unload "rx_p4f_lu";
-> "[det]" (* [sp] (* [u11] [u22]));
#3. Relative gain lambda11 and cross-check g11 times the 1,1 term of the inverse;
-> "[l11]" (/ (* (matrix get [G] 1 1) (matrix get [G] 2 2)) [det]);
-> "[l11b]" (* (matrix get [G] 1 1) (matrix get (matrix inverse [G]) 1 1));
concat "det G = " (math round [det] 4) " (matrix det " (math round (matrix det [G]) 4) ") | lambda11 = " (math round [l11] 4) " (check " (math round [l11b] 4) "), lambda12 = " (math round (- 1 [l11]) 4) " | " (if (> [l11] 0.5) {"pair reflux-xD and steam-xB"} {"pair reflux-xB and steam-xD"})

Result

mentdb

det G = -123.58 (matrix det -123.58) | lambda11 = 2.0094 (check 2.0094), lambda12 = -1.0094 | pair reflux-xD and steam-xB

Go further

Add the time delays of the model and compute the RGA at the crossover frequency with complex div
Check the Niederlinski index det G / (g11 g22) for the chosen pairing