Add SVD analysis

rotating_frame.org

@@ -417,14 +417,11 @@ Then, we compute the bode plot of the direct term and coupling term for multiple
 Then we compare the result between voice coil and piezoelectric actuators.
 
 *** Voice coil actuator
-
 #+begin_src matlab :exports none :results silent
   m = mlight;
   k = kvc;
-#+end_src
 
-#+begin_src matlab :exports none :results silent
+  ws = linspace(0, 2*pi, 5); % Rotation speed vector [rad/s]
-  ws = [0, 1, 10, 60]*2*pi/60; % [rmp]
 
   Gs = {zeros(1, length(ws))};
   Gcs = {zeros(1, length(ws))};
@@ -443,19 +440,18 @@ Then we compare the result between voice coil and piezoelectric actuators.
   ax1 = subaxis(2,1,1);
   hold on;
   for i = 1:length(ws)
-    plot(freqs, abs(squeeze(freqresp(Gs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
+    plot(freqs, abs(squeeze(freqresp(Gs{i}, freqs, 'Hz'))));
   end
   hold off;
   xlim([freqs(1), freqs(end)]);
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   set(gca, 'XTickLabel',[]);
   ylabel('Magnitude [m/N]');
-  legend('Location', 'southwest');
 
   ax2 = subaxis(2,1,2);
   hold on;
   for i = 1:length(ws)
-    plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}, freqs, 'Hz'))));
+    plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
   end
   hold off;
   yticks(-180:90:180);
@@ -463,6 +459,7 @@ Then we compare the result between voice coil and piezoelectric actuators.
   xlim([freqs(1), freqs(end)]);
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
   xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  legend('Location', 'southwest');
   linkaxes([ax1,ax2],'x');
 #+end_src
 
@@ -484,19 +481,18 @@ Then we compare the result between voice coil and piezoelectric actuators.
   ax1 = subaxis(2,1,1);
   hold on;
   for i = 1:length(ws)
-    plot(freqs, abs(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
+    plot(freqs, abs(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))));
   end
   hold off;
   xlim([freqs(1), freqs(end)]);
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   set(gca, 'XTickLabel',[]);
   ylabel('Magnitude [m/N]');
-  legend('Location', 'southwest');
 
   ax2 = subaxis(2,1,2);
   hold on;
   for i = 1:length(ws)
