Update the analysis

css/readtheorg.css

@@ -63,7 +63,7 @@ body{
     height:100%;
     margin-left:300px;
     /* margin:auto; */
-    max-width:800px;
+    max-width:1200px;
     min-height:100%;
     padding:1.618em 3.236em;
 }

index.org

@@ -1,3 +1,3 @@
 #+TITLE: Rotating Frame
 
-- [[file:rotating_frame.html][Rotating Frame]]
+- [[file:rotating_frame.org][Rotating Frame]]

rotating_frame.org

@@ -13,29 +13,38 @@
 #+LaTeX_HEADER: \newcommand{\authorLastName}{Dehaeze}
 #+LaTeX_HEADER: \newcommand{\authorEmail}{dehaeze.thomas@gmail.com}
 #+PROPERTY: header-args:matlab  :session *MATLAB*
-#+PROPERTY: header-args:matlab+ :tangle rotating_frame.m
 #+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
 
-* Goal of this note
-The control objective is to stabilize the position of a rotating object with respect to a non-rotating frame.
+* 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.
 
-The actuators are also rotating with the object.
+In section [[sec:system]], a simple system composed of a spindle and a translation stage is defined and the equations of motion are written.
+The rotation induces some coupling between the actuators and their displacement, and modifies the dynamics of the system.
+This is studied using the equations, and some numerical computations are used to compare the use of voice coil and piezoelectric actuators.
 
-We want to compare the two different approach:
+Then, in section [[sec:control_strategies]], two different control approach are compared where:
 - the measurement is made in the fixed frame
 - the measurement is made in the rotating frame
 
+In section [[sec:simscape]], the analytical study will be validated using a multi body model of the studied system.
+
+Finally, in section [[sec:control]], the control strategies are implemented using Simulink and Simscape and compared.
+
 * System
+  :PROPERTIES:
+  :HEADER-ARGS:matlab+: :tangle system_numerical_analysis.m
+  :END:
   <<sec:system>>
 ** System description
 The system consists of one 2 degree of freedom translation stage on top of a spindle (figure [[fig:rotating_frame_2dof]]).
 
-The control inputs are the forces applied in the translation stage ($F_u$ and $F_v$). As the translation stage is rotating around the Z axis due to the spindle, the forces are applied along $u$ and $v$.
+The control inputs are the forces applied by the actuators of the translation stage ($F_u$ and $F_v$).
+As the translation stage is rotating around the Z axis due to the spindle, the forces are applied along $u$ and $v$.
 
-The measurement is either the $x-y$ displacement of the object located on top of the translation stage or the $u-v$ displacement of the actuators.
+The measurement is either the $x-y$ displacement of the object located on top of the translation stage or the $u-v$ displacement of the sample with respect to a fixed reference frame.
 
 #+name: fig:rotating_frame_2dof
 #+caption: Schematic of the mecanical system
