Add cambell plot. Add plot about simscape analysis

rotating_frame.org

@@ -33,7 +33,7 @@ In section [[sec:simscape]], the analytical study will be validated using a mult
 
 Finally, in section [[sec:control]], the control strategies are implemented using Simulink and Simscape and compared.
 
-* System
+* System Description and Analysis
   :PROPERTIES:
   :HEADER-ARGS:matlab+: :tangle system_numerical_analysis.m
   :END:
@@ -290,6 +290,7 @@ Similarly we can obtain $d_v$ function of $F_u$ and $F_v$:
 
 The two previous equations can be written in a matrix form:
 #+begin_important
+#+NAME: eq:coupledplant
 \begin{equation}
 \begin{bmatrix} d_u \\ d_v \end{bmatrix} =
 \frac{1}{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2}
@@ -308,8 +309,7 @@ We plot on the same graph $\frac{|-m \omega^2 + (k - m {\omega_0}^2)|}{|2 m \ome
 - with the light sample (figure [[fig:coupling_light]]).
 - with the heavy sample (figure [[fig:coupling_heavy]]).
 
-#+HEADER: :exports none :results silent
+#+begin_src matlab :exports none :results silent
-#+begin_src matlab
   f = logspace(-1, 3, 1000);
 
   figure;
@@ -335,8 +335,7 @@ We plot on the same graph $\frac{|-m \omega^2 + (k - m {\omega_0}^2)|}{|2 m \ome
 #+RESULTS:
 [[file:Figures/coupling_light.png]]
 
-#+HEADER: :exports none :results silent
+#+begin_src matlab :exports none :results silent
-#+begin_src matlab
   figure;
   hold on;
   plot(f, abs(-mheavy*(2*pi*f).^2 + kvc - mheavy * wheavy^2)./abs(2*mheavy*wheavy*2*pi*f), 'DisplayName', 'Voice Coil')
@@ -405,6 +404,342 @@ From this analysis, we can determine the lowest practical stiffness that is poss
 |-------------+-------+-------|
 | k min [N/m] |  2199 |    89 |
 
+** Effect of rotation speed on the plant
+As shown in equation [[eq:coupledplant]], the plant changes with the rotation speed $\omega_0$.
+
+Then, we compute the bode plot of the direct term and coupling term for multiple rotating speed.
+
+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 = [0, 1, 10, 60]*2*pi/60; % [rmp]
+
+  Gs = {zeros(1, length(ws))};
+  Gcs = {zeros(1, length(ws))};
+
+  for i = 1:length(ws)
+    w = ws(i);
+    Gs(i) = {(m*s^2 + (k-m*w^2))/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)};
+    Gcs(i) = {(2*m*w*s)/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)};
+  end
+#+end_src
+
+#+begin_src matlab :exports none :results silent
+  freqs = logspace(-2, 1, 1000);
+
+  figure;
+  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));
+  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'))));
+  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]');
+  linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/G_ws_vc.png" :var figsize="wide-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:G_ws_vc
+#+CAPTION: Bode plot of the direct transfer function term (from $F_u$ to $D_u$) for multiple rotation speed - Voice coil
+#+RESULTS:
+[[file:Figures/G_ws_vc.png]]
+
+#+begin_src matlab :exports none :results silent
+  freqs = logspace(-2, 1, 1000);
+
+  figure;
+  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));
+  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'))));
+  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]');
+  linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/Gc_ws_vc.png" :var figsize="wide-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:Gc_ws_vc
+#+CAPTION: Bode plot of the coupling transfer function term (from $F_u$ to $D_v$) for multiple rotation speed - Voice coil
+#+RESULTS:
+[[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 = [0, 1, 10, 60]*2*pi/60; % [rmp]
+
+  Gs = {zeros(1, length(ws))};
+  Gcs = {zeros(1, length(ws))};
+
+  for i = 1:length(ws)
+    w = ws(i);
+    Gs(i) = {(m*s^2 + (k-m*w^2))/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)};
+    Gcs(i) = {(2*m*w*s)/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)};
+  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}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
+  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'))));
+  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]');
+  linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/G_ws_pz.png" :var figsize="wide-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:G_ws_pz
+#+CAPTION: Bode plot of the direct transfer function term (from $F_u$ to $D_u$) for multiple rotation speed - Piezoelectric actuator
+#+RESULTS:
+[[file:Figures/G_ws_pz.png]]
+
+#+begin_src matlab :exports none :results silent
+  figure;
+  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));
+  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'))));
+  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]');
+  linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/Gc_ws_pz.png" :var figsize="wide-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:Gc_ws_pz
+#+CAPTION: Bode plot of the coupling transfer function term (from $F_u$ to $D_v$) for multiple rotation speed - Piezoelectric actuator
+#+RESULTS:
+[[file:Figures/Gc_ws_pz.png]]
+
+*** Analysis
+When the rotation speed is null, the coupling terms are equal to zero and the diagonal terms corresponds to one degree of freedom mass spring system.
+
+When the rotation speed in not null, the resonance frequency is duplicated into two pairs of complex conjugate poles.
+
+As the rotation speed increases, one of the two resonant frequency goes to lower frequencies as the other one goes to higher frequencies.
+
+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.
+
+As shown in the previous figures, the system with voice coil is much more sensitive to rotation speed.
+
+*** Campbell diagram
+
+The poles of the system are computed for multiple values of the rotation frequency.
+
+#+begin_src matlab :results silent :exports code
+  m = mlight;
+  k = kvc;
+  c = 0.1*sqrt(k*m);
+
+  ws = linspace(0, 10, 100); % [rad/s]
+
+  polesvc = zeros(2, length(ws));
+
+  for i = 1:length(ws)
+    polei = pole(1/((m*s^2 + c*s + (k - m*ws(i)^2))^2 + (2*m*ws(i)*s)^2));
+    polesvc(:, i) = sort(polei(imag(polei) > 0));
+  end
+#+end_src
+
+#+begin_src matlab :results silent :exports code
+  m = mlight;
+  k = kpz;
+  c = 0.1*sqrt(k*m);
+
+  ws = linspace(0, 1000, 100); % [rad/s]
+
+  polespz = zeros(2, length(ws));
+
+  for i = 1:length(ws)
+    polei = pole(1/((m*s^2 + c*s + (k - m*ws(i)^2))^2 + (2*m*ws(i)*s)^2));
+    polespz(:, i) = sort(polei(imag(polei) > 0));
+  end
+#+end_src
+
+We then plot the real and imaginary part of the poles as a function of the rotation frequency (figures [[fig:poles_w_vc]] and [[fig:poles_w_pz]]).
+
+When the real part of one pole becomes positive, the system goes unstable.
+
+For the voice coil (figure [[fig:poles_w_vc]]), the system is unstable when the rotation speed is above 5 rad/s. The real and imaginary part of the poles of the system with piezoelectric actuators are changing much less (figure [[fig:poles_w_pz]]).
+
+#+begin_src matlab :results silent :exports none
+  figure;
+  % Amplitude
+  ax1 = subplot(1,2,1);
+  hold on;
+  plot(ws, real(polesvc(1, :)))
+  set(gca,'ColorOrderIndex',1)
+  plot(ws, real(polesvc(2, :)))
+  plot(ws, zeros(size(ws)), 'k--')
+  hold off;
+  xlabel('Rotation Frequency [rad/s]');
+  ylabel('Pole Real Part');
+  ax2 = subplot(1,2,2);
+  hold on;
+  plot(ws, imag(polesvc(1, :)))
+  set(gca,'ColorOrderIndex',1)
+  plot(ws, -imag(polesvc(1, :)))
+  set(gca,'ColorOrderIndex',1)
+  plot(ws, imag(polesvc(2, :)))
+  set(gca,'ColorOrderIndex',1)
+  plot(ws, -imag(polesvc(2, :)))
+  hold off;
+  xlabel('Rotation Frequency [rad/s]');
+  ylabel('Pole Imaginary Part');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/poles_w_vc.png" :var figsize="wide-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:poles_w_vc
+#+CAPTION: Real and Imaginary part of the poles of the system as a function of the rotation speed - Voice Coil and light sample
+#+RESULTS:
+[[file:Figures/poles_w_vc.png]]
+
+
+#+begin_src matlab :results silent :exports none
+  figure;
+  % Amplitude
+  ax1 = subplot(1,2,1);
+  hold on;
+  plot(ws, real(polespz(1, :)))
+  set(gca,'ColorOrderIndex',1)
+  plot(ws, real(polespz(2, :)))
+  plot(ws, zeros(size(ws)), 'k--')
+  hold off;
+  xlabel('Rotation Frequency [rad/s]');
+  ylabel('Pole Real Part');
+  ax2 = subplot(1,2,2);
+  hold on;
+  plot(ws, imag(polespz(1, :)))
+  set(gca,'ColorOrderIndex',1)
+  plot(ws, -imag(polespz(1, :)))
+  set(gca,'ColorOrderIndex',1)
+  plot(ws, imag(polespz(2, :)))
+  set(gca,'ColorOrderIndex',1)
+  plot(ws, -imag(polespz(2, :)))
+  hold off;
+  xlabel('Rotation Frequency [rad/s]');
+  ylabel('Pole Imaginary Part');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/poles_w_pz.png" :var figsize="wide-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:poles_w_pz
+#+CAPTION: Real and Imaginary part of the poles of the system as a function of the rotation speed - Voice Coil and light sample
+#+RESULTS:
+[[file:Figures/poles_w_pz.png]]
+
+
 * Control Strategies
   <<sec:control_strategies>>
 ** Measurement in the fixed reference frame
