Add figure about RGA

2020-11-23 18:01:13 +01:00
parent cd38200cc8
commit 87a0d98e01
11 changed files with 16173 additions and 14962 deletions
--- a/index.org
+++ b/index.org
@@ -442,7 +442,6 @@ This is computed over the following frequencies.
      set(gca,'ColorOrderIndex',3)
      plot(freqs, Gr_jacobian(:,in_i), 'HandleVisibility', 'off');
  end
-  plot(freqs, 0.5*ones(size(freqs)), 'k--', 'DisplayName', 'Limit')
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  hold off;
  xlabel('Frequency (Hz)'); ylabel('Gershgorin Radii')
@@ -1382,6 +1381,10 @@ The analysis of the SVD control applied to the Stewart platform is performed in
  addpath('STEP');
 #+end_src

+#+begin_src matlab
+  freqs = logspace(-1, 2, 1000);
+#+end_src
+
 ** Jacobian                                                        :noexport:
 First, the position of the "joints" (points of force application) are estimated and the Jacobian computed.
 #+begin_src matlab :tangle no
@@ -1413,7 +1416,7 @@ First, the position of the "joints" (points of force application) are estimated

  J = [As' , cross(Ab, As)'];

-  save('./jacobian.mat', 'Aa', 'Ab', 'As', 'l', 'J');
+  save('stewart_platform/jacobian.mat', 'Aa', 'Ab', 'As', 'l', 'J');
 #+end_src

 ** Simscape Model - Parameters
@@ -1433,6 +1436,11 @@ Definition of spring parameters:
  cz = 0.025;
 #+end_src

+We suppose the sensor is perfectly positioned.
+#+begin_src matlab
+  sens_pos_error = zeros(3,1);
+#+end_src
+
 Gravity:
 #+begin_src matlab
  g = 0;
@@ -1440,7 +1448,7 @@ Gravity:

 We load the Jacobian (previously computed from the geometry):
 #+begin_src matlab
-  load('./jacobian.mat', 'Aa', 'Ab', 'As', 'l', 'J');
+  load('jacobian.mat', 'Aa', 'Ab', 'As', 'l', 'J');
 #+end_src

 We initialize other parameters:
@@ -1524,8 +1532,6 @@ The elements of the transfer matrix $\bm{G}$ corresponding to the transfer funct
 One can easily see that the system is strongly coupled.

 #+begin_src matlab :exports none
-  freqs = logspace(-1, 2, 1000);
-
  figure;

  % Magnitude
@@ -1563,6 +1569,21 @@ One can easily see that the system is strongly coupled.
 Consider the control architecture shown in Figure [[fig:plant_decouple_jacobian]].
 The Jacobian matrix is used to transform forces/torques applied on the top platform to the equivalent forces applied by each actuator.

+The Jacobian matrix is computed from the geometry of the platform (position and orientation of the actuators).
+
+#+begin_src matlab :exports results :results value table replace :tangle no
+  data2orgtable(J, {}, {}, ' %.3f ');
+#+end_src
+
+#+caption: Computed Jacobian Matrix
+#+RESULTS:
+|  0.811 |    0.0 | 0.584 | -0.018 | -0.008 |  0.025 |
+| -0.406 | -0.703 | 0.584 | -0.016 | -0.012 | -0.025 |
+| -0.406 |  0.703 | 0.584 |  0.016 | -0.012 |  0.025 |
+|  0.811 |    0.0 | 0.584 |  0.018 | -0.008 | -0.025 |
+| -0.406 | -0.703 | 0.584 |  0.002 |  0.019 |  0.025 |
+| -0.406 |  0.703 | 0.584 | -0.002 |  0.019 | -0.025 |
+
 #+begin_src latex :file plant_decouple_jacobian.pdf :tangle no :exports results
  \begin{tikzpicture}
    \node[block] (G) {$G_u$};
@@ -1633,6 +1654,7 @@ This can be verified below where only the real value of $G_u(\omega_c)$ is shown
  data2orgtable(real(evalfr(Gu, j*wc)), {}, {}, ' %.1f ');
 #+end_src

+#+caption: Real part of $G$ at the decoupling frequency $\omega_c$
 #+RESULTS:
 |    4.4 |   -2.1 |   -2.1 |    4.4 |  -2.4 |   -2.4 |
 |   -0.2 |   -3.9 |    3.9 |    0.2 |  -3.8 |    3.8 |