-    plot(freqs, 180/pi*angle(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))));
+    plot(freqs, 180/pi*angle(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
   end
   hold off;
   yticks(-180:90:180);
@@ -504,6 +500,7 @@ Then we compare the result between voice coil and piezoelectric actuators.
   xlim([freqs(1), freqs(end)]);
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
   xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  legend('Location', 'southwest');
   linkaxes([ax1,ax2],'x');
 #+end_src
 
@@ -519,14 +516,11 @@ Then we compare the result between voice coil and piezoelectric actuators.
 [[file:Figures/Gc_ws_vc.png]]
 
 *** Piezoelectric actuator
-
 #+begin_src matlab :exports none :results silent
   m = mlight;
   k = kpz;
-#+end_src
 
-#+begin_src matlab :exports none :results silent
+  ws = linspace(0, 2*pi, 5); % Rotation speed vector [rad/s]
-  ws = [0, 1, 10, 60]*2*pi/60; % [rmp]
 
   Gs = {zeros(1, length(ws))};
   Gcs = {zeros(1, length(ws))};
@@ -545,19 +539,18 @@ Then we compare the result between voice coil and piezoelectric actuators.
   ax1 = subaxis(2,1,1);
   hold on;
   for i = 1:length(ws)
-    plot(freqs, abs(squeeze(freqresp(Gs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
+    plot(freqs, abs(squeeze(freqresp(Gs{i}, freqs, 'Hz'))));
   end
   hold off;
   xlim([freqs(1), freqs(end)]);
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   set(gca, 'XTickLabel',[]);
   ylabel('Magnitude [m/N]');
-  legend('Location', 'southwest');
 
   ax2 = subaxis(2,1,2);
   hold on;
   for i = 1:length(ws)
-    plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}, freqs, 'Hz'))));
+    plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
   end
   hold off;
   yticks(-180:90:180);
@@ -565,6 +558,7 @@ Then we compare the result between voice coil and piezoelectric actuators.
   xlim([freqs(1), freqs(end)]);
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
   xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  legend('Location', 'southwest');
   linkaxes([ax1,ax2],'x');
 #+end_src
 
@@ -584,19 +578,18 @@ Then we compare the result between voice coil and piezoelectric actuators.
   ax1 = subaxis(2,1,1);
   hold on;
   for i = 1:length(ws)
-    plot(freqs, abs(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
+    plot(freqs, abs(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))));
   end
   hold off;
   xlim([freqs(1), freqs(end)]);
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   set(gca, 'XTickLabel',[]);
   ylabel('Magnitude [m/N]');
-  legend('Location', 'southwest');
 
   ax2 = subaxis(2,1,2);
   hold on;
   for i = 1:length(ws)
-    plot(freqs, 180/pi*angle(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))));
+    plot(freqs, 180/pi*angle(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
   end
   hold off;
   yticks(-180:90:180);
@@ -604,6 +597,7 @@ Then we compare the result between voice coil and piezoelectric actuators.
   xlim([freqs(1), freqs(end)]);
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
   xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  legend('Location', 'southwest');
   linkaxes([ax1,ax2],'x');
 #+end_src
 
@@ -627,10 +621,11 @@ As the rotation speed increases, one of the two resonant frequency goes to lower
 
 The poles of the coupling terms are the same as the poles of the diagonal terms. The magnitude of the coupling terms are increasing with the rotation speed.
 
+#+begin_important
 As shown in the previous figures, the system with voice coil is much more sensitive to rotation speed.
+#+end_important
 
 *** Campbell diagram
-
 The poles of the system are computed for multiple values of the rotation frequency. To simplify the computation of the poles, we add some damping to the system.
 
 #+begin_src matlab :results silent :exports code
@@ -638,12 +633,12 @@ The poles of the system are computed for multiple values of the rotation frequen
   k = kvc;
   c = 0.1*sqrt(k*m);
 
-  ws = linspace(0, 10, 100); % [rad/s]
+  wsvc = linspace(0, 10, 100); % [rad/s]
 
-  polesvc = zeros(2, length(ws));
+  polesvc = zeros(2, length(wsvc));
 
-  for i = 1:length(ws)
+  for i = 1:length(wsvc)
-    polei = pole(1/((m*s^2 + c*s + (k - m*ws(i)^2))^2 + (2*m*ws(i)*s)^2));
+    polei = pole(1/((m*s^2 + c*s + (k - m*wsvc(i)^2))^2 + (2*m*wsvc(i)*s)^2));
     polesvc(:, i) = sort(polei(imag(polei) > 0));
   end
 #+end_src
@@ -653,12 +648,12 @@ The poles of the system are computed for multiple values of the rotation frequen
   k = kpz;
   c = 0.1*sqrt(k*m);
 
-  ws = linspace(0, 1000, 100); % [rad/s]
+  wspz = linspace(0, 1000, 100); % [rad/s]
 
-  polespz = zeros(2, length(ws));
+  polespz = zeros(2, length(wspz));
 
-  for i = 1:length(ws)
+  for i = 1:length(wspz)
-    polei = pole(1/((m*s^2 + c*s + (k - m*ws(i)^2))^2 + (2*m*ws(i)*s)^2));
+    polei = pole(1/((m*s^2 + c*s + (k - m*wspz(i)^2))^2 + (2*m*wspz(i)*s)^2));
     polespz(:, i) = sort(polei(imag(polei) > 0));
   end
 #+end_src
@@ -674,22 +669,22 @@ For the voice coil (figure [[fig:poles_w_vc]]), the system is unstable when the
   % Amplitude
   ax1 = subplot(1,2,1);
   hold on;
-  plot(ws, real(polesvc(1, :)))
+  plot(wsvc, real(polesvc(1, :)))
   set(gca,'ColorOrderIndex',1)
-  plot(ws, real(polesvc(2, :)))
+  plot(wsvc, real(polesvc(2, :)))
-  plot(ws, zeros(size(ws)), 'k--')
+  plot(wsvc, zeros(size(wsvc)), 'k--')
   hold off;
   xlabel('Rotation Frequency [rad/s]');
   ylabel('Pole Real Part');
   ax2 = subplot(1,2,2);
   hold on;
-  plot(ws, imag(polesvc(1, :)))
+  plot(wsvc, imag(polesvc(1, :)))
   set(gca,'ColorOrderIndex',1)
-  plot(ws, -imag(polesvc(1, :)))
+  plot(wsvc, -imag(polesvc(1, :)))
   set(gca,'ColorOrderIndex',1)
-  plot(ws, imag(polesvc(2, :)))
+  plot(wsvc, imag(polesvc(2, :)))
   set(gca,'ColorOrderIndex',1)
-  plot(ws, -imag(polesvc(2, :)))
+  plot(wsvc, -imag(polesvc(2, :)))
   hold off;
   xlabel('Rotation Frequency [rad/s]');
   ylabel('Pole Imaginary Part');
@@ -712,22 +707,22 @@ For the voice coil (figure [[fig:poles_w_vc]]), the system is unstable when the
   % Amplitude
   ax1 = subplot(1,2,1);
   hold on;
-  plot(ws, real(polespz(1, :)))
+  plot(wspz, real(polespz(1, :)))
   set(gca,'ColorOrderIndex',1)
-  plot(ws, real(polespz(2, :)))
+  plot(wspz, real(polespz(2, :)))
-  plot(ws, zeros(size(ws)), 'k--')
+  plot(wspz, zeros(size(wspz)), 'k--')
   hold off;
   xlabel('Rotation Frequency [rad/s]');
   ylabel('Pole Real Part');
   ax2 = subplot(1,2,2);
   hold on;
-  plot(ws, imag(polespz(1, :)))
+  plot(wspz, imag(polespz(1, :)))
   set(gca,'ColorOrderIndex',1)
-  plot(ws, -imag(polespz(1, :)))
+  plot(wspz, -imag(polespz(1, :)))
   set(gca,'ColorOrderIndex',1)
-  plot(ws, imag(polespz(2, :)))
+  plot(wspz, imag(polespz(2, :)))
   set(gca,'ColorOrderIndex',1)
-  plot(ws, -imag(polespz(2, :)))
+  plot(wspz, -imag(polespz(2, :)))
   hold off;
   xlabel('Rotation Frequency [rad/s]');
   ylabel('Pole Imaginary Part');
@@ -1058,6 +1053,8 @@ We can then plot the same loop gain with the rotating system using the same cont
   <<plt-matlab>>
 #+end_src
 
+#+NAME: fig:Gtvc_loop_gain
+#+CAPTION: Loop gain with the rotating system
 #+RESULTS:
 [[file:Figures/Gtvc_loop_gain.png]]
 
@@ -1157,7 +1154,7 @@ Then we compute the bode plot of the diagonal element (figure [[fig:Guu_ws]]) an
 As the rotation frequency increases:
 - one pole goes to lower frequencies while the other goes to higher frequencies
 - one zero appears between the two poles
-- [ ] the zero disappears when $\omega > \sqrt{\frac{k}{m}}$ and the low frequency pole becomes unstable (positive real part)
+- the zero disappears when $\omega > \sqrt{\frac{k}{m}}$ and the low frequency pole becomes unstable (positive real part)
 
 To stabilize the unstable pole, we need a control bandwidth of at least twice of frequency of the unstable pole.
 
@@ -1258,7 +1255,179 @@ To stabilize the unstable pole, we need a control bandwidth of at least twice of
 #+RESULTS:
 [[file:Figures/Guv_ws.png]]
 
-** TODO Plant Control - MIMO approach
+Then, we can look at the same plots for the piezoelectric actuator (figure [[fig:Guu_ws_pz]]). The effect of the rotation frequency has very little effect on the dynamics of the system to control.
+
+#+begin_src matlab :exports none :results silent
+  ws = linspace(0, 2*pi, 5); % Rotation speed vector [rad/s]
+  m = mlight; % mass of the sample [kg]
+
+  kTuv = kpz;
+  cTuv = 0.1*sqrt(kTuv*m);
+
+  Gs = {zeros(1, length(ws))};
+
+  for i = 1:length(ws)
+    w = ws(i);
+    Gs{i} = linearize(mdl, io, 0.1);
+  end
+#+end_src
+
+#+begin_src matlab :exports none :results silent
+  freqs = logspace(2, 3, 1000);
+
+  figure;
+  ax1 = subaxis(2,1,1);
+  hold on;
+  for i = 1:length(ws)
+    plot(freqs, abs(squeeze(freqresp(Gs{i}(1, 1), freqs, 'Hz'))));
+  end
+  hold off;
+  xlim([freqs(1), freqs(end)]);
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Magnitude [m/N]');
+
+  ax2 = subaxis(2,1,2);
+  hold on;
+  for i = 1:length(ws)
+    plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}(1, 1), freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
+  end
+  hold off;
+  yticks(-180:90:180);
+  ylim([-180 180]);
+  xlim([freqs(1), freqs(end)]);
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  legend('Location', 'northeast');
+  linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/Guu_ws_pz.png" :var figsize="full-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:Guu_ws_pz
+#+CAPTION: Diagonal term as a function of the rotation frequency
+#+RESULTS:
+[[file:Figures/Guu_ws_pz.png]]
+
+** Plant Control - MIMO approach
+*** TODO Analysis - SVD
+\[ G = U \Sigma V^H \]
+
+With:
+- $\Sigma$ is an $2 \times 2$ matrix with 2 non-negative *singular values* $\sigma_i$, arranged in descending order along its main diagonal
+- $U$ is an $2 \times 2$ unitary matrix. The columns vectors of $U$, denoted $u_i$, represent the *output directions* of the plant. They are orthonomal
+- $V$ is an $2 \times 2$ unitary matrix. The columns vectors of $V$, denoted $v_i$, represent the *input directions* of the plant. They are orthonomal
+
+We first look at the evolution of the singular values as a function of frequency (figure [[fig:G_sigma]]).
+#+begin_src matlab :exports none :results silent
+  ws = linspace(0, 2*pi, 5); % Rotation speed vector [rad/s]
+  m = mlight; % mass of the sample [kg]
+
+  kTuv = kvc;
+  cTuv = 0.1*sqrt(kTuv*m);
+
+  Gs = {zeros(1, length(ws))};
+
+  for i = 1:length(ws)
+    w = ws(i);
+    Gs{i} = linearize(mdl, io, 0.1);
+  end
+#+end_src
+
+#+begin_src matlab :results silent :exports none
+  freqs = logspace(-2, 2, 1000);
+
+  figure;
+  hold on;
+  for i = 1:length(ws)
+    plot(freqs, abs(squeeze(freqresp(Gs{i}(1,1), freqs, 'Hz'))), '-');
+  end
+  hold off;
+  set(gca,'xscale','log'); set(gca,'yscale','log');
+#+end_src
+
+#+begin_src matlab :results silent :exports none
+  freqs = logspace(-2, 1, 1000);
+
+  figure;
+  hold on;
+  for i = 1:length(ws)
+    sv = sigma(Gs{i}, 2*pi*freqs);
+    set(gca,'ColorOrderIndex',i)
+    plot(freqs, sv(1, :), 'DisplayName', sprintf('w = %.0f rpm', ws(i)*60/2/pi));
+    set(gca,'ColorOrderIndex',i)
+    plot(freqs, sv(2, :), '--', 'HandleVisibility', 'off');
+  end
+  hold off;
+  set(gca,'xscale','log'); set(gca,'yscale','log');
+  legend('location', 'southwest');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/G_sigma.png" :var figsize="full-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:G_sigma
+#+CAPTION: Evolution of the singular values with frequency
+#+RESULTS:
+[[file:Figures/G_sigma.png]]
+
+We compute
+#+begin_src matlab :results silent :exports code
+  [U,S,V] = svd(freqresp(Gtvc, 2*pi*10));
+#+end_src
+
+#+begin_src matlab :results output :exports results
+  U, S, V
+#+end_src
+
+#+RESULTS:
+#+begin_example
+U, S, V
+U =
+         -0.707101109012986 -   0.00283224868340902i         -0.707104254409621 -   0.00189034277692295i
+        0.00283224868340845 -     0.707101109012987i       -0.00189034277692242 +      0.70710425440962i
+S =
+      9.01532756059351e-06                         0
+                         0      6.01714794171208e-06
+V =
+          0.707106781186547 +                     0i          0.707106781186548 +                     0i
+      -1.57009245868378e-16 +     0.707106781186548i       1.57009245868377e-16 -     0.707106781186547i
+#+end_example
+
+The input and output directions are related through the singular values
+\[ G v_i = \sigma_i u_i \]
+
+So, if we consider an input in the direction $v_i$, then the output is in the direction $u_i$. Furthermore, since $\normtwo{v_i}=1$ and $\normtwo{u_i}=1$, we see that *the singular value $\sigma_i$ directly gives the gain of the matrix $G$ in this direction*.
+
+#+begin_src matlab
+  freqresp(Gtvc, 2*pi*10)*V(:, 1)
+#+end_src
+
+#+RESULTS:
+| -6.3747e-06-2.5534e-08i |
+| 2.5534e-08-6.3747e-06i  |
+
+#+begin_src matlab
+  S(1)*U(:, 1)
+#+end_src
+
+#+RESULTS:
+| -6.3747e-06-2.5534e-08i |
+| 2.5534e-08-6.3747e-06i  |
+
+*** Closed loop SVD
+#+begin_src matlab :results silent :exports none
+  figure;
+  sigma(Tvc, Ttvc)
+#+end_src
+
 ** test
 
 #+begin_src matlab :exports none :results silent
@@ -1549,4 +1718,4 @@ Finally, we run the linearization.
   <<sec:control>>
 ** Measurement in the fixed reference frame
 * Bibliography                                                       :ignore:
-#+BIBLIOGRAPHY: /home/tdehaeze/MEGA/These/Ressources/references.bib plain option:-a option:-noabstract option:-nokeywords option:-noheader option:-nofooter option:-nobibsource limit:t
+# #+BIBLIOGRAPHY: /home/tdehaeze/MEGA/These/Ressources/references.bib plain option:-a option:-noabstract option:-nokeywords option:-noheader option:-nofooter option:-nobibsource limit:t