@@ -445,6 +780,7 @@ The loop gain is $L = G K$.
   :END:
   <<sec:simscape>>
 
+** Initialize
 #+begin_src matlab :exports none :results silent :noweb yes
   <<matlab-init>>
   load('./mat/parameters.mat');
@@ -455,6 +791,7 @@ The loop gain is $L = G K$.
 #+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]
 
@@ -467,6 +804,7 @@ The loop gain is $L = G K$.
   cTuv = 0; % Damping of the Tuv stage [N/(m/s)]
 #+end_src
 
+Then, we defined parameters that will be used in the following analysis.
 #+begin_src matlab :exports code :results silent
   mlight = 5; % Mass for light sample [kg]
   mheavy = 55; % Mass for heavy sample [kg]
@@ -478,11 +816,15 @@ The loop gain is $L = G K$.
   kpz = 1e8; % Piezo Stiffness [N/m]
 
   d = 0.01; % Maximum excentricity from rotational axis [m]
+
+  freqs = logspace(-2, 3, 1000); % Frequency vector for analysis [Hz]
 #+end_src
 
 ** Identification in the rotating referenced frame
+We initialize the inputs and outputs of the system to identify:
+- Inputs: $f_u$ and $f_v$
+- Outputs: $d_u$ and $d_v$
 