@@ -1651,6 +1673,32 @@ First, the Singular Value Decomposition of $H_1$ is performed:
  [U,~,V] = svd(H1);
 #+end_src

+#+begin_src matlab :exports results :results value table replace :tangle no
+  data2orgtable(U, {}, {}, ' %.1g ');
+#+end_src
+
+#+caption: Obtained matrix $U$
+#+RESULTS:
+| -0.005 |  7e-06 |  6e-11 | -3e-06 |     -1 |    0.1 |
+| -7e-06 | -0.005 | -9e-09 | -5e-09 |   -0.1 |     -1 |
+|  4e-08 | -2e-10 | -6e-11 |     -1 |  3e-06 | -3e-07 |
+| -0.002 |     -1 | -5e-06 |  2e-10 | 0.0006 |  0.005 |
+|      1 | -0.002 | -1e-08 |  2e-08 | -0.005 | 0.0006 |
+| -4e-09 |  5e-06 |     -1 |  6e-11 | -2e-09 | -1e-08 |
+
+#+begin_src matlab :exports results :results value table replace :tangle no
+  data2orgtable(V, {}, {}, ' %.1g ');
+#+end_src
+
+#+caption: Obtained matrix $V$
+#+RESULTS:
+| -0.2 |   0.5 | -0.4 | -0.4 |   -0.6 | -0.2 |
+| -0.3 |   0.5 |  0.4 | -0.4 |    0.5 |  0.3 |
+| -0.3 |  -0.5 | -0.4 | -0.4 |    0.4 | -0.4 |
+| -0.2 |  -0.5 |  0.4 | -0.4 |   -0.5 |  0.3 |
+|  0.6 | -0.06 | -0.4 | -0.4 |    0.1 |  0.6 |
+|  0.6 |  0.06 |  0.4 | -0.4 | -0.006 | -0.6 |
+
 The obtained matrices $U$ and $V$ are used to decouple the system as shown in Figure [[fig:plant_decouple_svd]].

 #+begin_src latex :file plant_decouple_svd.pdf :tangle no :exports results
@@ -1694,10 +1742,6 @@ The "Gershgorin Radii" of a matrix $S$ is defined by:
 \[ \zeta_i(j\omega) = \frac{\sum\limits_{j\neq i}|S_{ij}(j\omega)|}{|S_{ii}(j\omega)|} \]

 This is computed over the following frequencies.
-#+begin_src matlab
-  freqs = logspace(-2, 2, 1000); % [Hz]
-#+end_src
-
 #+begin_src matlab :exports none
  % Gershgorin Radii for the coupled plant:
  Gr_coupled = zeros(length(freqs), size(Gu,2));
@@ -1735,7 +1779,6 @@ This is computed over the following frequencies.
      set(gca,'ColorOrderIndex',3)
      plot(freqs, Gr_jacobian(:,in_i), 'HandleVisibility', 'off');
  end
-  plot(freqs, 0.5*ones(size(freqs)), 'k--', 'DisplayName', 'Limit')
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  hold off;
  xlabel('Frequency (Hz)'); ylabel('Gershgorin Radii')
@@ -1752,14 +1795,105 @@ This is computed over the following frequencies.
 #+RESULTS:
 [[file:figs/simscape_model_gershgorin_radii.png]]

