309 lines
10 KiB
Matlab
309 lines
10 KiB
Matlab
%% 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
|
|
% <<sec:sensitivity_transfer_functions>>
|
|
|
|
% 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
|
|
% <<sec: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');
|