Add some analysis on the coupling. Add uv_to_xy simulink blocs

rotating_frame.org

@@ -779,7 +779,7 @@ The loop gain is $L = G K$.
   :END:
   <<sec:simscape>>
 
-** Initialize
+** Initialization
 #+begin_src matlab :exports none :results silent :noweb yes
   <<matlab-init>>
   load('./mat/parameters.mat');
@@ -789,7 +789,6 @@ The loop gain is $L = G K$.
   open rotating_frame.slx
 #+end_src
 
-** Parameter for the Simscape simulations
 First we define the parameters that must be defined in order to run the Simscape simulation.
 #+begin_src matlab :exports code :results silent
   w = 2*pi; % Rotation speed [rad/s]
@@ -873,8 +872,10 @@ Then we identify the system with an heavy mass and low speed.
 #+end_src
 
 ** Coupling ratio between $f_{uv}$ and $d_{uv}$
+In order to validate the equations written, we can compute the coupling ratio using the simscape model and compare with the equations.
+
+From the previous identification, we plot the coupling ratio in both case (figure [[fig:coupling_ratio_light_heavy]]).
 
-From the previous identification, we plot the coupling ratio in both case (figure [[fig:coupling_ration_light_heavy]]).
 We obtain the same result than the analytical case (figures [[fig:coupling_light]] and [[fig:coupling_heavy]]).
 #+begin_src matlab :results silent :exports none
   figure;
@@ -892,18 +893,17 @@ We obtain the same result than the analytical case (figures [[fig:coupling_light
 #+end_src
 
 #+HEADER: :tangle no :exports results :results file :noweb yes
-#+HEADER: :var filepath="Figures/coupling_ration_light_heavy.png" :var figsize="wide-tall"
+#+HEADER: :var filepath="Figures/coupling_ratio_light_heavy.png" :var figsize="wide-tall"
 #+begin_src matlab
   <<plt-matlab>>
 #+end_src
 
-#+NAME: fig:coupling_ration_light_heavy
+#+NAME: fig:coupling_ratio_light_heavy
 #+RESULTS:
-[[file:Figures/coupling_ration_light_heavy.png]]
+[[file:Figures/coupling_ratio_light_heavy.png]]
 
 ** Plant Control
-
-The goal is the study control problems due to the coupling that appears because of the rotation.
+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
   %% Options for Linearized
@@ -923,44 +923,210 @@ The goal is the study control problems due to the coupling that appears because
 
 First, we identify the system when the rotation speed is null and then when the rotation speed is equal to 60rpm.
 
-The actuators are voice coil with some damping.
+The actuators are voice coil with some damping added.
+
+The bode plot of the system not rotating and rotating at 60rpm is shown figure [[fig:Gvc_speed]].
 
 #+begin_src matlab :exports none :results silent
   w = 0; % Rotation speed [rad/s]
   m = mlight; % mass of the sample [kg]
   kTuv = kvc;
-  % cTuv = 0.1*sqrt(kTuv*m);
-  cTuv = 0;
+  cTuv = 0.1*sqrt(kTuv*m);
 
-  G = linearize(mdl, io, 0.1);
-  G.InputName  = {'Fu', 'Fv'};
-  G.OutputName = {'Du', 'Dv'};
+  Gvc = linearize(mdl, io, 0.1);
+  Gvc.InputName  = {'Fu', 'Fv'};
+  Gvc.OutputName = {'Du', 'Dv'};
 #+end_src
 
 #+begin_src matlab :exports none :results silent
-  w = 0.1; % Rotation speed [rad/s]
+  w = wlight; % Rotation speed [rad/s]
   m = mlight; % mass of the sample [kg]
   kTuv = kvc;
-  % cTuv = 0.1*sqrt(kTuv*m);
-  cTuv = 0;
+  cTuv = 0.1*sqrt(kTuv*m);
 
-  Gt = linearize(mdl, io, 0.1);
-  Gt.InputName  = {'Fu', 'Fv'};
-  Gt.OutputName = {'Du', 'Dv'};
+  Gtvc = linearize(mdl, io, 0.1);
+  Gtvc.InputName  = {'Fu', 'Fv'};
+  Gtvc.OutputName = {'Du', 'Dv'};
 #+end_src
 
 #+begin_src matlab :exports none :results silent
   figure;
-  bode(G, 'r-', Gt, 'b--')
-  legend({'G - w = 0', 'G - w = 60rpm'}, 'Location', 'southwest');
+  bode(Gvc, Gtvc)
+  legend({'Gvc - $\omega = 0$', 'Gvc - $\omega = 60$rpm'}, 'Location', 'southwest');
 #+end_src
 
 #+HEADER: :tangle no :exports results :results file :noweb yes
-#+HEADER: :var filepath="Figures/coupling_ration_light_heavy.png" :var figsize="wide-tall"
+#+HEADER: :var filepath="Figures/Gvc_speed.png" :var figsize="full-tall"
 #+begin_src matlab
   <<plt-matlab>>
 #+end_src
 
+#+NAME: fig:Gvc_speed
+#+CAPTION: Bode plot of the system not rotating and rotating at 60rmp - Voice coil and light sample
+#+RESULTS:
+[[file:Figures/Gvc_speed.png]]
+
+*** Controller design
+We design a controller based on the identification when the system is not rotating.
+
+#+begin_src matlab :results none :exports code
+  sisotool(Gvc('Du', 'fu'))
+#+end_src
+
+The controller is a lead-lag controller with the following transfer function.
+#+begin_src matlab :results none :exports code
+  Kll = 2.0698e09*(s+40.45)*(s+1.181)/(s*(s+198.4)*(s+2790));
+  K = [Kll 0;
+       0   Kll];
+#+end_src
+
+The loop gain is displayed figure [[fig:Gvc_loop_gain]].
+
+#+begin_src matlab :exports none :results silent
+  freqs = logspace(-2, 2, 1000);
+
+  figure;
+  % Amplitude
+  ax1 = subaxis(2,1,1);
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(Gvc('Du', 'fu')*Kll, freqs, 'Hz'))), '-');
+  set(gca,'xscale','log'); set(gca,'yscale','log');
+  ylabel('Amplitude [m/N]');
+  set(gca, 'XTickLabel',[]);
+  hold off;
+  % Phase
+  ax2 = subaxis(2,1,2);
+  hold on;
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gvc('Du', 'fu')*Kll, freqs, 'Hz')))), '-');
+  set(gca,'xscale','log');
+  yticks(-180:180:180);
+  ylim([-180 180]);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  hold off;
+  linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/Gvc_loop_gain.png" :var figsize="full-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:Gvc_loop_gain
+#+CAPTION: Loop gain obtained for a lead-lag controller on the system with a voice coil
+#+RESULTS:
+[[file:Figures/Gvc_loop_gain.png]]
+
+*** Controlling the rotating system
+We here want to see if the system is robust with respect to the rotation speed. We then use the controller based on the non-rotating system, and see if the system is stable and its dynamics.
+
+We can then plot the same loop gain with the rotating system using the same controller (figure [[fig:Gtvc_loop_gain]]). The result obtained is unstable.
+
+#+begin_src matlab :exports none :results silent
+  freqs = logspace(-2, 2, 1000);
+
+  figure;
+  % Amplitude
+  ax1 = subaxis(2,1,1);
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(Gtvc('Du', 'fu')*Kll, freqs, 'Hz'))), '-');
+  set(gca,'xscale','log'); set(gca,'yscale','log');
+  ylabel('Amplitude [m/N]');
+  set(gca, 'XTickLabel',[]);
+  hold off;
+  % Phase
+  ax2 = subaxis(2,1,2);
+  hold on;
+  plot(freqs, 180/pi*angle(squeeze(freqresp(Gtvc('Du', 'fu')*Kll, freqs, 'Hz'))), '-');
+  set(gca,'xscale','log');
+  yticks(-180:180:180);
+  ylim([-180 180]);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  hold off;
+  linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/Gtvc_loop_gain.png" :var figsize="full-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+RESULTS:
+[[file:Figures/Gtvc_loop_gain.png]]
+
+We can look at the poles of the system where we control only one direction ($u$ for instance). We obtain a pole with a positive real part.
+#+begin_src matlab :results table :exports both
+  pole(feedback(Gtvc, blkdiag(Kll, 0)))
+#+end_src
+
+#+RESULTS:
+|           -2798 |
+| -58.906+94.246i |
+| -58.906-94.246i |
+|         -71.654 |
+|          3.1648 |
+|         -3.3031 |
+|         -1.1902 |
+
+However, when we look at the poles of the closed loop with a diagonal controller, all the poles have negative real part and the system is stable.
+#+begin_src matlab :results table :exports both
+  pole(feedback(Gtvc, blkdiag(Kll, Kll)))
+#+end_src
+
+#+RESULTS:
+| -2798+0.035765i    |
+| -2798-0.035765i    |
+| -56.406+105.34i    |
+| -56.406-105.34i    |
+| -64.482+79.308i    |
+| -64.482-79.308i    |
+| -68.521+13.503i    |
+| -68.521-13.503i    |
+| -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
+
+#+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
+
+
+** test
 
 #+begin_src matlab :exports none :results silent
   figure;