+** Verification of the decoupling using the "Relative Gain Array"
+The relative gain array (RGA) is defined as:
+\begin{equation}
+  \Lambda\big(G(s)\big) = G(s) \times \big( G(s)^{-1} \big)^T
+\end{equation}
+where $\times$ denotes an element by element multiplication and $G(s)$ is an $n \times n$ square transfer matrix.
+
+The obtained RGA elements are shown in Figure [[fig:simscape_model_rga]].
+
+#+begin_src matlab :exports none
+  % Relative Gain Array for the coupled plant:
+  RGA_coupled = zeros(length(freqs), size(Gu,1), size(Gu,2));
+  Gu_inv = inv(Gu);
+  for f_i = 1:length(freqs)
+    RGA_coupled(f_i, :, :) = abs(evalfr(Gu, j*2*pi*freqs(f_i)).*evalfr(Gu_inv, j*2*pi*freqs(f_i))');
+  end
+
+  % Relative Gain Array for the decoupled plant using SVD:
+  RGA_svd = zeros(length(freqs), size(Gsvd,1), size(Gsvd,2));
+  Gsvd_inv = inv(Gsvd);
+  for f_i = 1:length(freqs)
+    RGA_svd(f_i, :, :) = abs(evalfr(Gsvd, j*2*pi*freqs(f_i)).*evalfr(Gsvd_inv, j*2*pi*freqs(f_i))');
+  end
+
+  % Relative Gain Array for the decoupled plant using the Jacobian:
+  RGA_x = zeros(length(freqs), size(Gx,1), size(Gx,2));
+  Gx_inv = inv(Gx);
+  for f_i = 1:length(freqs)
+    RGA_x(f_i, :, :) = abs(evalfr(Gx, j*2*pi*freqs(f_i)).*evalfr(Gx_inv, j*2*pi*freqs(f_i))');
+  end
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None');
+
+  ax1 = nexttile;
+  hold on;
+  for i_in = 1:6
+      for i_out = [1:i_in-1, i_in+1:6]
+          plot(freqs, RGA_svd(:, i_out, i_in), '--', 'color', [0 0 0 0.2], ...
+               'HandleVisibility', 'off');
+      end
+  end
+  plot(freqs, RGA_svd(:, 1, 2), '--', 'color', [0 0 0 0.2], ...
+       'DisplayName', '$RGA_{SVD}(i,j),\ i \neq j$');
+
+  plot(freqs, RGA_svd(:, 1, 1), 'k-', ...
+       'DisplayName', '$RGA_{SVD}(i,i)$');
+  for ch_i = 1:6
+    plot(freqs, RGA_svd(:, ch_i, ch_i), 'k-', ...
+         'HandleVisibility', 'off');
+  end
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Magnitude'); xlabel('Frequency [Hz]');
+  legend('location', 'southwest');
+
+  ax2 = nexttile;
+  hold on;
+  for i_in = 1:6
+      for i_out = [1:i_in-1, i_in+1:6]
+          plot(freqs, RGA_x(:, i_out, i_in), '--', 'color', [0 0 0 0.2], ...
+               'HandleVisibility', 'off');
+      end
+  end
+  plot(freqs, RGA_x(:, 1, 2), '--', 'color', [0 0 0 0.2], ...
+       'DisplayName', '$RGA_{X}(i,j),\ i \neq j$');
+
+  plot(freqs, RGA_x(:, 1, 1), 'k-', ...
+       'DisplayName', '$RGA_{X}(i,i)$');
+  for ch_i = 1:6
+    plot(freqs, RGA_x(:, ch_i, ch_i), 'k-', ...
+         'HandleVisibility', 'off');
+  end
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); set(gca, 'YTickLabel',[]);
+  legend('location', 'southwest');
+
+  linkaxes([ax1,ax2],'y');
+  ylim([1e-5, 1e1]);
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/simscape_model_rga.pdf', 'width', 'wide', 'height', 'tall');
+#+end_src
+
+#+name: fig:simscape_model_rga
+#+caption: Obtained norm of RGA elements for the SVD decoupled plant and the Jacobian decoupled plant
+#+RESULTS:
+[[file:figs/simscape_model_rga.png]]
+
 ** Obtained Decoupled Plants
 <<sec:stewart_decoupled_plant>>

 The bode plot of the diagonal and off-diagonal elements of $G_{SVD}$ are shown in Figure [[fig:simscape_model_decoupled_plant_svd]].

 #+begin_src matlab :exports none
-  freqs = logspace(-1, 2, 1000);
-
  figure;
  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');

@@ -1812,8 +1946,6 @@ The bode plot of the diagonal and off-diagonal elements of $G_{SVD}$ are shown i
 Similarly, the bode plots of the diagonal elements and off-diagonal elements of the decoupled plant $G_x(s)$ using the Jacobian are shown in Figure [[fig:simscape_model_decoupled_plant_jacobian]].

 #+begin_src matlab :exports none
-  freqs = logspace(-1, 2, 1000);
-
  figure;
  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');

@@ -1956,8 +2088,6 @@ $G_0$ is tuned such that the crossover frequency corresponding to the diagonal t
 The obtained diagonal elements of the loop gains are shown in Figure [[fig:stewart_comp_loop_gain_diagonal]].

 #+begin_src matlab :exports none
-  freqs = logspace(-1, 2, 1000);
-
  figure;
  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');

@@ -2038,8 +2168,6 @@ Let's first verify the stability of the closed-loop systems:
 The obtained transmissibility in Open-loop, for the centralized control as well as for the SVD control are shown in Figure [[fig:stewart_platform_simscape_cl_transmissibility]].

 #+begin_src matlab :exports results
-  freqs = logspace(-2, 2, 1000);
-
  figure;
  tiledlayout(2, 2, 'TileSpacing', 'None', 'Padding', 'None');

@@ -2102,6 +2230,159 @@ The obtained transmissibility in Open-loop, for the centralized control as well
 #+RESULTS:
 [[file:figs/stewart_platform_simscape_cl_transmissibility.png]]

+** Small error on the sensor location                             :no_export:
+Let's now consider a small position error of the sensor:
+#+begin_src matlab
+  sens_pos_error = [105 5 -1]*1e-3; % [m]
+#+end_src
+
+The system is identified again:
+#+begin_src matlab :exports none
+  %% Name of the Simulink File
+  mdl = 'drone_platform';
+
+  %% Input/Output definition
+  clear io; io_i = 1;
+  io(io_i) = linio([mdl, '/Dw'],              1, 'openinput');  io_i = io_i + 1; % Ground Motion
+  io(io_i) = linio([mdl, '/V-T'],             1, 'openinput');  io_i = io_i + 1; % Actuator Forces
+  io(io_i) = linio([mdl, '/Inertial Sensor'], 1, 'openoutput'); io_i = io_i + 1; % Top platform acceleration
+
+  G = linearize(mdl, io);
+  G.InputName  = {'Dwx', 'Dwy', 'Dwz', 'Rwx', 'Rwy', 'Rwz', ...
+                  'F1', 'F2', 'F3', 'F4', 'F5', 'F6'};
+  G.OutputName = {'Ax', 'Ay', 'Az', 'Arx', 'Ary', 'Arz'};
+
+  % Plant
+  Gu = G(:, {'F1', 'F2', 'F3', 'F4', 'F5', 'F6'});
+  % Disturbance dynamics
+  Gd = G(:, {'Dwx', 'Dwy', 'Dwz', 'Rwx', 'Rwy', 'Rwz'});
+#+end_src
+
+#+begin_src matlab
+  Gx = Gu*inv(J');
+  Gx.InputName  = {'Fx', 'Fy', 'Fz', 'Mx', 'My', 'Mz'};
+#+end_src
+
+#+begin_src matlab
+  Gsvd = inv(U)*Gu*inv(V');
+#+end_src
+
+#+begin_src matlab :exports none
+  % Gershgorin Radii for the coupled plant:
+  Gr_coupled = zeros(length(freqs), size(Gu,2));
+  H = abs(squeeze(freqresp(Gu, freqs, 'Hz')));
+  for out_i = 1:size(Gu,2)
+      Gr_coupled(:, out_i) = squeeze((sum(H(out_i,:,:)) - H(out_i,out_i,:))./H(out_i, out_i, :));
+  end
+
+  % Gershgorin Radii for the decoupled plant using SVD:
+  Gr_decoupled = zeros(length(freqs), size(Gsvd,2));
+  H = abs(squeeze(freqresp(Gsvd, freqs, 'Hz')));
+  for out_i = 1:size(Gsvd,2)
+      Gr_decoupled(:, out_i) = squeeze((sum(H(out_i,:,:)) - H(out_i,out_i,:))./H(out_i, out_i, :));
+  end
+
+  % Gershgorin Radii for the decoupled plant using the Jacobian:
+  Gr_jacobian = zeros(length(freqs), size(Gx,2));
+  H = abs(squeeze(freqresp(Gx, freqs, 'Hz')));
+  for out_i = 1:size(Gx,2)
+      Gr_jacobian(:, out_i) = squeeze((sum(H(out_i,:,:)) - H(out_i,out_i,:))./H(out_i, out_i, :));
+  end
+#+end_src
+
+#+begin_src matlab :exports results
+  figure;
+  hold on;
+  plot(freqs, Gr_coupled(:,1), 'DisplayName', 'Coupled');
+  plot(freqs, Gr_decoupled(:,1), 'DisplayName', 'SVD');
+  plot(freqs, Gr_jacobian(:,1), 'DisplayName', 'Jacobian');
+  for in_i = 2:6
+      set(gca,'ColorOrderIndex',1)
+      plot(freqs, Gr_coupled(:,in_i), 'HandleVisibility', 'off');
+      set(gca,'ColorOrderIndex',2)
+      plot(freqs, Gr_decoupled(:,in_i), 'HandleVisibility', 'off');
+      set(gca,'ColorOrderIndex',3)
+      plot(freqs, Gr_jacobian(:,in_i), 'HandleVisibility', 'off');
+  end
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  hold off;
+  xlabel('Frequency (Hz)'); ylabel('Gershgorin Radii')
+  legend('location', 'northwest');
+  ylim([1e-3, 1e3]);
+#+end_src
+
+#+begin_src matlab
+  L_cen = K_cen*Gx;
+  G_cen = feedback(G, pinv(J')*K_cen, [7:12], [1:6]);
+#+end_src
+
+#+begin_src matlab
+  L_svd = K_svd*Gsvd;
+  G_svd = feedback(G, inv(V')*K_svd*inv(U), [7:12], [1:6]);
+#+end_src
+
+#+begin_src matlab :results output replace text
+  isstable(G_cen)
+#+end_src
+
+#+begin_src matlab :results output replace text
+  isstable(G_svd)
+#+end_src
+
+#+begin_src matlab :exports results
+  figure;
+  tiledlayout(2, 2, 'TileSpacing', 'None', 'Padding', 'None');
+
+  ax1 = nexttile;
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G(    'Ax', 'Dwx')/s^2, freqs, 'Hz'))), 'DisplayName', 'Open-Loop');
+  plot(freqs, abs(squeeze(freqresp(G_cen('Ax', 'Dwx')/s^2, freqs, 'Hz'))), 'DisplayName', 'Centralized');
+  plot(freqs, abs(squeeze(freqresp(G_svd('Ax', 'Dwx')/s^2, freqs, 'Hz'))), '--', 'DisplayName', 'SVD');
+  set(gca,'ColorOrderIndex',1)
+  plot(freqs, abs(squeeze(freqresp(G(    'Ay', 'Dwy')/s^2, freqs, 'Hz'))), 'HandleVisibility', 'off');
+  plot(freqs, abs(squeeze(freqresp(G_cen('Ay', 'Dwy')/s^2, freqs, 'Hz'))), 'HandleVisibility', 'off');
+  plot(freqs, abs(squeeze(freqresp(G_svd('Ay', 'Dwy')/s^2, freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('$D_x/D_{w,x}$, $D_y/D_{w, y}$'); set(gca, 'XTickLabel',[]);
+  legend('location', 'southwest');
+
+  ax2 = nexttile;
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G(    'Az', 'Dwz')/s^2, freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_cen('Az', 'Dwz')/s^2, freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_svd('Az', 'Dwz')/s^2, freqs, 'Hz'))), '--');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('$D_z/D_{w,z}$'); set(gca, 'XTickLabel',[]);
+
+  ax3 = nexttile;
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G(    'Arx', 'Rwx')/s^2, freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_cen('Arx', 'Rwx')/s^2, freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_svd('Arx', 'Rwx')/s^2, freqs, 'Hz'))), '--');
+  set(gca,'ColorOrderIndex',1)
+  plot(freqs, abs(squeeze(freqresp(G(    'Ary', 'Rwy')/s^2, freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_cen('Ary', 'Rwy')/s^2, freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_svd('Ary', 'Rwy')/s^2, freqs, 'Hz'))), '--');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('$R_x/R_{w,x}$, $R_y/R_{w,y}$');  xlabel('Frequency [Hz]');
+
+  ax4 = nexttile;
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G(    'Arz', 'Rwz')/s^2, freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_cen('Arz', 'Rwz')/s^2, freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_svd('Arz', 'Rwz')/s^2, freqs, 'Hz'))), '--');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('$R_z/R_{w,z}$');  xlabel('Frequency [Hz]');
+
+  linkaxes([ax1,ax2,ax3,ax4],'xy');
+  xlim([freqs(1), freqs(end)]);
+  ylim([1e-3, 1e2]);
+#+end_src
+
 * Stewart Platform - Analytical Model                               :noexport:
 :PROPERTIES:
 :header-args:matlab+: :tangle stewart_platform/analytical_model.m