@@ -57,7 +66,6 @@ Indices $u$ and $v$ corresponds to signals in the rotating reference frame ($\ve
 
 ** Equations
   <<sec:equations>>
-
 Based on the figure [[fig:rotating_frame_2dof]], we can write the equations of motion of the system.
 
 Let's express the kinetic energy $T$ and the potential energy $V$ of the mass $m$:
@@ -131,97 +139,271 @@ Which is equivalent to:
 
 We can then subtract and add the previous equations to obtain the following equations:
 #+begin_important
-\begin{align*}
- m \ddot{d_u} + (k - m\dot{\theta}^2) d_u &= F_u + 2 m\dot{d_v}\dot{\theta} + m d_v\ddot{\theta} \\
- m \ddot{d_v} + (k - m\dot{\theta}^2) d_v &= F_v - 2 m\dot{d_u}\dot{\theta} - m d_u\ddot{\theta} \\
-\end{align*}
+#+NAME: eq:du_coupled
+\begin{equation}
+ m \ddot{d_u} + (k - m\dot{\theta}^2) d_u = F_u + 2 m\dot{d_v}\dot{\theta} + m d_v\ddot{\theta}
+\end{equation}
+#+NAME: eq:dv_coupled
+\begin{equation}
+ m \ddot{d_v} + (k - m\dot{\theta}^2) d_v = F_v - 2 m\dot{d_u}\dot{\theta} - m d_u\ddot{\theta}
+\end{equation}
 #+end_important
 
-** TODO Analysis
 We obtain two differential equations that are coupled through:
 - *Euler forces*: $m d_v \ddot{\theta}$
 - *Coriolis forces*: $2 m \dot{d_v} \dot{\theta}$
 
-Without the coupling terms, each equation is the equation of a one degree of freedom mass-spring system with mass $m$ and stiffness $k-d_u m\dot{\theta}^2$.
-Thus, the term $-d_u m\dot{\theta}^2$ acts like a negative stiffness (due to *centrifugal forces*).
+Without the coupling terms, each equation is the equation of a one degree of freedom mass-spring system with mass $m$ and stiffness $k- m\dot{\theta}^2$.
+Thus, the term $- m\dot{\theta}^2$ acts like a negative stiffness (due to *centrifugal forces*).
 
-*** Stiff actuators
-Let's say we use stiff actuators such that $m \ddot{d_u} + (k - m\dot{\theta}^2) d_u \approx k d_u$.
+The forces induced by the rotating reference frame are independent of the stiffness of the actuator.
+The resulting effect of those forces should then be higher when using softer actuators.
 
-Let's suppose that $F_u + 2 m\dot{d_v}\dot{\theta} + m d_v\ddot{\theta} \approx F_u$.
-
-Then we obtain $d_u = \frac{F_u}{k}$ that we can re inject in the other equation to obtain:
-\[ m \ddot{d_v} + (k - m\dot{\theta}^2) d_v &= F_v - 2 m\frac{\dot{F_u}}{k}\dot{\theta} - m \frac{F_u}{k}\ddot{\theta} \]
-
-*** Negative Stiffness
-If $\max{\dot{\theta}} \ll \sqrt{\frac{k}{m}}$, then the negative spring effect is negligible and $k - m\dot{\theta}^2 \approx k$.
-
-* Analytical Computation of forces for the NASS
-For the NASS, the Euler forces should be less of a problem as $\ddot{\theta}$ should be very small when conducting an experiment.
-
-** Parameters
+** Numerical Values for the NASS
 #+begin_src matlab :exports none :results silent :noweb yes
   <<matlab-init>>
 #+end_src
 
 Let's define the parameters for the NASS.
-#+begin_src matlab :exports code :results silent
-  mlight = 35; % [kg]
-  mheavy = 85; % [kg]
+#+begin_src matlab :exports none :results silent
+  mlight = 35; % Mass for light sample [kg]
+  mheavy = 85; % Mass for heavy sample [kg]
 
-  wlight = 2*pi; % [rad/s]
-  wheavy = 2*pi/60; % [rad/s]
+  wlight = 2*pi; % Max rot. speed for light sample [rad/s]
+  wheavy = 2*pi/60; % Max rot. speed for heavy sample [rad/s]
 
-  wdot = 1; % [rad/s2]
+  kvc = 1e3; % Voice Coil Stiffness [N/m]
+  kpz = 1e8; % Piezo Stiffness [N/m]
 
-  d = 0.1; % [m]
-  ddot = 0.2; % [m/s]
+  wdot = 1; % Maximum rotation acceleration [rad/s2]
+
+  d = 0.01; % Maximum excentricity from rotational axis [m]
+  ddot = 0.2; % Maximum Horizontal speed [m/s]
+
+  save('./mat/parameters.mat');
 #+end_src
 
-** Euler and Coriolis forces
+#+begin_src matlab :results table :exports results
+  labels = {'Light sample mass [kg]', ...
+            'Heavy sample mass [kg]', ...
+            'Max rot. speed - light [rpm]', ...
+            'Max rot. speed - heavy [rpm]', ...
+            'Voice Coil Stiffness [N/m]', ...
+            'Piezo Stiffness [N/m]', ...
+            'Max rot. acceleration [rad/s2]', ...
+            'Max mass excentricity [m]', ...
+            'Max Horizontal speed [m/s]'};
+  data = [mlight, mheavy, 60*wlight/2/pi, 60*wheavy/2/pi, kvc, kpz, wdot, d, ddot];
+  data2orgtable(data', labels, {}, ' %.1e ')
+#+end_src
+
+#+RESULTS:
+| Light sample mass [kg]         | 3.5e+01 |
+| Heavy sample mass [kg]         | 8.5e+01 |
+| Max rot. speed - light [rpm]   | 6.0e+01 |
+| Max rot. speed - heavy [rpm]   | 1.0e+00 |
+| Voice Coil Stiffness [N/m]     | 1.0e+03 |
+| Piezo Stiffness [N/m]          | 1.0e+08 |
+| Max rot. acceleration [rad/s2] | 1.0e+00 |
+| Max mass excentricity [m]      | 1.0e-02 |
+| Max Horizontal speed [m/s]     | 2.0e-01 |
+
+** Euler and Coriolis forces - Numerical Result
 First we will determine the value for Euler and Coriolis forces during regular experiment.
+- *Euler forces*: $m d_v \ddot{\theta}$
+- *Coriolis forces*: $2 m \dot{d_v} \dot{\theta}$
 
 #+begin_src matlab :exports none :results silent
   Felight = mlight*d*wdot;
   Feheavy = mheavy*d*wdot;
 
-  Fclight = 2*mlight*d*wlight;
-  Fcheavy = 2*mheavy*d*wheavy;
+  Fclight = 2*mlight*ddot*wlight;
+  Fcheavy = 2*mheavy*ddot*wheavy;
 #+end_src
 
-We then compute the corresponding values of the Coriolis and Euler forces, and the obtained values are displayed in table [[tab:euler_coriolis]].
+The obtained values are displayed in table [[tab:euler_coriolis]].
 
 #+begin_src matlab :results value table :exports results :post addhdr(*this*)
-  ans = sprintf(' | Light | Heavy | \n Coriolis | %.1f N | %.1f N | \n Euler | %.1f N | %.1f N', Fclight, Fcheavy, Felight, Feheavy)
+  data = [Fclight, Fcheavy ;
+          Felight, Feheavy];
+  data2orgtable(data, {'Coriolis', 'Euler'}, {'Light', 'Heavy'}, ' %.1fN ')
 #+end_src
 
 #+NAME: tab:euler_coriolis
 #+CAPTION: Euler and Coriolis forces for the NASS
 #+RESULTS:
-|          | Light  | Heavy |
-|----------+--------+-------|
-| Coriolis | 44.0 N | 1.8 N |
-| Euler    | 3.5 N  | 8.5 N |
+|          | Light | Heavy |
+|----------+-------+-------|
+| Coriolis | 88.0N | 3.6N  |
+| Euler    | 0.4N  | 0.8N  |
+
+** Negative Spring Effect - Numerical Result
+The negative stiffness due to the rotation is equal to $-m{\omega_0}^2$.
 
-** Negative Spring Effect
 #+begin_src matlab :exports none :results silent
-  Klight = mlight*d*wdot^2;
-  Kheavy = mheavy*d*wdot^2;
+  Klight = mlight*wlight^2;
+  Kheavy = mheavy*wheavy^2;
 #+end_src
 
 The values for the negative spring effect are displayed in table [[tab:negative_spring]].
+
 This is definitely negligible when using piezoelectric actuators. It may not be the case when using voice coil actuators.
 
 #+begin_src matlab :results value table :exports results :post addhdr(*this*)
-  ans = sprintf(' | Light | Heavy | \n Neg. Spring | %.1f N/m | %.1f N/m', Klight, Kheavy)
+  data = [Klight, Kheavy];
+  data2orgtable(data, {'Neg. Spring'}, {'Light', 'Heavy'}, ' %.1f[N/m] ')
 #+end_src
 
 #+NAME: tab:negative_spring
 #+CAPTION: Negative Spring effect
 #+RESULTS:
-|             | Light   | Heavy   |
-|-------------+---------+---------|
-| Neg. Spring | 3.5 N/m | 8.5 N/m |
+|             | Light       | Heavy    |
+|-------------+-------------+----------|
+| Neg. Spring | 1381.7[N/m] | 0.9[N/m] |
+
+** Limitations due to coupling
+To simplify, we consider a constant rotating speed $\dot{\theta} = {\omega_0}$ and thus $\ddot{\theta} = 0$.
+
+From equations [[eq:du_coupled]] and [[eq:dv_coupled]], we obtain:
+ \begin{align*}
+ (m s^2 + (k - m{\omega_0}^2)) d_u &= F_u + 2 m {\omega_0} s d_v \\
+ (m s^2 + (k - m{\omega_0}^2)) d_v &= F_v - 2 m {\omega_0} s d_u \\
+\end{align*}
+
+From second equation:
+\[ d_v = \frac{1}{m s^2 + (k - m{\omega_0}^2)} F_v - \frac{2 m {\omega_0} s}{m s^2 + (k - m{\omega_0}^2)} d_u \]
+
+And we re-inject $d_v$ into the first equation:
+\begin{equation*}
+(m s^2 + (k - m{\omega_0}^2)) d_u = F_u + \frac{2 m {\omega_0} s}{m s^2 + (k - m{\omega_0}^2)} F_v - \frac{(2 m {\omega_0} s)^2}{m s^2 + (k - m{\omega_0}^2)} d_u
+\end{equation*}
+
+\begin{equation*}
+\frac{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2}{m s^2 + (k - m{\omega_0}^2)} d_u = F_u + \frac{2 m {\omega_0} s}{m s^2 + (k - m{\omega_0}^2)} F_v
+\end{equation*}
+
+Finally we obtain $d_u$ function of $F_u$ and $F_v$.
+\[ d_u = \frac{m s^2 + (k - m{\omega_0}^2)}{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2} F_u + \frac{2 m {\omega_0} s}{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2} F_v \]
+
+Similarly we can obtain $d_v$ function of $F_u$ and $F_v$:
+\[ d_v = \frac{m s^2 + (k - m{\omega_0}^2)}{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2} F_v - \frac{2 m {\omega_0} s}{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2} F_u \]
+
+The two previous equations can be written in a matrix form:
+#+begin_important
+\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}
+\begin{bmatrix}
+  ms^2 + (k-m{\omega_0}^2) & 2 m \omega_0 s \\
+  -2 m \omega_0 s          & ms^2 + (k-m{\omega_0}^2) \\
+\end{bmatrix}
+\begin{bmatrix} F_u \\ F_v \end{bmatrix}
+\end{equation}
+#+end_important
+
+Then, coupling is negligible if $|-m \omega^2 + (k - m{\omega_0}^2)| \gg |2 m {\omega_0} \omega|$.
+
+*** Numerical Analysis
+We plot on the same graph $\frac{|-m \omega^2 + (k - m {\omega_0}^2)|}{|2 m \omega_0 \omega|}$ for the voice coil and the piezo:
+- with the light sample (figure [[fig:coupling_light]]).
+- with the heavy sample (figure [[fig:coupling_heavy]]).
+
+#+HEADER: :exports none :results silent
+#+begin_src matlab
+  f = logspace(-1, 2, 1000);
+
+  figure;
+  hold on;
+  plot(f, abs(-mlight*(2*pi*f).^2 + kvc - mlight * wlight^2)./abs(2*mlight*wlight*2*pi*f), 'DisplayName', 'Voice Coil')
+  plot(f, abs(-mlight*(2*pi*f).^2 + kpz - mlight * wlight^2)./abs(2*mlight*wlight*2*pi*f), 'DisplayName', 'Piezo')
+  plot(f, ones(1, length(f)), 'k--', 'HandleVisibility', 'off')
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]');
+  legend('Location', 'northeast');
+  hold off;
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/coupling_light.png" :var figsize="normal-normal"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:coupling_light
+#+CAPTION: Relative Coupling for light mass and high rotation speed
+#+RESULTS:
+[[file:Figures/coupling_light.png]]
+
+#+HEADER: :exports none :results silent
+#+begin_src matlab
+  f = logspace(-1, 2, 1000);
+
+  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')
+  plot(f, abs(-mheavy*(2*pi*f).^2 + kpz - mheavy * wheavy^2)./abs(2*mheavy*wheavy*2*pi*f), 'DisplayName', 'Piezo')
+  plot(f, ones(1, length(f)), 'k--', 'HandleVisibility', 'off')
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]');
+  legend('Location', 'northeast');
+  hold off;
+#+end_src
+
+#+HEADER: :tangle no :exports results :results file :noweb yes
+#+HEADER: :var filepath="Figures/coupling_heavy.png" :var figsize="normal-normal"
+#+begin_src matlab
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:coupling_heavy
+#+CAPTION: Relative Coupling for heavy mass and low rotation speed
+#+RESULTS:
+[[file:Figures/coupling_heavy.png]]
+
+#+begin_important
+Coupling is higher for actuators with small stiffness.
+#+end_important
+
+** Limitations due to negative stiffness effect
+If $\max{\dot{\theta}} \ll \sqrt{\frac{k}{m}}$, then the negative spring effect is negligible and $k - m\dot{\theta}^2 \approx k$.
+
+Let's estimate what is the maximum rotation speed for which the negative stiffness effect is still negligible ($\omega_\text{max} = 0.1 \sqrt{\frac{k}{m}}$). Results are shown table [[tab:negative_stiffness]].
+#+begin_src matlab :results table :exports results :post addhdr(*this*)
+  data = 0.1*60*(1/2/pi)*[sqrt(kvc/mlight), sqrt(kpz/mlight);
+                       sqrt(kvc/mheavy), sqrt(kpz/mheavy)];
+  data2orgtable(data, {'Light', 'Heavy'}, {'Voice Coil', 'Piezo'}, ' %.0f[rpm] ')
+#+end_src
+
+#+NAME: tab:negative_stiffness
+#+CAPTION: Maximum rotation speed at which negative stiffness is negligible ($0.1\sqrt{\frac{k}{m}}$)
+#+RESULTS:
+|       | Voice Coil | Piezo     |
+|-------+------------+-----------|
+| Light | 5[rpm]     | 1614[rpm] |
+| Heavy | 3[rpm]     | 1036[rpm] |
+
+The negative spring effect is proportional to the rotational speed $\omega$.
+The system dynamics will be much more affected when using soft actuator.
+
+#+begin_important
+  Negative stiffness effect has very important effect when using soft actuators.
+#+end_important
+
+The system can even goes unstable when $m \omega^2 > k$, that is when the centrifugal forces are higher than the forces due to stiffness.
+
+From this analysis, we can determine the lowest practical stiffness that is possible to use: $k_\text{min} = 10 m \omega^2$ (table [[tab:min_k]])
+
+#+begin_src matlab :results table :exports results :post addhdr(*this*)
+  data = 10*[mlight*2*pi, mheavy*2*pi/60]
+  data2orgtable(data, {'k min [N/m]'}, {'Light', 'Heavy'}, ' %.0f ')
+#+end_src
+
+#+NAME: tab:min_k
+#+CAPTION: Minimum possible stiffness
+#+RESULTS:
+|             | Light | Heavy |
+|-------------+-------+-------|
+| k min [N/m] |  2199 |    89 |
 
 * Control Strategies
   <<sec:control_strategies>>
