Add rotation speed analysis

rotating_frame.org

@@ -1,22 +1,27 @@
 #+TITLE: Control in a rotating frame
+:DRAWER:
+#+STARTUP: overview
+
 #+HTML_HEAD: <link rel="stylesheet" type="text/css" href="css/htmlize.css"/>
 #+HTML_HEAD: <link rel="stylesheet" type="text/css" href="css/readtheorg.css"/>
 #+HTML_HEAD: <script src="js/jquery.min.js"></script>
 #+HTML_HEAD: <script src="js/bootstrap.min.js"></script>
 #+HTML_HEAD: <script type="text/javascript" src="js/jquery.stickytableheaders.min.js"></script>
 #+HTML_HEAD: <script type="text/javascript" src="js/readtheorg.js"></script>
+
 #+LATEX_CLASS: cleanreport
 #+LaTeX_CLASS_OPTIONS: [tocnp, secbreak, minted]
-#+STARTUP: overview
 #+LaTeX_HEADER: \usepackage{svg}
 #+LaTeX_HEADER: \newcommand{\authorFirstName}{Thomas}
 #+LaTeX_HEADER: \newcommand{\authorLastName}{Dehaeze}
 #+LaTeX_HEADER: \newcommand{\authorEmail}{dehaeze.thomas@gmail.com}
+
 #+PROPERTY: header-args:matlab  :session *MATLAB*
 #+PROPERTY: header-args:matlab+ :comments org
 #+PROPERTY: header-args:matlab+ :exports both
 #+PROPERTY: header-args:matlab+ :eval no-export
 #+PROPERTY: header-args:matlab+ :output-dir Figures
+:END:
 
 * Introduction
 The objective of this note it to highlight some control problems that arises when controlling the position of an object using actuators that are rotating with respect to a fixed reference frame.
@@ -902,7 +907,8 @@ We obtain the same result than the analytical case (figures [[fig:coupling_light
 #+RESULTS:
 [[file:Figures/coupling_ratio_light_heavy.png]]
 
-** Plant Control
+** Plant Control - SISO approach
+*** Plant identification
 The goal is to study the control problems due to the coupling that appears because of the rotation.
 
 #+begin_src matlab :exports none :results silent
@@ -1086,45 +1092,128 @@ However, when we look at the poles of the closed loop with a diagonal controller
 | -1.1837+0.0041777i |
 | -1.1837-0.0041777i |
 
-***
-
-
-
-
-
-
-
-
-
-
-
-#+begin_src matlab :results none :exports code
-  figure;
-  bode(Kll*Gvc('Du', 'fu'))
-#+end_src
-
-#+begin_src matlab :results none :exports code
-  sisotool(Gtvc('Du', 'fu'), Kll)
-#+end_src
-
-#+begin_src matlab :results none :exports code
-  figure;
-  bode(Kll*Gtvc('Du', 'fu'))
-#+end_src
-
+*** TODO Close loop performance
+First we compute the close loop transfer functions.
 #+begin_src matlab :results none :exports code
   Gcl = feedback(Gvc, K);
-#+end_src
-
-#+begin_src matlab :results none :exports code
   Gtcl = feedback(Gtvc, K);
 #+end_src
 
+
 #+begin_src matlab :results none :exports code
   figure;
   bode(Gcl, Gtcl)
 #+end_src
 
+*** Effect of rotation speed
+We first identify the system (voice coil and light mass) for multiple rotation speed.
+Then we compute the bode plot of the diagonal element (figure [[fig:Guu_ws]]) and of the coupling element (figure [[fig:Guv_ws]]).
+
+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)
+
+To stabilize the unstable pole, we need a control bandwidth of at least twice of frequency of the unstable pole.
+
+#+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 :exports none :results silent
+  freqs = logspace(-2, 2, 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.png" :var figsize="full-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:Guu_ws
+#+CAPTION: Diagonal term as a function of the rotation frequency
+#+RESULTS:
+[[file:Figures/Guu_ws.png]]
+
+#+begin_src matlab :exports none :results silent
+  freqs = logspace(-2, 2, 1000);
+
+  figure;
+  ax1 = subaxis(2,1,1);
+  hold on;
+  for i = 1:length(ws)
+    plot(freqs, abs(squeeze(freqresp(Gs{i}(1, 2), 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, 2), 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/Guv_ws.png" :var figsize="full-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:Guv_ws
+#+CAPTION: Couplin term as a function of the rotation frequency
+#+RESULTS:
+[[file:Figures/Guv_ws.png]]
+
+** TODO Plant Control - MIMO approach
 
 ** test