@@ -1004,10 +1170,6 @@ The actuators are voice coil with some damping.
   linkaxes([ax1,ax2],'x');
 #+end_src
 
-
-
-
-
 #+begin_src matlab :exports none :results silent
   figure;
   bode(Gpz_light, Gvc_light);
@@ -1072,13 +1234,13 @@ Plot the ratio between the main transfer function and the coupling term:
 #+end_src
 
 #+HEADER: :tangle no :exports results :results file :noweb yes
-#+HEADER: :var filepath="Figures/coupling_ration_simscape_light.png" :var figsize="wide-tall"
+#+HEADER: :var filepath="Figures/coupling_ratio_simscape_light.png" :var figsize="wide-tall"
 #+begin_src matlab
   <<plt-matlab>>
 #+end_src
 
 #+RESULTS:
-[[file:Figures/coupling_ration_simscape_light.png]]
+[[file:Figures/coupling_ratio_simscape_light.png]]
 
 #+begin_src matlab :results silent :exports none
   freqs = logspace(-2, 3, 1000);
@@ -1095,13 +1257,13 @@ Plot the ratio between the main transfer function and the coupling term:
 #+end_src
 
 #+HEADER: :tangle no :exports results :results file :noweb yes
-#+HEADER: :var filepath="Figures/coupling_ration_simscape_heavy.png" :var figsize="wide-tall"
+#+HEADER: :var filepath="Figures/coupling_ratio_simscape_heavy.png" :var figsize="wide-tall"
 #+begin_src matlab
   <<plt-matlab>>
 #+end_src
 
 #+RESULTS:
-[[file:Figures/coupling_ration_simscape_heavy.png]]
+[[file:Figures/coupling_ratio_simscape_heavy.png]]
 
 *** Low rotation speed and High rotation speed
 #+begin_src matlab :exports code :results silent