%% Clear Workspace and Close figures clear; close all; clc; %% Intialize Laplace variable s = zpk('s'); % Let's define the system parameters on Matlab. 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)$. G = 1/(m*s^2 + c*s + k); % Plant Gd = (c*s + k)/(m*s^2 + c*s + k); % Disturbance % The Sensitivity Function % <> % The sensitivity function is indisputably the most important closed-loop transfer function of a feedback system. % In this section, we will see how the shape of the sensitivity function will impact the performances of the closed-loop system. % Suppose we have developed a "/reference/" controller $K_r(s)$ and made three small changes to obtained three controllers $K_1(s)$, $K_2(s)$ and $K_3(s)$. % The obtained sensitivity functions for these four controllers are shown in Figure [[fig:sensitivity_shape_effect]] and the corresponding step responses are shown in Figure [[fig:sensitivity_shape_effect_step]]. % The comparison of the sensitivity functions shapes and their effect on the step response is summarized in Table [[tab:compare_sensitivity_shapes]]. % #+name: tab:compare_sensitivity_shapes % #+caption: Comparison of the sensitivity function shape and the corresponding step response for the three controller variations % #+attr_latex: :environment tabularx :width 0.9\linewidth :align lXX % | *Controller* | *Sensitivity Function Shape* | *Change of the Step Response* | % |--------------+----------------------------------------------------+----------------------------------| % | $K_1(s)$ | Larger bandwidth $\omega_b$ | Faster rise time | % | $K_2(s)$ | Larger peak value $\Vert S\Vert_\infty$ | Large overshoot and oscillations | % | $K_3(s)$ | Larger low frequency gain $\vert S(j\cdot 0)\vert$ | Larger static error | wc = 2*pi*1; L_w = 8; wi = 2*pi*0.02; Kr = 1/((s + wi)^2) * ... % Double Integrator (1 + s/(wc/sqrt(L_w)))/(1 + s/(wc*sqrt(L_w))); % Lead Kr = Kr/abs(evalfr(Kr*G, 1j*wc)); wc = 2*pi*2; L_w = 8; wi = 2*pi*0.045; K1 = 1/((s + wi)^2) * ... % Double Integrator (1 + s/(wc/sqrt(L_w)))/(1 + s/(wc*sqrt(L_w))); % Lead K1 = K1/abs(evalfr(K1*G, 1j*wc)); wc = 2*pi*1; L_w = 2; wi = 2*pi*0.03; K2 = 1/((s + wi)^2) * ... % Double Integrator (1 + s/(wc/sqrt(L_w)))/(1 + s/(wc*sqrt(L_w))); % Lead K2 = K2/abs(evalfr(K2*G, 1j*wc)); wc = 2*pi*1; L_w = 8; wi = 2*pi*0.2; K3 = 1/((s + wi)^2) * ... % Double Integrator (1 + s/(wc/sqrt(L_w)))/(1 + s/(wc*sqrt(L_w))); % Lead K3 = K3/abs(evalfr(K3*G, 1j*wc)); Sr = 1/(1 + Kr*G); S1 = 1/(1 + K1*G); S2 = 1/(1 + K2*G); S3 = 1/(1 + K3*G); Tr = Kr*G/(1 + Kr*G); T1 = K1*G/(1 + K1*G); T2 = K2*G/(1 + K2*G); T3 = K3*G/(1 + K3*G); freqs = logspace(-2, 2, 1000); figure; hold on; plot(freqs, abs(squeeze(freqresp(S1, freqs, 'Hz'))), 'DisplayName', '$K_1(s)$'); plot(freqs, abs(squeeze(freqresp(S2, freqs, 'Hz'))), 'DisplayName', '$K_2(s)$'); plot(freqs, abs(squeeze(freqresp(S3, freqs, 'Hz'))), 'DisplayName', '$K_3(s)$'); plot(freqs, abs(squeeze(freqresp(Sr, freqs, 'Hz'))), 'k-', 'DisplayName', '$K_r(s)$'); hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); xlabel('Frequency [Hz]'); ylabel('Sensitivity Magnitude'); hold off; legend('location', 'southeast'); ylim([1e-3, 3]); % #+name: fig:sensitivity_shape_effect % #+caption: Sensitivity function magnitude $|S(j\omega)|$ corresponding to the reference controller $K_r(s)$ and the three modified controllers $K_i(s)$ % #+RESULTS: % [[file:figs/sensitivity_shape_effect.png]] t = linspace(0, 5, 1000); figure; hold on; plot(t, step(T1, t), 'DisplayName', '$K_1(s)$') plot(t, step(T2, t), 'DisplayName', '$K_2(s)$') plot(t, step(T3, t), 'DisplayName', '$K_3(s)$') plot(t, step(Tr, t), 'k-', 'DisplayName', '$K_r(s)$') hold off xlabel('Time [s]'); ylabel('Step Response'); legend('location', 'northeast'); ylim([0, 1.7]) % Robustness: Module Margin % <> % Let's start this section by an example demonstrating why the phase and gain margins might not be good indicators of robustness. % Will follow a discussion about the module margin, a robustness indicator that can be linked to the $\mathcal{H}_\infty$ norm of $S$ and that will prove to be very useful. % #+begin_exampl % Consider the following plant $G_t(s)$: 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); % Let's say we have designed a controller $K_t(s)$ that gives the loop gain shown in Figure [[fig:phase_gain_margin_model_plant]]. Kt = 1.2e6*(s + w0)/s; % The following characteristics can be determined from the Loop gain in Figure [[fig:phase_gain_margin_model_plant]]: % - Control bandwidth of $\approx 10\, \text{Hz}$ % - Infinite gain margin (the phase of the loop-gain never reaches $-180^o$) % - More than $90^o$ of phase margin % This clearly indicate very good robustness of the closed-loop system! % Or does it? Let's find out. freqs = logspace(0, 3, 1000); figure; tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); ax1 = nexttile([2,1]); hold on; plot(freqs, abs(squeeze(freqresp(Gt*Kt, freqs, 'Hz'))), 'DisplayName', '$L(s)$'); hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); ylabel('Magnitude'); set(gca, 'XTickLabel',[]); hold off; legend('location', 'northeast'); ax2 = nexttile; hold on; plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gt*Kt, freqs, 'Hz'))))); set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); yticks(-360:90:360); ylim([-200, 0]); xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); hold off; linkaxes([ax1,ax2],'x'); xlim([freqs(1), freqs(end)]); % #+name: fig:phase_gain_margin_model_plant % #+caption: Bode plot of the obtained Loop Gain $L(s)$ % #+attr_latex: :float nil % #+RESULTS: % [[file:figs/phase_gain_margin_model_plant.png]] % 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: xi = 0.03; Gr = 1/k*(s/w0/4 + 1)/(s^2/w0^2 + 2*xi*s/w0 + 1); % The obtained "real" loop gain is shown in Figure [[fig:phase_gain_margin_real_plant]]. % At a frequency little bit above 100Hz, the phase of the loop gain reaches -180 degrees while its magnitude is more than one which indicates instability. % It is confirmed by checking the stability of the closed loop system: isstable(feedback(Gr,Kt)) % #+RESULTS: % : 0 freqs = logspace(0, 3, 1000); figure; tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); ax1 = nexttile([2,1]); hold on; plot(freqs, abs(squeeze(freqresp(Gt*Kt, freqs, 'Hz'))), 'DisplayName', '$L(s)$'); plot(freqs, abs(squeeze(freqresp(Gr*Kt, freqs, 'Hz'))), 'DisplayName', '$L_r(s)$'); hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); ylabel('Magnitude'); set(gca, 'XTickLabel',[]); hold off; legend('location', 'northeast'); ax2 = nexttile; hold on; plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gt*Kt, freqs, 'Hz'))))); plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gr*Kt, freqs, 'Hz'))))); set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); yticks(-360:90:360); ylim([-200, 0]); xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); hold off; linkaxes([ax1,ax2],'x'); xlim([freqs(1), freqs(end)]); % #+name: fig:phase_gain_margin_real_plant % #+caption: Bode plots of $L(s)$ (loop gain corresponding the nominal plant) and $L_r(s)$ (loop gain corresponding to the real plant) % #+attr_latex: :float nil % #+RESULTS: % [[file:figs/phase_gain_margin_real_plant.png]] % Therefore, even a small change of the plant parameter renders the system unstable even though both the gain margin and the phase margin for the nominal plant are excellent. % This is due to the fact that the gain and phase margin are robustness indicators corresponding a *pure* change or gain or a *pure* change of phase but not a combination of both. % #+end_exampl % Let's now determine a new robustness indicator based on the Nyquist Stability Criteria. % #+begin_definition % - Nyquist Stability Criteria (for stable systems) :: % If the open-loop transfer function $L(s)$ is stable, then the closed-loop system will be unstable for any encirclement of the point $-1$ on the Nyquist plot. % - Nyquist Plot :: % The Nyquist plot shows the evolution of $L(j\omega)$ in the complex plane from $\omega = 0 \to \infty$. % #+end_definition % #+begin_seealso % For more information about the /general/ Nyquist Stability Criteria, you may want to look at [[https://www.youtube.com/watch?v=sof3meN96MA][this]] video. % #+end_seealso % From the Nyquist stability criteria, it is clear that we want $L(j\omega)$ to be as far as possible from the $-1$ point (called the /unstable point/) in the complex plane. % This minimum distance is called the *module margin*. % #+begin_definition % - Module Margin :: % The Module Margin $\Delta M$ is defined as the *minimum distance* between the point $-1$ and the loop gain $L(j\omega)$ in the complex plane. % #+end_definition % #+begin_exampl % A typical Nyquist plot is shown in Figure [[fig:module_margin_example]]. % The gain, phase and module margins are graphically shown to have an idea of what they represent. % Example Plant k = 1e6; % Stiffness [N/m] c = 4e2; % Damping [N/(m/s)] m = 10; % Mass [kg] G = 1/(m*s^2 + c*s + k); % Plant % Example Controller K = 14e8 * ... % Gain 1/(s^2) * ... % Double Integrator (1 + s/(2*pi*10/sqrt(8)))/(1 + s/(2*pi*10*sqrt(8))); % Lead L = G*K; L_resp = squeeze(freqresp(L, freqs, 'Hz')); % Module Margin Dm = min(abs(1 + L_resp)); % Phase Gain Margin [Gm, Pm, Wcg, Wcp] = margin(L); freqs = logspace(0, 3, 1000); figure; hold on; % Gain Margin plot([-1, -1/Gm], [0, 0], '-', 'DisplayName', sprintf('$\\Delta G = %.1f$', Gm)) % Phase Margin theta = -pi:0.01:-pi+Pm*pi/180; plot(cos(theta), sin(theta), '-', 'DisplayName', sprintf('$\\Delta \\phi = %.1f^o$', Pm)); % Module Margin theta = 0 : 0.01 : 2*pi; plot(Dm*cos(theta)-1, Dm*sin(theta), '-', 'DisplayName', sprintf('$\\Delta M = %.1f$', Dm)); % Nyquist Plot plot(real(L_resp), imag(L_resp), 'k-', 'DisplayName', '$L(j\omega)$') plot(-1, 0, 'k*', 'HandleVisibility', 'off'); hold off; xlabel('Real Axis'); ylabel('Imaginary Axis'); xlim([-1.5, 0.5]); ylim([-1, 1]); axis equal; legend('location', 'southeast');