@@ -257,50 +439,264 @@ The corresponding block diagram is shown figure [[fig:control_measure_rotating_2
 
 The loop gain is $L = G K$.
 
-* Effect of the rotating Speed
+* Multi Body Model - Simscape
+  :PROPERTIES:
+  :HEADER-ARGS:matlab+: :tangle simscape_analysis.m
+  :END:
+  <<sec:simscape>>
+
+#+begin_src matlab :exports none :results silent :noweb yes
+  <<matlab-init>>
+  load('./mat/parameters.mat');
+#+end_src
+
+#+begin_src matlab :exports none :results silent
+  open rotating_frame.slx
+#+end_src
+
+** Identification in the rotating referenced frame
+We initialize the inputs and outputs of the system to identify.
+#+begin_src matlab :exports code :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
+
+*** Piezo and Voice coil
+We start we identify the transfer functions at high speed with the light sample.
+#+begin_src matlab :exports code :results silent
+  rot_speed = wlight;
+  angle_e = 0;
+  m = mlight;
+
+  k = kpz;
+  c = 1e3;
+  Gpz_light = linearize(mdl, io, 0.1);
+
+  k = kvc;
+  c = 1e3;
+  Gvc_light = linearize(mdl, io, 0.1);
+
+  Gpz_light.InputName  = {'Fu', 'Fv'};
+  Gpz_light.OutputName = {'Du', 'Dv'};
+  Gvc_light.InputName  = {'Fu', 'Fv'};
+  Gvc_light.OutputName = {'Du', 'Dv'};
+#+end_src
+
+#+begin_src matlab :exports none :results silent
+  figure;
+  bode(Gpz_light, Gvc_light);
+#+end_src
+
+And then with the heavy sample.
+#+begin_src matlab :exports code :results silent
+  rot_speed = wheavy;
+  angle_e = 0;
+  m = mheavy;
+
+  k = kpz;
+  c = 1e3;
+  Gpz_heavy = linearize(mdl, io, 0.1);
+
+  k = kvc;
+  c = 1e3;
+  Gvc_heavy = linearize(mdl, io, 0.1);
+
+  Gpz_heavy.InputName  = {'Fu', 'Fv'};
+  Gpz_heavy.OutputName = {'Du', 'Dv'};
+  Gvc_heavy.InputName  = {'Fu', 'Fv'};
+  Gvc_heavy.OutputName = {'Du', 'Dv'};
+#+end_src
+
+#+begin_src matlab :exports none :results silent
+  figure;
+  bode(Gpz_heavy, Gvc_heavy);
+#+end_src
+
+Plot the ratio between the main transfer function and the coupling term:
+#+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')));
+  plot(freqs, abs(squeeze(freqresp(Gpz_light('Du', 'Fu'), freqs, 'Hz'))))./abs(squeeze(freqresp(Gpz_light('Dv', 'Fu'), freqs, 'Hz')));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Coupling ratio');
+  legend({'Voice Coil', 'Piezoelectric'})
+#+end_src
+
+#+begin_src matlab :results silent :exports none
+  freqs = logspace(-2, 3, 1000);
+
+  figure;
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(Gvc_heavy('Du', 'Fu'), freqs, 'Hz'))))./abs(squeeze(freqresp(Gvc_heavy('Dv', 'Fu'), freqs, 'Hz')));
+  plot(freqs, abs(squeeze(freqresp(Gpz_heavy('Du', 'Fu'), freqs, 'Hz'))))./abs(squeeze(freqresp(Gpz_heavy('Dv', 'Fu'), freqs, 'Hz')));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Coupling ratio');
+  legend({'Voice Coil', 'Piezoelectric'})
+#+end_src
+
+*** Low rotation speed and High rotation speed
+#+begin_src matlab :exports code :results silent
+  rot_speed = 2*pi/60; angle_e = 0;
+  G_low = linearize(mdl, io, 0.1);
+
+  rot_speed = 2*pi; angle_e = 0;
+  G_high = linearize(mdl, io, 0.1);
+
+  G_low.InputName  = {'Fu', 'Fv'};
+  G_low.OutputName = {'Du', 'Dv'};
+  G_high.InputName  = {'Fu', 'Fv'};
+  G_high.OutputName = {'Du', 'Dv'};
+#+end_src
+
+#+begin_src matlab :results silent
+  figure;
+  bode(G_low, G_high);
+#+end_src
+
+** Identification in the fixed frame
+Let's define some options as well as the inputs and outputs for linearization.
+#+begin_src matlab :exports code :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, '/fx'], 1, 'input');
+  io(2) = linio([mdl, '/fy'], 1, 'input');
+
+  io(3) = linio([mdl, '/dx'], 1, 'output');
+  io(4) = linio([mdl, '/dy'], 1, 'output');
+#+end_src
+
+We then define the error estimation of the error and the rotational speed.
+#+begin_src matlab :exports code :results silent
+  %% Run the linearization
+  angle_e = 0;
+  rot_speed = 0;
+#+end_src
+
+Finally, we run the linearization.
+#+begin_src matlab :exports code :results silent
+  G = linearize(mdl, io, 0);
+
+  %% Input/Output names
+  G.InputName  = {'Fx', 'Fy'};
+  G.OutputName = {'Dx', 'Dy'};
+#+end_src
+
+#+begin_src matlab :exports code :results silent
+  %% Run the linearization
+  angle_e = 0;
+  rot_speed = 2*pi;
+  Gr = linearize(mdl, io, 0);
+
+  %% Input/Output names
+  Gr.InputName  = {'Fx', 'Fy'};
+  Gr.OutputName = {'Dx', 'Dy'};
+#+end_src
+
+#+begin_src matlab :exports code :results silent
+  %% Run the linearization
+  angle_e = 1*2*pi/180;
+  rot_speed = 2*pi;
+  Ge = linearize(mdl, io, 0);
+
+  %% Input/Output names
+  Ge.InputName  = {'Fx', 'Fy'};
+  Ge.OutputName = {'Dx', 'Dy'};
+#+end_src
+
+#+begin_src matlab :exports code :results silent
+  figure;
+  bode(G);
+  % exportFig('G_x_y', 'wide-tall');
+
+  figure;
+  bode(Ge);
+  % exportFig('G_x_y_e', 'normal-normal');
+#+end_src
+
+** Identification from actuator forces to displacement in the fixed frame
+#+begin_src matlab :exports code :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, '/dx'], 1, 'output');
+  io(4) = linio([mdl, '/dy'], 1, 'output');
+#+end_src
+
+#+begin_src matlab :exports code :results silent
+  rot_speed = 2*pi;
+  angle_e = 0;
+  G = linearize(mdl, io, 0.0);
+
+  G.InputName  = {'Fu', 'Fv'};
+  G.OutputName = {'Dx', 'Dy'};
+#+end_src
+
+#+begin_src matlab :exports code :results silent
+  rot_speed = 2*pi;
+  angle_e = 0;
+  G1 = linearize(mdl, io, 0.4);
+
+  G1.InputName  = {'Fu', 'Fv'};
+  G1.OutputName = {'Dx', 'Dy'};
+#+end_src
+
+#+begin_src matlab :exports code :results silent
+  rot_speed = 2*pi;
+  angle_e = 0;
+  G2 = linearize(mdl, io, 0.8);
+
+  G2.InputName  = {'Fu', 'Fv'};
+  G2.OutputName = {'Dx', 'Dy'};
+#+end_src
+
+#+begin_src matlab :exports code :results silent
+  figure;
+  bode(G, G1, G2);
+  exportFig('G_u_v_to_x_y', 'wide-tall');
+#+end_src
+
+** Effect of the rotating Speed
   <<sec:effect_rot_speed>>
 
 #+begin_src matlab :exports none :results silent :noweb yes
   <<matlab-init>>
 #+end_src
 
