diff --git a/content/research/2020-h-infinity-lecture/_index.md b/content/research/2020-h-infinity-lecture/_index.md index fd3715f..4f5d721 100644 --- a/content/research/2020-h-infinity-lecture/_index.md +++ b/content/research/2020-h-infinity-lecture/_index.md @@ -191,16 +191,16 @@ Having obtained \\(G(s)\\) and \\(G\_d(s)\\), we can transform the system shown Let's define the system parameters on Matlab. ```matlab { linenos=true, linenostart=1 } - k = 1e6; % Stiffness [N/m] - c = 4e2; % Damping [N/(m/s)] - m = 10; % Mass [kg] +k = 1e6; % Stiffness [N/m] +c = 4e2; % Damping [N/(m/s)] +m = 10; % Mass [kg] ``` And now the system dynamics \\(G(s)\\) and \\(G\_d(s)\\). ```matlab { linenos=true, linenostart=4 } - G = 1/(m*s^2 + c*s + k); % Plant - Gd = (c*s + k)/(m*s^2 + c*s + k); % Disturbance +G = 1/(m*s^2 + c*s + k); % Plant +Gd = (c*s + k)/(m*s^2 + c*s + k); % Disturbance ``` The Bode plots of \\(G(s)\\) and \\(G\_d(s)\\) are shown in Figures [Figure 4](#figure--fig:bode-plot-example-afm) and [Figure 5](#figure--fig:bode-plot-example-Gd). @@ -296,7 +296,7 @@ The specifications are: Using `SISOTOOL`, design a controller that fulfills the specifications. ```matlab - sisotool(G) +sisotool(G) ```
Hint @@ -318,16 +318,16 @@ You can follow this procedure: Let's say we came up with the following controller. ```matlab - K = 14e8 * ... % Gain - 1/(s^2) * ... % Double Integrator - 1/(1 + s/2/pi/40) * ... % Low Pass Filter - (1 + s/(2*pi*10/sqrt(8)))/(1 + s/(2*pi*10*sqrt(8))); % Lead +K = 14e8 * ... % Gain + 1/(s^2) * ... % Double Integrator + 1/(1 + s/2/pi/40) * ... % Low Pass Filter + (1 + s/(2*pi*10/sqrt(8)))/(1 + s/(2*pi*10*sqrt(8))); % Lead ``` The bode plot of the Loop Gain is shown in Figure [Figure 8](#figure--fig:loop-gain-manual-afm) and we can verify that we have the wanted stability margins using the `margin` command: ```matlab - [Gm, Pm, ~, Wc] = margin(G*K) +[Gm, Pm, ~, Wc] = margin(G*K); ``` | Requirements | Manual Method | @@ -350,7 +350,7 @@ The synthesis of controllers based on the Loop Shaping method can be automated u Using Matlab, it can be easily performed using the `loopsyn` command: ```matlab - K = loopsyn(G, Lw); +K = loopsyn(G, Lw); ``` where: @@ -399,15 +399,15 @@ Then the gain is adjusted such that the \\(|L\_w(j2 \pi 10)| = 1\\). Using Matlab, we have: ```matlab - Lw = 2.3e3 * ... - 1/(s^2) * ... % Double Integrator - (1 + s/(2*pi*10/sqrt(3)))/(1 + s/(2*pi*10*sqrt(3))); % Lead +Lw = 2.3e3 * ... + 1/(s^2) * ... % Double Integrator + (1 + s/(2*pi*10/sqrt(3)))/(1 + s/(2*pi*10*sqrt(3))); % Lead ``` The \\(\mathcal{H}\_\infty\\) open loop shaping synthesis is then performed using the `loopsyn` command: ```matlab - [K, ~, GAM] = loopsyn(G, Lw); +[K, ~, GAM] = loopsyn(G, Lw); ``` The obtained Loop Gain is shown in Figure [Figure 9](#figure--fig:open-loop-shaping-hinf-L) and matches the specified one by a factor \\(\gamma \approx 2\\). @@ -488,7 +488,7 @@ For a single-input single-output system \\(G(s)\\), it is simply the peak value Let's compute the \\(\mathcal{H}\_\infty\\) norm of our test plant \\(G(s)\\) using the `hinfnorm` function: ```matlab - hinfnorm(G) +hinfnorm(G) ``` ```text @@ -604,7 +604,7 @@ The obtained controller \\(K\\) and the generalized plant are connected as shown Using Matlab, the \\(\mathcal{H}\_\infty\\) Synthesis applied on a Generalized plant can be applied using the `hinfsyn` command ([documentation](https://www.mathworks.com/help/robust/ref/hinfsyn.html)): ```matlab - K = hinfsyn(P, nmeas, ncont); +K = hinfsyn(P, nmeas, ncont); ``` where: @@ -663,11 +663,11 @@ The obtained generalized plant shown in Figure [Figure 15](#figure--fig:mixed-se Using Matlab, the generalized plant can be defined as follows: ```matlab - P = [1 -G; - 0 1; - 1 -G] - P.InputName = {'w', 'u'}; - P.OutputName = {'e', 'u', 'v'}; +P = [1 -G; + 0 1; + 1 -G] +P.InputName = {'w', 'u'}; +P.OutputName = {'e', 'u', 'v'}; ```
@@ -868,17 +868,17 @@ Will follow a discussion about the module margin, a robustness indicator that ca Consider the following plant \\(G\_t(s)\\): ```matlab - w0 = 2*pi*100; - xi = 0.1; - k = 1e7; +w0 = 2*pi*100; +xi = 0.1; +k = 1e7; - Gt = 1/k*(s/w0/4 + 1)/(s^2/w0^2 + 2*xi*s/w0 + 1); +Gt = 1/k*(s/w0/4 + 1)/(s^2/w0^2 + 2*xi*s/w0 + 1); ``` Let's say we have designed a controller \\(K\_t(s)\\) that gives the loop gain shown in Figure [Figure 20](#figure--fig:phase-gain-margin-model-plant). ```matlab - Kt = 1.2e6*(s + w0)/s; +Kt = 1.2e6*(s + w0)/s; ``` The following characteristics can be determined from the Loop gain in Figure [Figure 20](#figure--fig:phase-gain-margin-model-plant): @@ -897,7 +897,7 @@ Or does it? Let's find out. Now let's suppose the controller is implemented in practice, and the "real" plant \\(G\_r(s)\\) as a slightly lower damping factor than the one estimated for the model: ```matlab - xi = 0.03; +xi = 0.03; ``` The obtained "real" loop gain is shown in Figure [Figure 21](#figure--fig:phase-gain-margin-real-plant). @@ -906,7 +906,7 @@ At a frequency little bit above 100Hz, the phase of the loop gain reaches -180 d It is confirmed by checking the stability of the closed loop system: ```matlab - isstable(feedback(Gr,K)) +isstable(feedback(Gr,Kt)) ``` ```text @@ -1116,16 +1116,16 @@ The weighted generalized plant can be defined in Matlab using two techniques: The two solutions below can be used. ```matlab - Pw = [Ws -Ws*G; - 1 -G]; +Pw = [Ws -Ws*G; + 1 -G]; ``` The second solution is however more general, and can also be used when weights are added at the inputs by post-multiplying instead of pre-multiplying. ```matlab - P = [1 -G; - 1 -G]; - Pw = blkdiag(Ws, 1)*P; +P = [1 -G; + 1 -G]; +Pw = blkdiag(Ws, 1)*P; ``` @@ -1157,7 +1157,7 @@ Good guidelines for design of weighting function are given in ( @@ -1210,34 +1210,34 @@ The parameters permit to specify: A Matlab function implementing Equation is shown below: ```matlab - function [W] = generateWeight(args) - arguments - args.G0 (1,1) double {mustBeNumeric, mustBePositive} = 0.1 - args.G1 (1,1) double {mustBeNumeric, mustBePositive} = 10 - args.Gc (1,1) double {mustBeNumeric, mustBePositive} = 1 - args.wc (1,1) double {mustBeNumeric, mustBePositive} = 2*pi - args.n (1,1) double {mustBeInteger, mustBePositive} = 1 - end +function [W] = generateWeight(args) + arguments + args.G0 (1,1) double {mustBeNumeric, mustBePositive} = 0.1 + args.G1 (1,1) double {mustBeNumeric, mustBePositive} = 10 + args.Gc (1,1) double {mustBeNumeric, mustBePositive} = 1 + args.wc (1,1) double {mustBeNumeric, mustBePositive} = 2*pi + args.n (1,1) double {mustBeInteger, mustBePositive} = 1 + end - if (args.Gc <= args.G0 && args.Gc <= args.G1) || (args.Gc >= args.G0 && args.Gc >= args.G1) - eid = 'value:range'; - msg = 'Gc must be between G0 and G1'; - throwAsCaller(MException(eid,msg)) - end + if (args.Gc <= args.G0 && args.Gc <= args.G1) || (args.Gc >= args.G0 && args.Gc >= args.G1) + eid = 'value:range'; + msg = 'Gc must be between G0 and G1'; + throwAsCaller(MException(eid,msg)) + end - s = zpk('s'); + s = zpk('s'); - W = (((1/args.wc) * sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.G1)^(2/args.n)))*s + (args.G0/args.Gc)^(1/args.n)) / ((1/args.G1)^(1/args.n) * (1/args.wc) * sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.G1)^(2/args.n)))*s + (1/args.Gc)^(1/args.n)))^args.n; + W = (((1/args.wc) * sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.G1)^(2/args.n)))*s + (args.G0/args.Gc)^(1/args.n)) / ((1/args.G1)^(1/args.n) * (1/args.wc) * sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.G1)^(2/args.n)))*s + (1/args.Gc)^(1/args.n)))^args.n; - end +end ``` Let's use this function to generate three weights with the same high and low frequency gains, but but different slopes. ```matlab - W1 = generateWeight('G0', 1e2, 'G1', 1/2, 'Gc', 1, 'wc', 2*pi*10, 'n', 1); - W2 = generateWeight('G0', 1e2, 'G1', 1/2, 'Gc', 1, 'wc', 2*pi*10, 'n', 2); - W3 = generateWeight('G0', 1e2, 'G1', 1/2, 'Gc', 1, 'wc', 2*pi*10, 'n', 3); +W1 = generateWeight('G0', 1e2, 'G1', 1/2, 'Gc', 1, 'wc', 2*pi*10, 'n', 1); +W2 = generateWeight('G0', 1e2, 'G1', 1/2, 'Gc', 1, 'wc', 2*pi*10, 'n', 2); +W3 = generateWeight('G0', 1e2, 'G1', 1/2, 'Gc', 1, 'wc', 2*pi*10, 'n', 3); ``` The obtained shapes are shown in Figure [Figure 26](#figure--fig:high-order-weight). @@ -1295,16 +1295,16 @@ We want to design the weighting function \\(W\_s(s)\\) such that: Using Matlab, such weighting function can be generated using the `makeweight` function as shown below: ```matlab - Ws = makeweight(1e3, [2*pi*10, sqrt(2)], 1/2); +Ws = makeweight(1e3, [2*pi*10, sqrt(2)], 1/2); ``` Or using the `generateWeight` function: ```matlab - Ws = generateWeight('G0', 1e3, ... - 'G1', 1/2, ... - 'Gc', sqrt(2), 'wc', 2*pi*10, ... - 'n', 2); +Ws = generateWeight('G0', 1e3, ... + 'G1', 1/2, ... + 'Gc', sqrt(2), 'wc', 2*pi*10, ... + 'n', 2); ``` @@ -1314,24 +1314,24 @@ Or using the `generateWeight` function: Let's say we came up with the following weighting function: ```matlab - Ws = generateWeight('G0', 1e3, ... - 'G1', 1/2, ... - 'Gc', sqrt(2), 'wc', 2*pi*10, ... - 'n', 2); +Ws = generateWeight('G0', 1e3, ... + 'G1', 1/2, ... + 'Gc', sqrt(2), 'wc', 2*pi*10, ... + 'n', 2); ``` The weighting function is then added to the generalized plant. ```matlab - P = [1 -G; - 1 -G]; - Pw = blkdiag(Ws, 1)*P; +P = [1 -G; + 1 -G]; +Pw = blkdiag(Ws, 1)*P; ``` And the \\(\mathcal{H}\_\infty\\) synthesis is performed on the _weighted_ generalized plant. ```matlab - K = hinfsyn(Pw, 1, 1, 'Display', 'on'); +K = hinfsyn(Pw, 1, 1, 'Display', 'on'); ``` ```text @@ -1419,10 +1419,10 @@ Weighting functions: ```matlab - P = [1 -G - 0 1 - 1 -G]; - Pw = blkdiag(W1, W2, 1)*P; +P = [1 -G + 0 1 + 1 -G]; +Pw = blkdiag(W1, W2, 1)*P; ``` @@ -1441,10 +1441,10 @@ Weighting functions: ```matlab - P = [1 -G - 0 G - 1 -G]; - Pw = blkdiag(W1, W2, 1)*P; +P = [1 -G + 0 G + 1 -G]; +Pw = blkdiag(W1, W2, 1)*P; ``` @@ -1463,10 +1463,10 @@ Weighting functions: ```matlab - P = [1 -1 - G -G - G -G]; - Pw = blkdiag(W1, W2, 1)*P; +P = [1 -1 + G -G + G -G]; +Pw = blkdiag(W1, W2, 1)*P; ``` @@ -1486,11 +1486,11 @@ Weighting functions: ```matlab - P = [1 -G - 0 1 - 0 G - 1 -G]; - Pw = blkdiag(W1, W2, W3, 1)*P; +P = [1 -G + 0 1 + 0 G + 1 -G]; +Pw = blkdiag(W1, W2, W3, 1)*P; ``` @@ -1510,11 +1510,11 @@ Weighting functions: ```matlab - P = [1 -1 - G -G - 0 1 - G -G]; - Pw = blkdiag(W1, W2, W3, 1)*P; +P = [1 -1 + G -G + 0 1 + G -G]; +Pw = blkdiag(W1, W2, W3, 1)*P; ``` @@ -1535,10 +1535,10 @@ Weighting functions: ```matlab - P = [ 1 -G -G - 0 0 1 - 1 -G -G]; - Pw = blkdiag(W1, W2, 1)*P*blkdiag(1, W3, 1); +P = [ 1 -G -G + 0 0 1 + 1 -G -G]; +Pw = blkdiag(W1, W2, 1)*P*blkdiag(1, W3, 1); ``` @@ -1692,37 +1692,37 @@ Let's first convert the system of Figure [Figure 36](#figure--fig:ex-test-system The two transfer functions present in the system are derived and defined below: ```matlab { linenos=true, linenostart=1 } - k = 1e6; % Stiffness [N/m] - c = 4e2; % Damping [N/(m/s)] - m = 10; % Mass [kg] +k = 1e6; % Stiffness [N/m] +c = 4e2; % Damping [N/(m/s)] +m = 10; % Mass [kg] - % Control Plant - G = 1/(m*s^2 + c*s + k); - % Disturbance dynamics - Gd = (c*s + k)/(m*s^2 + c*s + k); +% Control Plant +G = 1/(m*s^2 + c*s + k); +% Disturbance dynamics +Gd = (c*s + k)/(m*s^2 + c*s + k); ``` We also define the inputs signals that will be used for time domain simulations. They are graphically shown in Figure [Figure 38](#figure--fig:ex-inputs-signals). ```matlab { linenos=true, linenostart=9 } - % Time Vector - t = 0:1e-4:0.5; +% Time Vector +t = 0:1e-4:0.5; - % Reference Input - r = zeros(size(t)); - r(t>0.05 & t<=0.25) = 1e-4*(t(t>0.05 & t<=0.25)-0.05); - r(t>0.25) = 2e-5; +% Reference Input +r = zeros(size(t)); +r(t>0.05 & t<=0.25) = 1e-4*(t(t>0.05 & t<=0.25)-0.05); +r(t>0.25) = 2e-5; - % Measurement Noise - Fs = 1e3; % Sampling Frequency [Hz] - Ts = 1/Fs; % Sampling Time [s] - n = sqrt(Fs/2)*randn(1, length(t)); % Signal with an ASD equal to one - n = lsim(1e-6*(s + 2*pi*1e2)^2/(s + 2*pi*1e3)^2/(1+s/2/pi/500), n, t)'; % Shaped noise +% Measurement Noise +Fs = 1e3; % Sampling Frequency [Hz] +Ts = 1/Fs; % Sampling Time [s] +n = sqrt(Fs/2)*randn(1, length(t)); % Signal with an ASD equal to one +n = lsim(1e-6*(s + 2*pi*1e2)^2/(s + 2*pi*1e3)^2/(1+s/2/pi/500), n, t)'; % Shaped noise - % Disturbance - d = zeros(size(t)); - d(t>0.3) = 5e-6; +% Disturbance +d = zeros(size(t)); +d(t>0.3) = 5e-6; ``` @@ -1738,25 +1738,25 @@ We also define the generalized plant corresponding to the system and that will b The Generalized plant of Figure [Figure 39](#figure--fig:ex-general-plant-sim) is defined on Matlab as follows: ```matlab { linenos=true, linenostart=26 } - Psim = [0 0 Gd G - 0 0 0 1 - 1 -1 -Gd -G]; +Psim = [0 0 Gd G + 0 0 0 1 + 1 -1 -Gd -G]; - Psim.InputName = {'r', 'n', 'd', 'u'}; - Psim.OutputName = {'y', 'u', 'e'}; +Psim.InputName = {'r', 'n', 'd', 'u'}; +Psim.OutputName = {'y', 'u', 'e'}; ``` Time domain simulations will be performed by first computing the closed-loop system using the `lft` command and then using the `lsim` command on the closed-loop system: ```matlab { linenos=true, linenostart=32 } - % Compute the closed-Loop System, K is the controller - P_CL = lft(Psim, K); +% Compute the closed-Loop System, K is the controller +P_CL = lft(Psim, K); - % Time simulation of the closed-loop system with specified inputs - z = lsim(P_CL, [r; n; d], t); - % The two outputs are - y = z(:,1); % Output Motion [m] - u = z(:,2); % Input usage [N] +% Time simulation of the closed-loop system with specified inputs +z = lsim(P_CL, [r; n; d], t); +% The two outputs are +y = z(:,1); % Output Motion [m] +u = z(:,2); % Input usage [N] ``` @@ -1804,10 +1804,10 @@ To do so, we use to generalized plant shown in Figure [Figure 40](#figure--fig:e This generalized plant is defined on Matlab as follows: ```matlab { linenos=true, linenostart=40 } - P = [1 -1 - G -G - 0 1 - G -G]; +P = [1 -1 + G -G + 0 1 + G -G]; ``` @@ -1818,7 +1818,7 @@ However, to performed the \\(\mathcal{H}\_\infty\\) loop shaping, we have to inc We obtain the **weighted generalized plant** in Figure [Figure 41](#figure--fig:ex-general-weighted-plant), and that is computed using Matlab as follows: ```matlab { linenos=true, linenostart=44 } - Pw = blkdiag(W1, W2, W3, 1)*P; +Pw = blkdiag(W1, W2, W3, 1)*P; ``` @@ -1828,7 +1828,7 @@ We obtain the **weighted generalized plant** in Figure [Figure 41](#figure--fig: Finlay, performing the \\(\mathcal{H}\_\infty\\) Shaping of \\(S\\), \\(GS\\) and \\(T\\) is as simple as ruining the `hinfsyn` command: ```matlab { linenos=true, linenostart=45 } - K = hinfsyn(Pw, 1, 1); +K = hinfsyn(Pw, 1, 1); ``` Now let's shape the three closed-loop transfer functions sequentially: diff --git a/themes/clean-even b/themes/clean-even index 365fae8..9bbe44d 160000 --- a/themes/clean-even +++ b/themes/clean-even @@ -1 +1 @@ -Subproject commit 365fae817ca92db8b0c55def67dc2bea7e932a70 +Subproject commit 9bbe44d6d2a988c4f733d531cbe05fca1418727b