-We initialize the inputs and outputs of the system to identify.
 #+begin_src matlab :exports code :results silent
   %% Options for Linearized
   options = linearizeOptions;
@@ -531,11 +873,11 @@ Then we identify the system with an heavy mass and low speed.
   Gvc_heavy.OutputName = {'Du', 'Dv'};
 #+end_src
 
-Finally, we plot the coupling ratio in both case (figure [[fig:coupling_ration_light_heavy]]).
+** Coupling ratio between $f_{uv}$ and $d_{uv}$
+
+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
-  freqs = logspace(-2, 3, 1000);
-
   figure;
   hold on;
   plot(freqs, abs(squeeze(freqresp(Gvc_light('Du', 'Fu'), freqs, 'Hz')))./abs(squeeze(freqresp(Gvc_light('Dv', 'Fu'), freqs, 'Hz'))));
@@ -560,6 +902,108 @@ We obtain the same result than the analytical case (figures [[fig:coupling_light
 #+RESULTS:
 [[file:Figures/coupling_ration_light_heavy.png]]
 
+** Plant Control
+
+The goal is the study control problems due to the coupling that appears because of the rotation.
+
+#+begin_src matlab :exports none :results silent
+  %% Options for Linearized
+  options = linearizeOptions;
+  options.SampleTime = 0;
+
+  %% Name of the Simulink File
+  mdl = 'rotating_frame';
+
+  %% Input/Output definition
+  io(1) = linio([mdl, '/fu'], 1, 'input');
+  io(2) = linio([mdl, '/fv'], 1, 'input');
+
+  io(3) = linio([mdl, '/du'], 1, 'output');
+  io(4) = linio([mdl, '/dv'], 1, 'output');
+#+end_src
+
+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.
+
+#+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;
+
+  G = linearize(mdl, io, 0.1);
+  G.InputName  = {'Fu', 'Fv'};
+  G.OutputName = {'Du', 'Dv'};
+#+end_src
+
+#+begin_src matlab :exports none :results silent
+  w = 0.1; % Rotation speed [rad/s]
+  m = mlight; % mass of the sample [kg]
+  kTuv = kvc;
+  % cTuv = 0.1*sqrt(kTuv*m);
+  cTuv = 0;
+
+  Gt = linearize(mdl, io, 0.1);
+  Gt.InputName  = {'Fu', 'Fv'};
+  Gt.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');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/coupling_ration_light_heavy.png" :var figsize="wide-tall"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+
+#+begin_src matlab :exports none :results silent
+  figure;
+  % Amplitude
+  ax1 = subaxis(2,1,1);
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(Gvc_light('du', 'fu'), 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_light('du', 'fu'), freqs, 'Hz')))), '-');
+  set(gca,'xscale','log');
+  yticks(-180:180:180);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  hold off;
+  linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+begin_src matlab :exports none :results silent
+  figure;
+  % Amplitude
+  ax1 = subaxis(2,1,1);
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(-Gvc_light('dv', 'fv'), 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_light('dv', 'fv'), freqs, 'Hz')))), '-');
+  set(gca,'xscale','log');
+  yticks(-180:180:180);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  hold off;
+  linkaxes([ax1,ax2],'x');
+#+end_src