-** TODO Use realistic parameters for the mass of the sample and stiffness of the X-Y stage
-
-** TODO Check if the plant is changing a lot when we are not turning to when we are turning at the maximum speed (60rpm)
-
-* Effect of the X-Y stage stiffness
+*** TODO Use realistic parameters for the mass of the sample and stiffness of the X-Y stage
+*** TODO Check if the plant is changing a lot when we are not turning to when we are turning at the maximum speed (60rpm)
+** Effect of the X-Y stage stiffness
   <<sec:effect_stiffness>>
-** TODO At full speed, check how the coupling changes with the stiffness of the actuators
-* Noweb Snippets                                                   :noexport:
-  :PROPERTIES:
-  :HEADER-ARGS:matlab+: :results none :exports none :tangle no
-  :HEADER-ARGS:emacs-lisp+: :tangle no :eval no-export
-  :END:
-
-** Matlab Init
-#+NAME: matlab-init
-#+BEGIN_SRC matlab :results none :exports none
-  clear; close all; clc;
-
-  %% Add path with some functions
-  addpath('./src/');
-
-  %% Intialize Laplace variable
-  s = tf('s');
-
-  %% Initialize ans with org-babel
-  ans = 0;
-#+END_SRC
-
-** Matlab Export Figure
-#+NAME: matlab-export-figure
-#+BEGIN_SRC matlab :var fn="figure" ft="png" size="normal-normal" :exports none
-  exportFig(fn, size);
-  ans = [fn, ".", ft];
-#+END_SRC
-
-** Add hline to org table
-#+name: addhdr
-#+begin_src emacs-lisp :var tbl=""
-(cons (car tbl) (cons 'hline (cdr tbl)))
-#+end_src
+*** TODO At full speed, check how the coupling changes with the stiffness of the actuators
+* Control Implementation
+  <<sec:control>>
+** Measurement in the fixed reference frame

rotating_frame_old.org

@@ -1,54 +0,0 @@
-#+TITLE:     Control in a rotating frame
-#+AUTHOR:    Dehaeze Thomas
-#+EMAIL:     dehaeze.thomas@gmail.com
-#+DATE:      {{{time(%Y-%m-%d)}}}
-#+DESCRIPTION: Course about Optimal and Robust Control
-#+KEYWORDS:  control optimal robust lqr lqg hinfinity
-#+LANGUAGE:  en
-#+STARTUP: beamer
-#+LaTeX_CLASS: clean-beamer
-#+LaTeX_CLASS_OPTIONS: [t]
-#+OPTIONS:   H:2
-#+OPTIONS:   num:t toc:t ::t |:t ^:{} -:t f:t *:t <:t
-#+SELECT_TAGS: export
-#+EXCLUDE_TAGS: noexport
-#+COLUMNS: %20ITEM %13BEAMER_env(Env) %6BEAMER_envargs(Args) %4BEAMER_col(Col) %7BEAMER_extra(Extra)
-#+latex_header: \AtBeginSection[]{\begin{frame}<beamer>\frametitle{Topic}\tableofcontents[currentsection]\end{frame}}
-
-* Simscape Model
-#+name: fig:simscape
-#+caption: Screen of simscape multibody model
-#+attr_latex: :float t :width 0.7\linewidth
-[[./Figures/simscape.png]]
-
-* Simscape Model
-#+name: fig:simulink_ctrl
-#+caption: Simulink Blocks
-#+attr_latex: :float t :width 0.9\linewidth
-[[./Figures/simulink_ctrl.png]]
-
-* Identification in the rotating frame
-The rotation speed increase the coupling between the rotating actuators and sensors (figure ref:fig:G_u_v).
-
-#+name: fig:G_u_v
-#+caption: Transfer function from forces to displacement in the rotating frame
-#+attr_latex: :float t
-[[./Figures/G_u_v.pdf]]
-
-* Identification in the cartesian frame
-#+name: fig:G_x_y
-#+caption: Plant from force to displacement in the cartesian frame
-#+attr_latex: :float t :width 0.8\linewidth
-[[./Figures/G_x_y.pdf]]
-
-* Identification in the cartesian frame
-#+name: fig:G_x_y_e
-#+caption: Plant from force to displacement in the cartesian frame with small angle estimation error
-#+attr_latex: :float t :width 0.8\linewidth
-[[./Figures/G_x_y_e.pdf]]
-
-* Control result
-#+name: fig:D_x_y
-#+caption: Control result with a simple PID
-#+attr_latex: :float t :width 0.8\linewidth
-[[./Figures/D_x_y.pdf]]