nass-simscape/active_damping/index.org

#+TITLE: Active Damping applied on the Simscape Model
:DRAWER:
#+STARTUP: overview

#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas

#+HTML_LINK_HOME: ../index.html
#+HTML_LINK_UP: ../index.html

#+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: <link rel="stylesheet" type="text/css" href="../css/zenburn.css"/>
#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" 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>

#+HTML_MATHJAX: align: center tagside: right font: TeX

#+PROPERTY: header-args:matlab  :session *MATLAB*
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :results none
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :output-dir figs
#+PROPERTY: header-args:matlab+ :tangle no
#+PROPERTY: header-args:matlab+ :mkdirp yes

#+PROPERTY: header-args:shell  :eval no-export

#+PROPERTY: header-args:latex  :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/Cloud/thesis/latex/}{config.tex}")
#+PROPERTY: header-args:latex+ :imagemagick t :fit yes
#+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150
#+PROPERTY: header-args:latex+ :imoutoptions -quality 100
#+PROPERTY: header-args:latex+ :results raw replace :buffer no
#+PROPERTY: header-args:latex+ :eval no-export
#+PROPERTY: header-args:latex+ :exports both
#+PROPERTY: header-args:latex+ :mkdirp yes
#+PROPERTY: header-args:latex+ :output-dir figs
:END:

* Introduction                                                       :ignore:
First, in section [[sec:undamped_system]], we will looked at the undamped system.

Then, we will compare three active damping techniques:
- In section [[sec:iff]]: the integral force feedback is used
- In section [[sec:dvf]]: the direct velocity feedback is used
- In section [[sec:ine]]: inertial control is used

For each of the active damping technique, we will:
- Look at the damped plant
- Simulate tomography experiments
- Compare the sensitivity from disturbances

The disturbances are:
- Ground motion
- Motion errors of all the stages

* Undamped System
:PROPERTIES:
:header-args:matlab+: :tangle matlab/undamped_system.m
:header-args:matlab+: :comments none :mkdirp yes
:END:
<<sec:undamped_system>>

** Introduction                                                      :ignore:
In this section, we identify the dynamic of the system from forces applied in the nano-hexapod legs to the various sensors included in the nano-hexapod that could be use for Active Damping, namely:
- A relative motion sensor, measuring the relative displacement of each of the leg
- A force sensor measuring the total force transmitted to the top part of the leg in the direction of the leg
- A absolute velocity sensor measuring the absolute velocity of the top part of the leg in the direction of the leg

After that, a tomography experiment is simulation without any active damping techniques.

** Matlab Init                                              :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
  <<matlab-dir>>
#+end_src

#+begin_src matlab :exports none :results silent :noweb yes
  <<matlab-init>>
#+end_src

#+begin_src matlab :tangle no
  simulinkproject('../');
#+end_src

#+begin_src matlab
  addpath('active_damping/src/');
#+end_src

#+begin_src matlab
  open('active_damping/matlab/sim_nass_active_damping.slx')
#+end_src

** Identification of the dynamics for Active Damping
*** Initialize the Simulation
We initialize all the stages with the default parameters.
#+begin_src matlab
  initializeGround();
  initializeGranite();
  initializeTy();
  initializeRy();
  initializeRz();
  initializeMicroHexapod();
  initializeAxisc();
  initializeMirror();
#+end_src

The nano-hexapod is a piezoelectric hexapod and the sample has a mass of 50kg.
#+begin_src matlab
  initializeNanoHexapod('actuator', 'piezo');
  initializeSample('mass', 50);
#+end_src

We set the references to zero.
#+begin_src matlab
  initializeReferences();
#+end_src

No disturbance is included in the system.
#+begin_src matlab
  initializeDisturbances('enable', false);
#+end_src

And all the controllers are set to 0.
#+begin_src matlab
  K = tf(zeros(6));
  save('./mat/controllers.mat', 'K', '-append');
  K_ine = tf(zeros(6));
  save('./mat/controllers.mat', 'K_ine', '-append');
  K_iff = tf(zeros(6));
  save('./mat/controllers.mat', 'K_iff', '-append');
  K_dvf = tf(zeros(6));
  save('./mat/controllers.mat', 'K_dvf', '-append');
#+end_src

*** Identification
First, we identify the dynamics of the system using the =linearize= function.
#+begin_src matlab
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'sim_nass_active_damping';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fnl'],           1, 'openinput');               io_i = io_i + 1; % Actuator Inputs
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Dnlm');  io_i = io_i + 1; % Relative Motion Outputs
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Fnlm');  io_i = io_i + 1; % Force Sensors
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Vlm');   io_i = io_i + 1; % Absolute Velocity Outputs
  io(io_i) = linio([mdl, '/Compute Error in NASS base'], 2, 'openoutput'); io_i = io_i + 1; % Metrology Outputs

  %% Run the linearization
  G = linearize(mdl, io, 0.5, options);
  G.InputName  = {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'};
  G.OutputName = {'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6', ...
                  'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6', ...
                  'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6', ...
                  'Dxn', 'Dyn', 'Dzn', 'Rxn', 'Ryn', 'Rzn'};
#+end_src

We then create transfer functions corresponding to the active damping plants.
#+begin_src matlab
  G_iff = minreal(G({'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}));
  G_dvf = minreal(G({'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}));
  G_ine = minreal(G({'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}));
#+end_src

#+begin_src matlab
  load('mat/stages.mat', 'nano_hexapod');
  G_cart = minreal(G({'Dxn', 'Dyn', 'Dzn', 'Rxn', 'Ryn', 'Rzn'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))*inv(nano_hexapod.J');
  G_cart.InputName = {'Fnx', 'Fny', 'Fnz', 'Mnx', 'Mny', 'Mnz'};
#+end_src

And we save them for further analysis.
#+begin_src matlab
  save('./active_damping/mat/undamped_plants.mat', 'G_iff', 'G_dvf', 'G_ine');
  save('./active_damping/mat/cart_plants.mat', 'G_cart');
#+end_src

*** Obtained Plants for Active Damping
#+begin_src matlab
  load('./active_damping/mat/undamped_plants.mat', 'G_iff', 'G_dvf', 'G_ine');
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:6
    plot(freqs, abs(squeeze(freqresp(G_iff(['Fnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:6
    plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff(['Fnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_active_damping_iff_plant.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_active_damping_iff_plant
#+CAPTION: =G_iff=: Transfer functions from forces applied in the actuators to the force sensor in each actuator ([[./figs/nass_active_damping_iff_plant.png][png]], [[./figs/nass_active_damping_iff_plant.pdf][pdf]])
[[file:figs/nass_active_damping_iff_plant.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:6
    plot(freqs, abs(squeeze(freqresp(G_dvf(['Dnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:6
    plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf(['Dnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_active_damping_dvf_plant.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_active_damping_dvf_plant
#+CAPTION: =G_dvf=: Transfer functions from forces applied in the actuators to the relative motion sensor in each actuator ([[./figs/nass_active_damping_dvf_plant.png][png]], [[./figs/nass_active_damping_dvf_plant.pdf][pdf]])
[[file:figs/nass_active_damping_dvf_plant.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:6
    plot(freqs, abs(squeeze(freqresp(G_ine(['Vnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [$\frac{m/s}{N}$]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:6
    plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine(['Vnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_active_damping_inertial_plant.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_active_damping_inertial_plant
#+CAPTION: =G_ine=: Transfer functions from forces applied in the actuators to the geophone located in each leg measuring the absolute velocity of the top part of the leg in the direction of the leg ([[./figs/nass_active_damping_inertial_plant.png][png]], [[./figs/nass_active_damping_inertial_plant.pdf][pdf]])
[[file:figs/nass_active_damping_inertial_plant.png]]

** Tomography Experiment
*** Simulation
We initialize elements for the tomography experiment.
#+begin_src matlab
  prepareTomographyExperiment();
#+end_src

We change the simulation stop time.
#+begin_src matlab
  load('mat/conf_simscape.mat');
  set_param(conf_simscape, 'StopTime', '3');
#+end_src

And we simulate the system.
#+begin_src matlab
  sim('sim_nass_active_damping');
#+end_src

Finally, we save the simulation results for further analysis
#+begin_src matlab
  save('./active_damping/mat/tomo_exp.mat', 'En', 'Eg', '-append');
#+end_src

*** Results
We load the results of tomography experiments.
#+begin_src matlab
  load('./active_damping/mat/tomo_exp.mat', 'En');
  t = linspace(0, 3, length(En(:,1)));
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(t, En(:,1), 'DisplayName', '$\epsilon_{x}$')
  plot(t, En(:,2), 'DisplayName', '$\epsilon_{y}$')
  plot(t, En(:,3), 'DisplayName', '$\epsilon_{z}$')
  hold off;
  legend();
  xlabel('Time [s]'); ylabel('Position Error [m]');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_act_damp_undamped_sim_tomo_trans.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_act_damp_undamped_sim_tomo_trans
#+CAPTION: Position Error during tomography experiment - Translations ([[./figs/nass_act_damp_undamped_sim_tomo_trans.png][png]], [[./figs/nass_act_damp_undamped_sim_tomo_trans.pdf][pdf]])
[[file:figs/nass_act_damp_undamped_sim_tomo_trans.png]]

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(t, En(:,4), 'DisplayName', '$\epsilon_{\theta_x}$')
  plot(t, En(:,5), 'DisplayName', '$\epsilon_{\theta_y}$')
  plot(t, En(:,6), 'DisplayName', '$\epsilon_{\theta_z}$')
  hold off;
  xlim([0.5,inf]);
  legend();
  xlabel('Time [s]'); ylabel('Position Error [rad]');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_act_damp_undamped_sim_tomo_rot.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_act_damp_undamped_sim_tomo_rot
#+CAPTION: Position Error during tomography experiment - Rotations ([[./figs/nass_act_damp_undamped_sim_tomo_rot.png][png]], [[./figs/nass_act_damp_undamped_sim_tomo_rot.pdf][pdf]])
[[file:figs/nass_act_damp_undamped_sim_tomo_rot.png]]

* Variability of the system dynamics for Active Damping
:PROPERTIES:
:header-args:matlab+: :tangle matlab/act_damp_variability_plant.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
<<sec:act_damp_variability_plant>>

** Introduction                                                      :ignore:
The goal of this section is to study how the dynamics of the Active Damping plants are changing with the experimental conditions.
These experimental conditions are:
- The mass of the sample (section [[sec:variability_sample_mass]])
- The spindle angle with a null rotating speed (section [[sec:variability_spindle_angle]])
- The spindle rotation speed (section [[sec:variability_rotation_speed]])
- The tilt angle (section [[sec:variability_tilt_angle]])
- The scans of the translation stage (section [[sec:variability_ty_scans]])

For the identification of the dynamics, the system is simulation for $\approx 0.5s$ before the linearization is performed.
This is done in order for the transient phase to be over.

** Matlab Init                                              :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
  <<matlab-dir>>
#+end_src

#+begin_src matlab :exports none :results silent :noweb yes
  <<matlab-init>>
#+end_src

#+begin_src matlab :tangle no
  simulinkproject('../');
  addpath('active_damping/src/');
#+end_src

#+begin_src matlab
  open('active_damping/matlab/sim_nass_active_damping.slx')
  load('mat/conf_simscape.mat');
#+end_src

** Variation of the Sample Mass
<<sec:variability_sample_mass>>
*** Introduction                                                    :ignore:
For all the identifications, the disturbances are disabled and no controller are used.

*** Initialize the Simulation                                     :noexport:
We initialize all the stages with the default parameters.
#+begin_src matlab
  initializeGround();
  initializeGranite();
  initializeTy();
  initializeRy();
  initializeRz();
  initializeMicroHexapod();
  initializeAxisc();
  initializeMirror();
#+end_src

No disturbances.
#+begin_src matlab
  initializeDisturbances('enable', false);
#+end_src

The nano-hexapod is a piezoelectric hexapod.
#+begin_src matlab
  initializeNanoHexapod('actuator', 'piezo');
#+end_src

We set the references to zero.
#+begin_src matlab
  initializeReferences();
#+end_src

And all the controllers are set to 0.
#+begin_src matlab
  K = tf(zeros(6));
  save('./mat/controllers.mat', 'K', '-append');
  K_ine = tf(zeros(6));
  save('./mat/controllers.mat', 'K_ine', '-append');
  K_iff = tf(zeros(6));
  save('./mat/controllers.mat', 'K_iff', '-append');
  K_dvf = tf(zeros(6));
  save('./mat/controllers.mat', 'K_dvf', '-append');
#+end_src

*** Identification                                                  :ignore:
#+begin_src matlab :exports none
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'sim_nass_active_damping';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fnl'],           1, 'openinput');              io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Dnlm'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Fnlm'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Vlm');  io_i = io_i + 1;
#+end_src

We identify the dynamics for the following sample mass.
#+begin_src matlab
  masses = [1, 10, 50]; % [kg]
#+end_src

#+begin_src matlab :exports none
  Gm     = {zeros(length(masses))};
  Gm_iff = {zeros(length(masses))};
  Gm_dvf = {zeros(length(masses))};
  Gm_ine = {zeros(length(masses))};

  for i = 1:length(masses)
    initializeSample('mass', masses(i));

    %% Run the linearization
    G = linearize(mdl, io, 0.3, options);
    G.InputName  = {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'};
    G.OutputName = {'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6', ...
                    'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6', ...
                    'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'};
    Gm(i) = {G};
    Gm_iff(i) = {minreal(G({'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))};
    Gm_dvf(i) = {minreal(G({'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))};
    Gm_ine(i) = {minreal(G({'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))};
  end
#+end_src

#+begin_src matlab :exports none
  save('./active_damping/mat/plants_variable.mat', 'masses', 'Gm_iff', 'Gm_dvf', 'Gm_ine', '-append');
#+end_src

*** Plots                                                           :ignore:
#+begin_src matlab :exports none
  load('./active_damping/mat/plants_variable.mat', 'masses', 'Gm_iff', 'Gm_dvf', 'Gm_ine');
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:length(Gm_iff)
    plot(freqs, abs(squeeze(freqresp(Gm_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Gm_iff)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gm_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))), ...
         'DisplayName', sprintf('$M = %.0f$ [kg]', masses(i)));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_iff_sample_mass.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_iff_sample_mass
#+CAPTION: Variability of the IFF plant with the Sample Mass ([[./figs/act_damp_variability_iff_sample_mass.png][png]], [[./figs/act_damp_variability_iff_sample_mass.pdf][pdf]])
[[file:figs/act_damp_variability_iff_sample_mass.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;

  for i = 1:length(Gm_dvf)
    plot(freqs, abs(squeeze(freqresp(Gm_dvf{i}('Dnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Gm_dvf)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gm_dvf{i}('Dnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$M = %.0f$ [kg]', masses(i)));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_dvf_sample_mass.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_dvf_sample_mass
#+CAPTION: Variability of the DVF plant with the Sample Mass ([[./figs/act_damp_variability_dvf_sample_mass.png][png]], [[./figs/act_damp_variability_dvf_sample_mass.pdf][pdf]])
[[file:figs/act_damp_variability_dvf_sample_mass.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:length(Gm_ine)
    plot(freqs, abs(squeeze(freqresp(Gm_ine{i}('Vnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [$\frac{m/s}{N}$]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Gm_ine)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gm_ine{i}('Vnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$M = %.0f$ [kg]', masses(i)));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_ine_sample_mass.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_ine_sample_mass
#+CAPTION: Variability of the Inertial plant with the Sample Mass ([[./figs/act_damp_variability_ine_sample_mass.png][png]], [[./figs/act_damp_variability_ine_sample_mass.pdf][pdf]])
[[file:figs/act_damp_variability_ine_sample_mass.png]]

** Variation of the Spindle Angle
<<sec:variability_spindle_angle>>
*** Introduction                                                    :ignore:
*** Initialize the Simulation                                     :noexport:
We initialize all the stages with the default parameters.
#+begin_src matlab
  initializeGround();
  initializeGranite();
  initializeTy();
  initializeRy();
  initializeRz();
  initializeMicroHexapod();
  initializeAxisc();
  initializeMirror();
#+end_src

No disturbances.
#+begin_src matlab
  initializeDisturbances('enable', false);
#+end_src

The nano-hexapod is a piezoelectric hexapod.
#+begin_src matlab
  initializeNanoHexapod('actuator', 'piezo');
  initializeSample('mass', 50);
#+end_src

And all the controllers are set to 0.
#+begin_src matlab
  K = tf(zeros(6));
  save('./mat/controllers.mat', 'K', '-append');
  K_ine = tf(zeros(6));
  save('./mat/controllers.mat', 'K_ine', '-append');
  K_iff = tf(zeros(6));
  save('./mat/controllers.mat', 'K_iff', '-append');
  K_dvf = tf(zeros(6));
  save('./mat/controllers.mat', 'K_dvf', '-append');
#+end_src

*** Identification                                                  :ignore:
#+begin_src matlab :exports none
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'sim_nass_active_damping';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fnl'],           1, 'openinput');              io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Dnlm'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Fnlm'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Vlm');  io_i = io_i + 1;
#+end_src

We identify the dynamics for the following Spindle angles.
#+begin_src matlab
  Rz_amplitudes = [0, pi/4, pi/2, pi]; % [rad]
#+end_src

#+begin_src matlab :exports none
  Ga     = {zeros(length(Rz_amplitudes))};
  Ga_iff = {zeros(length(Rz_amplitudes))};
  Ga_dvf = {zeros(length(Rz_amplitudes))};
  Ga_ine = {zeros(length(Rz_amplitudes))};

  for i = 1:length(Rz_amplitudes)
    initializeReferences('Rz_type', 'constant', 'Rz_amplitude', Rz_amplitudes(i))

    %% Run the linearization
    G = linearize(mdl, io, 0.3, options);
    G.InputName  = {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'};
    G.OutputName = {'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6', ...
                    'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6', ...
                    'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'};
    Ga(i) = {G};
    Ga_iff(i) = {minreal(G({'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))};
    Ga_dvf(i) = {minreal(G({'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))};
    Ga_ine(i) = {minreal(G({'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))};
  end
#+end_src

#+begin_src matlab :exports none
  save('./active_damping/mat/plants_variable.mat', 'Rz_amplitudes', 'Ga_iff', 'Ga_dvf', 'Ga_ine', '-append');
#+end_src

*** Plots                                                           :ignore:
#+begin_src matlab :exports none
  load('./active_damping/mat/plants_variable.mat', 'Rz_amplitudes', 'Ga_iff', 'Ga_dvf', 'Ga_ine');
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:length(Ga_iff)
    plot(freqs, abs(squeeze(freqresp(Ga_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Ga_iff)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Ga_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$Rz = %.0f$ [deg]', Rz_amplitudes(i)*180/pi));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_iff_spindle_angle.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_iff_spindle_angle
#+CAPTION: Variability of the IFF plant with the Spindle Angle ([[./figs/act_damp_variability_iff_spindle_angle.png][png]], [[./figs/act_damp_variability_iff_spindle_angle.pdf][pdf]])
[[file:figs/act_damp_variability_iff_spindle_angle.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;

  for i = 1:length(Ga_dvf)
    plot(freqs, abs(squeeze(freqresp(Ga_dvf{i}('Dnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Ga_dvf)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Ga_dvf{i}('Dnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$Rz = %.0f$ [deg]', Rz_amplitudes(i)*180/pi));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_dvf_spindle_angle.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_dvf_spindle_angle
#+CAPTION: Variability of the DVF plant with the Spindle Angle ([[./figs/act_damp_variability_dvf_spindle_angle.png][png]], [[./figs/act_damp_variability_dvf_spindle_angle.pdf][pdf]])
[[file:figs/act_damp_variability_dvf_spindle_angle.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:length(Ga_ine)
    plot(freqs, abs(squeeze(freqresp(Ga_ine{i}('Vnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [$\frac{m/s}{N}$]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Ga_ine)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Ga_ine{i}('Vnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$Rz = %.0f$ [deg]', Rz_amplitudes(i)*180/pi));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_ine_spindle_angle.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_ine_spindle_angle
#+CAPTION: Variability of the Inertial plant with the Spindle Angle ([[./figs/act_damp_variability_ine_spindle_angle.png][png]], [[./figs/act_damp_variability_ine_spindle_angle.pdf][pdf]])
[[file:figs/act_damp_variability_ine_spindle_angle.png]]

** Variation of the Spindle Rotation Speed
<<sec:variability_rotation_speed>>
*** Introduction                                                    :ignore:
*** Initialize the Simulation                                     :noexport:
We initialize all the stages with the default parameters.
#+begin_src matlab
  initializeGround();
  initializeGranite();
  initializeTy();
  initializeRy();
  initializeRz();
  initializeMicroHexapod();
  initializeAxisc();
  initializeMirror();
#+end_src

No disturbances.
#+begin_src matlab
  initializeDisturbances('enable', false);
#+end_src

The nano-hexapod is a piezoelectric hexapod.
#+begin_src matlab
  initializeNanoHexapod('actuator', 'piezo');
  initializeSample('mass', 50);
#+end_src

And all the controllers are set to 0.
#+begin_src matlab
  K = tf(zeros(6));
  save('./mat/controllers.mat', 'K', '-append');
  K_ine = tf(zeros(6));
  save('./mat/controllers.mat', 'K_ine', '-append');
  K_iff = tf(zeros(6));
  save('./mat/controllers.mat', 'K_iff', '-append');
  K_dvf = tf(zeros(6));
  save('./mat/controllers.mat', 'K_dvf', '-append');
#+end_src

*** Identification                                                  :ignore:
#+begin_src matlab :exports none
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'sim_nass_active_damping';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fnl'],           1, 'openinput');              io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Dnlm'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Fnlm'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Vlm');  io_i = io_i + 1;
#+end_src

We identify the dynamics for the following Spindle rotation periods.
#+begin_src matlab
  Rz_periods = [60, 6, 2, 1]; % [s]
#+end_src

The identification of the dynamics is done at the same Spindle angle position.

#+begin_src matlab :exports none
  Gw     = {zeros(length(Rz_periods))};
  Gw_iff = {zeros(length(Rz_periods))};
  Gw_dvf = {zeros(length(Rz_periods))};
  Gw_ine = {zeros(length(Rz_periods))};

  for i = 1:length(Rz_periods)
    initializeReferences('Rz_type', 'rotating', ...
                         'Rz_period', Rz_periods(i), ... % Rotation period [s]
                         'Rz_amplitude', -0.5*(2*pi/Rz_periods(i))); % Angle offset [rad]

    load('mat/nass_references.mat', 'Rz'); % We load the reference for the Spindle
    [~, i_end] = min(abs(Rz.signals.values)); % Obtain the indice where the spindle angle is zero
    t_sim = Rz.time(i_end) % Simulation time before identification [s]

    %% Run the linearization
    G = linearize(mdl, io, t_sim, options);
    G.InputName  = {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'};
    G.OutputName = {'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6', ...
                    'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6', ...
                    'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'};

    Gw(i) = {G};
    Gw_iff(i) = {minreal(G({'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))};
    Gw_dvf(i) = {minreal(G({'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))};
    Gw_ine(i) = {minreal(G({'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))};
  end
#+end_src

#+begin_src matlab :exports none
  save('./active_damping/mat/plants_variable.mat', 'Rz_periods', 'Gw_iff', 'Gw_dvf', 'Gw_ine', '-append');
#+end_src

*** Dynamics of the Active Damping plants
#+begin_src matlab :exports none
  load('./active_damping/mat/plants_variable.mat', 'Rz_periods', 'Gw_iff', 'Gw_dvf', 'Gw_ine');
  load('./active_damping/mat/undamped_plants.mat', 'G_iff', 'G_dvf', 'G_ine');
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 10000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:length(Gw_iff)
    plot(freqs, abs(squeeze(freqresp(Gw_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  plot(freqs, abs(squeeze(freqresp(G_iff('Fnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Gw_iff)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gw_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$Rz = %.0f$ [rpm]', 60/Rz_periods(i)));
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff('Fnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--', 'DisplayName', 'No Rotation');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_iff_spindle_speed.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_iff_spindle_speed
#+CAPTION: Variability of the IFF plant with the Spindle rotation speed ([[./figs/act_damp_variability_iff_spindle_speed.png][png]], [[./figs/act_damp_variability_iff_spindle_speed.pdf][pdf]])
[[file:figs/act_damp_variability_iff_spindle_speed.png]]

#+begin_src matlab :exports none
  xlim([20, 30]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_iff_spindle_speed_zoom.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_iff_spindle_speed_zoom
#+CAPTION: Variability of the IFF plant with the Spindle rotation speed ([[./figs/act_damp_variability_iff_spindle_speed_zoom.png][png]], [[./figs/act_damp_variability_iff_spindle_speed_zoom.pdf][pdf]])
[[file:figs/act_damp_variability_iff_spindle_speed_zoom.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 5000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;

  for i = 1:length(Gw_dvf)
    plot(freqs, abs(squeeze(freqresp(Gw_dvf{i}('Dnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  plot(freqs, abs(squeeze(freqresp(G_dvf('Dnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Gw_dvf)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gw_dvf{i}('Dnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$Rz = %.0f$ [rpm]', 60/Rz_periods(i)));
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf('Dnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--', 'DisplayName', 'No Rotation');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_dvf_spindle_speed.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_dvf_spindle_speed
#+CAPTION: Variability of the DVF plant with the Spindle rotation speed ([[./figs/act_damp_variability_dvf_spindle_speed.png][png]], [[./figs/act_damp_variability_dvf_spindle_speed.pdf][pdf]])
[[file:figs/act_damp_variability_dvf_spindle_speed.png]]

#+begin_src matlab :exports none
  xlim([20, 30]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_dvf_spindle_speed_zoom.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_dvf_spindle_speed_zoom
#+CAPTION: Variability of the DVF plant with the Spindle rotation speed ([[./figs/act_damp_variability_dvf_spindle_speed_zoom.png][png]], [[./figs/act_damp_variability_dvf_spindle_speed_zoom.pdf][pdf]])
[[file:figs/act_damp_variability_dvf_spindle_speed_zoom.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 5000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:length(Gw_ine)
    plot(freqs, abs(squeeze(freqresp(Gw_ine{i}('Vnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  plot(freqs, abs(squeeze(freqresp(G_ine('Vnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [$\frac{m/s}{N}$]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Gw_ine)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gw_ine{i}('Vnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$Rz = %.0f$ [rpm]', 60/Rz_periods(i)));
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine('Vnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--', 'DisplayName', 'No Rotation');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_ine_spindle_speed.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_ine_spindle_speed
#+CAPTION: Variability of the Inertial plant with the Spindle rotation speed ([[./figs/act_damp_variability_ine_spindle_speed.png][png]], [[./figs/act_damp_variability_ine_spindle_speed.pdf][pdf]])
[[file:figs/act_damp_variability_ine_spindle_speed.png]]

#+begin_src matlab :exports none
  xlim([20, 30]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_ine_spindle_speed_zoom.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_ine_spindle_speed_zoom
#+CAPTION: Variability of the Inertial plant with the Spindle rotation speed ([[./figs/act_damp_variability_ine_spindle_speed_zoom.png][png]], [[./figs/act_damp_variability_ine_spindle_speed_zoom.pdf][pdf]])
[[file:figs/act_damp_variability_ine_spindle_speed_zoom.png]]

*** Variation of the poles and zeros with the Spindle rotation frequency
#+begin_src matlab :exports none
  load('./active_damping/mat/plants_variable.mat', 'Rz_periods', 'Gw_iff', 'Gw_dvf', 'Gw_ine');
  load('./active_damping/mat/undamped_plants.mat', 'G_iff', 'G_dvf', 'G_ine');
#+end_src


#+begin_src matlab :exports none
  figure;

  subplot(1,2,1);
  hold on;
  for i = 1:length(Gw_iff)
    G_poles = pole(Gw_iff{i}('Fnlm1', 'Fnl1'));
    plot(1/Rz_periods(i), real(G_poles(imag(G_poles)<2*pi*30 & imag(G_poles)>2*pi*22)), 'kx');
  end
  G_poles = pole(G_iff('Fnlm1', 'Fnl1'));
  plot(0, real(G_poles(imag(G_poles)<2*pi*30 & imag(G_poles)>2*pi*22)), 'kx');
  hold off;
  ylim([-inf, 0]);
  xlabel('Rotation Speed [Hz]');
  ylabel('Real Part');

  subplot(1,2,2);
  hold on;
  for i = 1:length(Gw_iff)
    G_poles = pole(Gw_iff{i}('Fnlm1', 'Fnl1'));
    plot(1/Rz_periods(i), imag(G_poles(imag(G_poles)<2*pi*30 & imag(G_poles)>2*pi*22)), 'kx');
  end
  G_poles = pole(G_iff('Fnlm1', 'Fnl1'));
  plot(0, imag(G_poles(imag(G_poles)<2*pi*30 & imag(G_poles)>2*pi*22)), 'kx');
  hold off;
  ylim([0, inf]);
  xlabel('Rotation Speed [Hz]');
  ylabel('Imaginary Part');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/campbell_diagram_spindle_rotation.pdf" :var figsize="full-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:campbell_diagram_spindle_rotation
#+CAPTION: Evolution of the pole with respect to the spindle rotation speed ([[./figs/campbell_diagram_spindle_rotation.png][png]], [[./figs/campbell_diagram_spindle_rotation.pdf][pdf]])
[[file:figs/campbell_diagram_spindle_rotation.png]]

#+begin_src matlab :exports none
  figure;

  subplot(1,2,1);
  hold on;
  for i = 1:length(Gw_ine)
    set(gca,'ColorOrderIndex',1);
    G_zeros = zero(Gw_ine{i}('Vnlm1', 'Fnl1'));
    plot(1/Rz_periods(i), real(G_zeros(imag(G_zeros)<2*pi*25 & imag(G_zeros)>2*pi*22)), 'o');

    set(gca,'ColorOrderIndex',2);
    G_zeros = zero(Gw_iff{i}('Fnlm1', 'Fnl1'));
    plot(1/Rz_periods(i), real(G_zeros(imag(G_zeros)<2*pi*25 & imag(G_zeros)>2*pi*22)), 'o');

    set(gca,'ColorOrderIndex',3);
    G_zeros = zero(Gw_dvf{i}('Dnlm1', 'Fnl1'));
    plot(1/Rz_periods(i), real(G_zeros(imag(G_zeros)<2*pi*25 & imag(G_zeros)>2*pi*22)), 'o');
  end
  hold off;
  xlabel('Rotation Speed [Hz]');
  ylabel('Real Part');

  subplot(1,2,2);
  hold on;
  for i = 1:length(Gw_ine)
    set(gca,'ColorOrderIndex',1);
    G_zeros = zero(Gw_ine{i}('Vnlm1', 'Fnl1'));
    p_ine = plot(1/Rz_periods(i), imag(G_zeros(imag(G_zeros)<2*pi*25 & imag(G_zeros)>2*pi*22)), 'o');

    set(gca,'ColorOrderIndex',2);
    G_zeros = zero(Gw_iff{i}('Fnlm1', 'Fnl1'));
    p_iff = plot(1/Rz_periods(i), imag(G_zeros(imag(G_zeros)<2*pi*25 & imag(G_zeros)>2*pi*22)), 'o');

    set(gca,'ColorOrderIndex',3);
    G_zeros = zero(Gw_dvf{i}('Dnlm1', 'Fnl1'));
    p_dvf = plot(1/Rz_periods(i), imag(G_zeros(imag(G_zeros)<2*pi*25 & imag(G_zeros)>2*pi*22)), 'o');
  end
  hold off;
  xlabel('Rotation Speed [Hz]');
  ylabel('Imaginary Part');
  legend([p_ine p_iff p_dvf],{'Inertial Sensor','Force Sensor', 'Relative Motion Sensor'}, 'location', 'southwest');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/variation_zeros_active_damping_plants.pdf" :var figsize="full-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:variation_zeros_active_damping_plants
#+CAPTION: Evolution of the zero with respect to the spindle rotation speed ([[./figs/variation_zeros_active_damping_plants.png][png]], [[./figs/variation_zeros_active_damping_plants.pdf][pdf]])
[[file:figs/variation_zeros_active_damping_plants.png]]

** Variation of the Tilt Angle
<<sec:variability_tilt_angle>>
*** Introduction                                                    :ignore:
*** Initialize the Simulation                                     :noexport:
We initialize all the stages with the default parameters.
#+begin_src matlab
  initializeGround();
  initializeGranite();
  initializeTy();
  initializeRy();
  initializeRz();
  initializeMicroHexapod();
  initializeAxisc();
  initializeMirror();
#+end_src

No disturbances.
#+begin_src matlab
  initializeDisturbances('enable', false);
#+end_src

The nano-hexapod is a piezoelectric hexapod.
#+begin_src matlab
  initializeNanoHexapod('actuator', 'piezo');
  initializeSample('mass', 50);
#+end_src

And all the controllers are set to 0.
#+begin_src matlab
  K = tf(zeros(6));
  save('./mat/controllers.mat', 'K', '-append');
  K_ine = tf(zeros(6));
  save('./mat/controllers.mat', 'K_ine', '-append');
  K_iff = tf(zeros(6));
  save('./mat/controllers.mat', 'K_iff', '-append');
  K_dvf = tf(zeros(6));
  save('./mat/controllers.mat', 'K_dvf', '-append');
#+end_src

*** Identification                                                  :ignore:
#+begin_src matlab :exports none
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'sim_nass_active_damping';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fnl'],           1, 'openinput');              io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Dnlm'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Fnlm'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Vlm');  io_i = io_i + 1;
#+end_src

We identify the dynamics for the following Tilt stage angles.
#+begin_src matlab
  Ry_amplitudes = [-3*pi/180, 3*pi/180]; % [rad]
#+end_src

#+begin_src matlab :exports none
  Gy     = {zeros(length(Ry_amplitudes))};
  Gy_iff = {zeros(length(Ry_amplitudes))};
  Gy_dvf = {zeros(length(Ry_amplitudes))};
  Gy_ine = {zeros(length(Ry_amplitudes))};

  for i = 1:length(Ry_amplitudes)
    initializeReferences('Ry_type', 'constant', 'Ry_amplitude', Ry_amplitudes(i))

    %% Run the linearization
    G = linearize(mdl, io, 0.3, options);
    G.InputName  = {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'};
    G.OutputName = {'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6', ...
                    'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6', ...
                    'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'};
    Gy(i) = {G};
    Gy_iff(i) = {minreal(G({'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))};
    Gy_dvf(i) = {minreal(G({'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))};
    Gy_ine(i) = {minreal(G({'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))};
  end
#+end_src

#+begin_src matlab :exports none
  save('./active_damping/mat/plants_variable.mat', 'Ry_amplitudes', 'Gy_iff', 'Gy_dvf', 'Gy_ine', '-append');
#+end_src

*** Plots                                                           :ignore:
#+begin_src matlab :exports none
  load('./active_damping/mat/plants_variable.mat', 'Ry_amplitudes', 'Gy_iff', 'Gy_dvf', 'Gy_ine');
  load('./active_damping/mat/undamped_plants.mat', 'G_iff', 'G_dvf', 'G_ine');
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:length(Gy_iff)
    plot(freqs, abs(squeeze(freqresp(Gy_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  plot(freqs, abs(squeeze(freqresp(G_iff('Fnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Gy_iff)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gy_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$Ry = %.0f$ [deg]', Ry_amplitudes(i)*180/pi));
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff('Fnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--', 'DisplayName', '$Ry = 0$ [deg]');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_iff_tilt_angle.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_iff_tilt_angle
#+CAPTION: Variability of the IFF plant with the Tilt stage Angle ([[./figs/act_damp_variability_iff_tilt_angle.png][png]], [[./figs/act_damp_variability_iff_tilt_angle.pdf][pdf]])
[[file:figs/act_damp_variability_iff_tilt_angle.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;

  for i = 1:length(Gy_dvf)
    plot(freqs, abs(squeeze(freqresp(Gy_dvf{i}('Dnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  plot(freqs, abs(squeeze(freqresp(G_dvf('Dnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Gy_dvf)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gy_dvf{i}('Dnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$Ry = %.0f$ [deg]', Ry_amplitudes(i)*180/pi));
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf('Dnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--', 'DisplayName', '$Ry = 0$ [deg]');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_dvf_tilt_angle.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_dvf_tilt_angle
#+CAPTION: Variability of the DVF plant with the Tilt Angle ([[./figs/act_damp_variability_dvf_tilt_angle.png][png]], [[./figs/act_damp_variability_dvf_tilt_angle.pdf][pdf]])
[[file:figs/act_damp_variability_dvf_tilt_angle.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:length(Gy_ine)
    plot(freqs, abs(squeeze(freqresp(Gy_ine{i}('Vnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  plot(freqs, abs(squeeze(freqresp(G_ine('Vnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [$\frac{m/s}{N}$]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Gy_ine)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gy_ine{i}('Vnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$Ry = %.0f$ [deg]', Ry_amplitudes(i)*180/pi));
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine('Vnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--', 'DisplayName', '$Ry = 0$ [deg]');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_ine_tilt_angle.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_ine_tilt_angle
#+CAPTION: Variability of the Inertial plant with the Tilt Angle ([[./figs/act_damp_variability_ine_tilt_angle.png][png]], [[./figs/act_damp_variability_ine_tilt_angle.pdf][pdf]])
[[file:figs/act_damp_variability_ine_tilt_angle.png]]

** Scans of the Translation Stage
<<sec:variability_ty_scans>>
*** Introduction                                                    :ignore:
We want here to verify if the dynamics used for Active damping is varying when using the translation stage for scans.

*** Initialize the Simulation                                     :noexport:
We initialize all the stages with the default parameters.
#+begin_src matlab
  initializeGround();
  initializeGranite();
  initializeTy();
  initializeRy();
  initializeRz();
  initializeMicroHexapod();
  initializeAxisc();
  initializeMirror();
#+end_src

No disturbances.
#+begin_src matlab
  initializeDisturbances('enable', false);
#+end_src

The nano-hexapod is a piezoelectric hexapod and the mass of the sample is set to 50kg.
#+begin_src matlab
  initializeNanoHexapod('actuator', 'piezo');
  initializeSample('mass', 50);
#+end_src

And all the controllers are set to 0.
#+begin_src matlab
  K = tf(zeros(6));
  save('./mat/controllers.mat', 'K', '-append');
  K_ine = tf(zeros(6));
  save('./mat/controllers.mat', 'K_ine', '-append');
  K_iff = tf(zeros(6));
  save('./mat/controllers.mat', 'K_iff', '-append');
  K_dvf = tf(zeros(6));
  save('./mat/controllers.mat', 'K_dvf', '-append');
#+end_src

*** Identification                                                  :ignore:
#+begin_src matlab :exports none
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'sim_nass_active_damping';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fnl'],           1, 'openinput');              io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Dnlm'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Fnlm'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Vlm');  io_i = io_i + 1;
#+end_src

We initialize the translation stage reference to be a sinus with an amplitude of 5mm and a period of 1s (Figure [[fig:ty_scanning_reference_sinus]]).
#+begin_src matlab
  initializeReferences('Dy_type', 'sinusoidal', ...
                       'Dy_amplitude', 5e-3, ... % [m]
                       'Dy_period', 1); % [s]
#+end_src

#+begin_src matlab :exports none
  load('mat/nass_references.mat', 'Dy');
  figure;
  plot(Dy.time, Dy.signals.values);
  xlabel('Time [s]'); ylabel('Dy - Position [m]');
  xlim([0, 2]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/ty_scanning_reference_sinus.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:ty_scanning_reference_sinus
#+CAPTION: Reference path for the translation stage ([[./figs/ty_scanning_reference_sinus.png][png]], [[./figs/ty_scanning_reference_sinus.pdf][pdf]])
[[file:figs/ty_scanning_reference_sinus.png]]

We identify the dynamics at different positions (times) when scanning with the Translation stage.
#+begin_src matlab
  t_lin = [0.5, 0.75, 1, 1.25];
#+end_src

#+begin_src matlab :exports none
  Gty     = {zeros(length(t_lin))};
  Gty_iff = {zeros(length(t_lin))};
  Gty_dvf = {zeros(length(t_lin))};
  Gty_ine = {zeros(length(t_lin))};

  %% Run the linearization
  G = linearize(mdl, io, t_lin, options);
  G.InputName  = {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'};
  G.OutputName = {'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6', ...
                  'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6', ...
                  'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'};

  for i = 1:length(t_lin)
    Gty(i) = {G(:,:,i)};
    Gty_iff(i) = {minreal(G({'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}, i))};
    Gty_dvf(i) = {minreal(G({'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}, i))};
    Gty_ine(i) = {minreal(G({'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}, i))};
  end
#+end_src

#+begin_src matlab :exports none
  Gty_tlin = t_lin;
  save('./active_damping/mat/plants_variable.mat', 'Gty_tlin', 'Dy', 'Gty_iff', 'Gty_dvf', 'Gty_ine', '-append');
#+end_src

*** Plots                                                           :ignore:
#+begin_src matlab :exports none
  load('./active_damping/mat/plants_variable.mat', 'Gty_tlin', 'Dy', 'Gty_iff', 'Gty_dvf', 'Gty_ine');
  load('./active_damping/mat/undamped_plants.mat', 'G_iff', 'G_dvf', 'G_ine');
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:length(Gty_iff)
    plot(freqs, abs(squeeze(freqresp(Gty_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  plot(freqs, abs(squeeze(freqresp(G_iff('Fnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Gty_iff)
    [~, i_t] = min(abs(Dy.time - Gty_tlin(i)));
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gty_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$Dy = %.0f$ [mm]', 1e3*Dy.signals.values(i_t)));
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff('Fnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--', 'DisplayName', '$Ry = 0$ [deg]');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_iff_ty_scans.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_iff_ty_scans
#+CAPTION: Variability of the IFF plant at different Ty scan positions ([[./figs/act_damp_variability_iff_ty_scans.png][png]], [[./figs/act_damp_variability_iff_ty_scans.pdf][pdf]])
[[file:figs/act_damp_variability_iff_ty_scans.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;

  for i = 1:length(Gty_dvf)
    plot(freqs, abs(squeeze(freqresp(Gty_dvf{i}('Dnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  plot(freqs, abs(squeeze(freqresp(G_dvf('Dnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Gty_dvf)
    [~, i_t] = min(abs(Dy.time - Gty_tlin(i)));
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gty_dvf{i}('Dnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$Dy = %.0f$ [mm]', 1e3*Dy.signals.values(i_t)));
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf('Dnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--', 'DisplayName', '$Ry = 0$ [deg]');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_dvf_ty_scans.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_dvf_ty_scans
#+CAPTION: Variability of the DVF plant at different Ty scan positions ([[./figs/act_damp_variability_dvf_ty_scans.png][png]], [[./figs/act_damp_variability_dvf_ty_scans.pdf][pdf]])
[[file:figs/act_damp_variability_dvf_ty_scans.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i = 1:length(Gty_ine)
    plot(freqs, abs(squeeze(freqresp(Gty_ine{i}('Vnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  plot(freqs, abs(squeeze(freqresp(G_ine('Vnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [$\frac{m/s}{N}$]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i = 1:length(Gty_ine)
    [~, i_t] = min(abs(Dy.time - Gty_tlin(i)));
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gty_ine{i}('Vnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$Dy = %.0f$ [mm]', 1e3*Dy.signals.values(i_t)));
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine('Vnlm1', 'Fnl1'), freqs, 'Hz'))), 'k--', 'DisplayName', '$Ry = 0$ [deg]');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_variability_ine_ty_scans.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_variability_ine_ty_scans
#+CAPTION: Variability of the Inertial plant at different Ty scan positions ([[./figs/act_damp_variability_ine_ty_scans.png][png]], [[./figs/act_damp_variability_ine_ty_scans.pdf][pdf]])
[[file:figs/act_damp_variability_ine_ty_scans.png]]

** Conclusion

#+name: tab:active_damping_plant_variability
#+caption: Conclusion on the variability of the system dynamics for active damping
|            <c>            |                    <c>                     |
|                           |            *Change of Dynamics*            |
|---------------------------+--------------------------------------------|
|       *Sample Mass*       |                   Large                    |
|      *Spindle Angle*      |                   Small                    |
| *Spindle Rotation Speed*  | First resonance is split in two resonances |
|       *Tilt Angle*        |                 Negligible                 |
| *Translation Stage scans* |                 Negligible                 |

Also, for the Inertial Sensor, a RHP zero may appear when the spindle is rotating fast.

#+begin_important
  When using either a force sensor or a relative motion sensor for active damping, the only "parameter" that induce a large change for the dynamics to be controlled is the *sample mass*.
  Thus, the developed damping techniques should be robust to variations of the sample mass.
#+end_important

* Integral Force Feedback
:PROPERTIES:
:header-args:matlab+: :tangle matlab/iff.m
:header-args:matlab+: :comments none :mkdirp yes
:END:
<<sec:iff>>

** ZIP file containing the data and matlab files                     :ignore:
#+begin_src bash :exports none :results none
  if [ matlab/iff.m -nt data/iff.zip ]; then
    cp matlab/iff.m iff.m;
    zip data/iff \
        mat/plant.mat \
        iff.m
    rm iff.m;
  fi
#+end_src

#+begin_note
  All the files (data and Matlab scripts) are accessible [[file:data/iff.zip][here]].
#+end_note

** Introduction                                                      :ignore:
Here, we study the use of *Integral Force Feedback* (IFF) to actively damp the resonances.

The IFF control is applied in a decentralized way: there is on controller for each leg.

The control architecture is represented in figure [[fig:iff_1dof]] where one of the 6 nano-hexapod legs is represented.

#+name: fig:iff_1dof
#+caption: Integral Force Feedback applied to a 1dof system
#+RESULTS:
[[file:figs/iff_1dof.png]]

** Matlab Init                                              :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
  <<matlab-dir>>
#+end_src

#+begin_src matlab :exports none :results silent :noweb yes
  <<matlab-init>>
#+end_src

#+begin_src matlab :tangle no
  simulinkproject('../');
#+end_src

#+begin_src matlab
  addpath('active_damping/src/');
#+end_src

#+begin_src matlab
  open('active_damping/matlab/sim_nass_active_damping.slx')
#+end_src

** Control Design
*** Plant
Let's load the previously identified undamped plant:
#+begin_src matlab
  load('./active_damping/mat/undamped_plants.mat', 'G_iff');
  load('./active_damping/mat/plants_variable.mat', 'masses', 'Gm_iff');
#+end_src

Let's look at the transfer function from actuator forces in the nano-hexapod to the force sensor in the nano-hexapod legs for all 6 pairs of actuator/sensor (figure [[fig:iff_plant]]).

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i=1:length(masses)
    plot(freqs, abs(squeeze(freqresp(-Gm_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i=1:length(masses)
    plot(freqs, 180/pi*angle(squeeze(freqresp(-Gm_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))), ...
         'DisplayName', sprintf('$M = %.0f$ [kg]', masses(i)));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/iff_plant.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:iff_plant
#+CAPTION: Transfer function from forces applied in the legs to force sensor ([[./figs/iff_plant.png][png]], [[./figs/iff_plant.pdf][pdf]])
[[file:figs/iff_plant.png]]

*** Control Design
The controller for each pair of actuator/sensor is:
#+begin_src matlab
  w0 = 2*pi*50;
  K_iff = -5000/s * (s/w0)/(1 + s/w0);
#+end_src

The corresponding loop gains are shown in figure [[fig:iff_open_loop]].

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i=1:length(masses)
    plot(freqs, abs(squeeze(freqresp(K_iff*Gm_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i=1:length(masses)
    plot(freqs, 180/pi*angle(squeeze(freqresp(K_iff*Gm_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))), ...
         'DisplayName', sprintf('$M = %.0f$ [kg]', masses(i)));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/iff_open_loop.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:iff_open_loop
#+CAPTION: Loop Gain for the Integral Force Feedback ([[./figs/iff_open_loop.png][png]], [[./figs/iff_open_loop.pdf][pdf]])
[[file:figs/iff_open_loop.png]]

*** Diagonal Controller
We create the diagonal controller and we add a minus sign as we have a positive
feedback architecture.
#+begin_src matlab
  K_iff = -K_iff*eye(6);
#+end_src

We save the controller for further analysis.
#+begin_src matlab
  save('./active_damping/mat/K_iff.mat', 'K_iff');
#+end_src

** TODO Identification of the damped plant                         :noexport:
*** Initialize the Simulation                                     :noexport:
We initialize all the stages with the default parameters.
#+begin_src matlab
  initializeGround();
  initializeGranite();
  initializeTy();
  initializeRy();
  initializeRz();
  initializeMicroHexapod();
  initializeAxisc();
  initializeMirror();
#+end_src

No disturbances.
#+begin_src matlab
  initializeDisturbances('enable', false);
#+end_src

The nano-hexapod is a piezoelectric hexapod.
#+begin_src matlab
  initializeNanoHexapod('actuator', 'piezo');
  initializeSample('mass', 50);
#+end_src

And all the controllers are set to 0.
#+begin_src matlab
  K = tf(zeros(6));
  save('./mat/controllers.mat', 'K', '-append');
  K_ine = tf(zeros(6));
  save('./mat/controllers.mat', 'K_ine', '-append');
  load('./active_damping/mat/K_iff.mat', 'K_iff');
  save('./mat/controllers.mat', 'K_iff', '-append');
  K_dvf = tf(zeros(6));
  save('./mat/controllers.mat', 'K_dvf', '-append');
#+end_src

*** Identification
First, we identify the dynamics of the system using the =linearize= function.
#+begin_src matlab
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'sim_nass_active_damping';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fnl'],           1, 'openinput');               io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Compute Error in NASS base'], 2, 'openoutput'); io_i = io_i + 1;

  %% Run the linearization
  G = linearize(mdl, io, 0.5, options);
  G.InputName  = {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'};
  G.OutputName = {'Dxn', 'Dyn', 'Dzn', 'Rxn', 'Ryn', 'Rzn'};
#+end_src

#+begin_src matlab
  load('mat/stages.mat', 'nano_hexapod');
  G_cart_iff = G*inv(nano_hexapod.J');
  G_cart_iff.InputName = {'Fnx', 'Fny', 'Fnz', 'Mnx', 'Mny', 'Mnz'};
#+end_src

And we save them for further analysis.
#+begin_src matlab
  save('./active_damping/mat/cart_plants.mat', 'G_cart_iff', '-append');
#+end_src

*** TODO Sensitivity to disturbances
As shown on figure [[fig:sensitivity_dist_iff]]:
- The top platform of the nano-hexapod how behaves as a "free-mass".
- The transfer function from direct forces $F_s$ to the relative displacement $D$ is equivalent to the one of an isolated mass.
- The transfer function from ground motion $D_g$ to the relative displacement $D$ tends to the transfer function from $D_g$ to the displacement of the granite (the sample is being isolated thanks to IFF).
  However, as the goal is to make the relative displacement $D$ as small as possible (e.g. to make the sample motion follows the granite motion), this is not a good thing.

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  subplot(2, 1, 1);
  title('$D_g$ to $D$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_gm('Dx', 'Dgx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / D_{g,x}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_gm('Dy', 'Dgy'), freqs, 'Hz'))), 'DisplayName', '$\left|D_y / D_{g,y}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_gm('Dz', 'Dgz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / D_{g,z}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_iff.G_gm('Dx', 'Dgx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_iff.G_gm('Dy', 'Dgy'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_iff.G_gm('Dz', 'Dgz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');

  subplot(2, 1, 2);
  title('$F_s$ to $D$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_fs('Dx', 'Fsx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / F_{s,x}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_fs('Dy', 'Fsy'), freqs, 'Hz'))), 'DisplayName', '$\left|D_y / F_{s,y}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_fs('Dz', 'Fsz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{s,z}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_iff.G_fs('Dx', 'Fsx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_iff.G_fs('Dy', 'Fsy'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_iff.G_fs('Dz', 'Fsz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_dist_iff.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_dist_iff
#+CAPTION: Sensitivity to disturbance once the IFF controller is applied to the system ([[./figs/sensitivity_dist_iff.png][png]], [[./figs/sensitivity_dist_iff.pdf][pdf]])
[[file:figs/sensitivity_dist_iff.png]]

#+begin_warning
  The order of the models are very high and thus the plots may be wrong.
  For instance, the plots are not the same when using =minreal=.
#+end_warning

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_dist('Dz', 'Frzz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{rz, z}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_dist('Dz', 'Ftyz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{ty, z}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_dist('Dx', 'Ftyx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / F_{ty, x}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(minreal(prescale(G_iff.G_dist('Dz', 'Frzz'), {2*pi, 2*pi*1e3})), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(minreal(G_iff.G_dist('Dz', 'Ftyz')), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(minreal(G_iff.G_dist('Dx', 'Ftyx')), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_dist_stages_iff.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_dist_stages_iff
#+CAPTION: Sensitivity to force disturbances in various stages when IFF is applied ([[./figs/sensitivity_dist_stages_iff.png][png]], [[./figs/sensitivity_dist_stages_iff.pdf][pdf]])
[[file:figs/sensitivity_dist_stages_iff.png]]

*** Damped Plant
Now, look at the new damped plant to control.
#+begin_src matlab
  load('./active_damping/mat/cart_plants.mat', 'G_cart', 'G_cart_iff');
#+end_src

It damps the plant (resonance of the nano hexapod as well as other resonances) as shown in figure [[fig:plant_iff_damped]].

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 2, 1);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_cart('Dxn', 'Fnx'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G_cart('Dyn', 'Fny'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G_cart('Dzn', 'Fnz'), freqs, 'Hz'))));
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_cart_iff('Dxn', 'Fnx'), freqs, 'Hz'))), '--');
  plot(freqs, abs(squeeze(freqresp(G_cart_iff('Dyn', 'Fny'), freqs, 'Hz'))), '--');
  plot(freqs, abs(squeeze(freqresp(G_cart_iff('Dzn', 'Fnz'), freqs, 'Hz'))), '--');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');

  ax2 = subplot(2, 2, 2);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_cart('Rxn', 'Mnx'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G_cart('Ryn', 'Mny'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G_cart('Rzn', 'Mnz'), freqs, 'Hz'))));
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_cart_iff('Rxn', 'Mnx'), freqs, 'Hz'))), '--');
  plot(freqs, abs(squeeze(freqresp(G_cart_iff('Ryn', 'Mny'), freqs, 'Hz'))), '--');
  plot(freqs, abs(squeeze(freqresp(G_cart_iff('Rzn', 'Mnz'), freqs, 'Hz'))), '--');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [rad/(Nm)]'); xlabel('Frequency [Hz]');

  ax3 = subplot(2, 2, 3);
  hold on;
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart('Dxn', 'Fnx'), freqs, 'Hz')))), 'DisplayName', '$\left|D_x / F_{n,x}\right|$');
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart('Dyn', 'Fny'), freqs, 'Hz')))), 'DisplayName', '$\left|D_y / F_{n,y}\right|$');
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart('Dzn', 'Fnz'), freqs, 'Hz')))), 'DisplayName', '$\left|D_z / F_{n,z}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart_iff('Dxn', 'Fnx'), freqs, 'Hz')))), '--', 'HandleVisibility', 'off');
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart_iff('Dyn', 'Fny'), freqs, 'Hz')))), '--', 'HandleVisibility', 'off');
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart_iff('Dzn', 'Fnz'), freqs, 'Hz')))), '--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  yticks([-540:180:540]);
  legend('location', 'southwest');

  ax4 = subplot(2, 2, 4);
  hold on;
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart('Rxn', 'Mnx'), freqs, 'Hz')))), 'DisplayName', '$\left|R_x / M_{n,x}\right|$');
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart('Ryn', 'Mny'), freqs, 'Hz')))), 'DisplayName', '$\left|R_y / M_{n,y}\right|$');
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart('Rzn', 'Mnz'), freqs, 'Hz')))), 'DisplayName', '$\left|R_z / M_{n,z}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart_iff('Rxn', 'Mnx'), freqs, 'Hz')))), '--', 'HandleVisibility', 'off');
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart_iff('Ryn', 'Mny'), freqs, 'Hz')))), '--', 'HandleVisibility', 'off');
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart_iff('Rzn', 'Mnz'), freqs, 'Hz')))), '--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  yticks([-540:180:540]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2,ax3,ax4],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/plant_iff_damped.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:plant_iff_damped
#+CAPTION: Damped Plant after IFF is applied ([[./figs/plant_iff_damped.png][png]], [[./figs/plant_iff_damped.pdf][pdf]])
[[file:figs/plant_iff_damped.png]]

However, it increases coupling at low frequency (figure [[fig:plant_iff_coupling]]).
#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  for ix = 1:6
    for iy = 1:6
      subplot(6, 6, (ix-1)*6 + iy);
      hold on;
      plot(freqs, abs(squeeze(freqresp(G_cart(ix, iy), freqs, 'Hz'))), 'k-');
      plot(freqs, abs(squeeze(freqresp(G_cart_iff(ix, iy), freqs, 'Hz'))), 'k--');
      set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
      ylim([1e-13, 1e-4]);
      xticks([])
      yticks([])
    end
  end
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/plant_iff_coupling.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:plant_iff_coupling
#+CAPTION: Coupling induced by IFF ([[./figs/plant_iff_coupling.png][png]], [[./figs/plant_iff_coupling.pdf][pdf]])
[[file:figs/plant_iff_coupling.png]]

** Tomography Experiment
*** Simulation with IFF Controller
We initialize elements for the tomography experiment.
#+begin_src matlab
  prepareTomographyExperiment();
#+end_src

We set the IFF controller.
#+begin_src matlab
  load('./active_damping/mat/K_iff.mat', 'K_iff');
  save('./mat/controllers.mat', 'K_iff', '-append');
#+end_src

We change the simulation stop time.
#+begin_src matlab
  load('mat/conf_simscape.mat');
  set_param(conf_simscape, 'StopTime', '3');
#+end_src

And we simulate the system.
#+begin_src matlab
  sim('sim_nass_active_damping');
#+end_src

Finally, we save the simulation results for further analysis
#+begin_src matlab
  En_iff = En;
  Eg_iff = Eg;
  save('./active_damping/mat/tomo_exp.mat', 'En_iff', 'Eg_iff', '-append');
#+end_src

*** Simulation with IFF Controller with added High Pass Filter
We initialize elements for the tomography experiment.
#+begin_src matlab
  prepareTomographyExperiment();
#+end_src

We set the IFF controller with the High Pass Filter.
#+begin_src matlab
  load('./active_damping/mat/K_iff_hpf.mat', 'K_iff');
  save('./mat/controllers.mat', 'K_iff', '-append');
#+end_src

We change the simulation stop time.
#+begin_src matlab
  load('mat/conf_simscape.mat');
  set_param(conf_simscape, 'StopTime', '3');
#+end_src

And we simulate the system.
#+begin_src matlab
  sim('sim_nass_active_damping');
#+end_src

Finally, we save the simulation results for further analysis
#+begin_src matlab
  En_iff_hpf = En;
  Eg_iff_hpf = Eg;
  save('./active_damping/mat/tomo_exp.mat', 'En_iff_hpf', 'Eg_iff_hpf', '-append');
#+end_src

*** Compare with Undamped system
We load the results of tomography experiments.
#+begin_src matlab
  load('./active_damping/mat/tomo_exp.mat', 'En', 'En_iff', 'En_iff_hpf');
  t = linspace(0, 3, length(En(:,1)));
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(En(:,1),         En(:,2),         'DisplayName', '$\epsilon_{x,y}$ - OL')
  plot(En_iff(:,1),     En_iff(:,2),     'DisplayName', '$\epsilon_{x,y}$ - IFF')
  plot(En_iff_hpf(:,1), En_iff_hpf(:,2), 'DisplayName', '$\epsilon_{x,y}$ - IFF + HPF')
  xlabel('X Motion [m]'); ylabel('Y Motion [m]');
  legend('location', 'northwest');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_act_damp_iff_sim_tomo_xy.pdf" :var figsize="small-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_act_damp_iff_sim_tomo_xy
#+CAPTION: Position Error during tomography experiment - XY Motion ([[./figs/nass_act_damp_iff_sim_tomo_xy.png][png]], [[./figs/nass_act_damp_iff_sim_tomo_xy.pdf][pdf]])
[[file:figs/nass_act_damp_iff_sim_tomo_xy.png]]

#+begin_src matlab :exports none
  figure;
  ax1 = subplot(3, 1, 1);
  hold on;
  plot(t, En(:,1), 'DisplayName', '$\epsilon_{x}$')
  plot(t, En_iff(:,1), 'DisplayName', '$\epsilon_{x}$ - IFF')
  plot(t, En_iff_hpf(:,1), 'DisplayName', '$\epsilon_{x}$ - IFF + HPF')
  legend('location', 'southwest');

  ax2 = subplot(3, 1, 2);
  hold on;
  plot(t, En(:,2), 'DisplayName', '$\epsilon_{y}$')
  plot(t, En_iff(:,2), 'DisplayName', '$\epsilon_{y}$ - IFF')
  plot(t, En_iff_hpf(:,2), 'DisplayName', '$\epsilon_{y}$ - IFF + HPF')
  legend('location', 'southwest');
  ylabel('Position Error [m]');

  ax3 = subplot(3, 1, 3);
  hold on;
  plot(t, En(:,3), 'DisplayName', '$\epsilon_{z}$')
  plot(t, En_iff(:,3), 'DisplayName', '$\epsilon_{z}$ - IFF')
  plot(t, En_iff_hpf(:,3), 'DisplayName', '$\epsilon_{z}$ - IFF + HPF')
  legend('location', 'northwest');
  xlabel('Time [s]');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_act_damp_iff_sim_tomo_trans.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_act_damp_iff_sim_tomo_trans
#+CAPTION: Position Error during tomography experiment - Translations ([[./figs/nass_act_damp_iff_sim_tomo_trans.png][png]], [[./figs/nass_act_damp_iff_sim_tomo_trans.pdf][pdf]])
[[file:figs/nass_act_damp_iff_sim_tomo_trans.png]]

#+begin_src matlab :exports none
  figure;
  ax1 = subplot(3, 1, 1);
  hold on;
  plot(t, En(:,4), 'DisplayName', '$\epsilon_{\theta_x}$')
  plot(t, En_iff(:,4), 'DisplayName', '$\epsilon_{\theta_x}$ - IFF')
  plot(t, En_iff_hpf(:,4), 'DisplayName', '$\epsilon_{\theta_x}$ - IFF + HPF')
  legend('location', 'northwest');

  ax2 = subplot(3, 1, 2);
  hold on;
  plot(t, En(:,5), 'DisplayName', '$\epsilon_{\theta_y}$')
  plot(t, En_iff(:,5), 'DisplayName', '$\epsilon_{\theta_y}$ - IFF')
  plot(t, En_iff_hpf(:,5), 'DisplayName', '$\epsilon_{\theta_y}$ - IFF + HPF')
  legend('location', 'southwest');
  ylabel('Position Error [rad]');

  ax3 = subplot(3, 1, 3);
  hold on;
  plot(t, En(:,6), 'DisplayName', '$\epsilon_{\theta_z}$')
  plot(t, En_iff(:,6), 'DisplayName', '$\epsilon_{\theta_z}$ - IFF')
  plot(t, En_iff_hpf(:,6), 'DisplayName', '$\epsilon_{\theta_z}$ - IFF + HPF')
  legend();
  xlabel('Time [s]');

  linkaxes([ax1,ax2,ax3],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_act_damp_iff_sim_tomo_rot.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_act_damp_iff_sim_tomo_rot
#+CAPTION: Position Error during tomography experiment - Rotations ([[./figs/nass_act_damp_iff_sim_tomo_rot.png][png]], [[./figs/nass_act_damp_iff_sim_tomo_rot.pdf][pdf]])
[[file:figs/nass_act_damp_iff_sim_tomo_rot.png]]

** Conclusion
#+begin_important
Integral Force Feedback:
- Robust (guaranteed stability)
- Acceptable Damping
- Increase the sensitivity to disturbances at low frequencies
#+end_important

* Direct Velocity Feedback
:PROPERTIES:
:header-args:matlab+: :tangle matlab/dvf.m
:header-args:matlab+: :comments none :mkdirp yes
:END:
<<sec:dvf>>

** ZIP file containing the data and matlab files                     :ignore:
#+begin_src bash :exports none :results none
  if [ matlab/dvf.m -nt data/dvf.zip ]; then
    cp matlab/dvf.m dvf.m;
    zip data/dvf \
        mat/plant.mat \
        dvf.m
    rm dvf.m;
  fi
#+end_src

#+begin_note
  All the files (data and Matlab scripts) are accessible [[file:data/dvf.zip][here]].
#+end_note

** Introduction                                                      :ignore:
In the Direct Velocity Feedback (DVF), a derivative feedback is applied between the measured actuator displacement to the actuator force input.
The actuator displacement can be measured with a capacitive sensor for instance.

** Matlab Init                                              :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
  <<matlab-dir>>
#+end_src

#+begin_src matlab :exports none :results silent :noweb yes
  <<matlab-init>>
#+end_src

#+begin_src matlab :tangle no
  simulinkproject('../');
#+end_src

#+begin_src matlab
  addpath('active_damping/src/');
#+end_src

#+begin_src matlab
  open('active_damping/matlab/sim_nass_active_damping.slx')
#+end_src

** Control Design
*** Plant
Let's load the undamped plant:
#+begin_src matlab
  load('./active_damping/mat/undamped_plants.mat', 'G_dvf');
#+end_src

Let's look at the transfer function from actuator forces in the nano-hexapod to the measured displacement of the actuator for all 6 pairs of actuator/sensor (figure [[fig:dvf_plant]]).

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i=1:6
    plot(freqs, abs(squeeze(freqresp(G_dvf(['Dnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i=1:6
    plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf(['Dnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/dvf_plant.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:dvf_plant
#+CAPTION: Transfer function from forces applied in the legs to leg displacement sensor ([[./figs/dvf_plant.png][png]], [[./figs/dvf_plant.pdf][pdf]])
[[file:figs/dvf_plant.png]]

*** Control Design
The Direct Velocity Feedback is defined below.
A Low pass Filter is added to make the controller transfer function proper.
#+begin_src matlab
  K_dvf = s*20000/(1 + s/2/pi/10000);
#+end_src

The obtained loop gains are shown in figure [[fig:dvf_open_loop]].

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i=1:6
    plot(freqs, abs(squeeze(freqresp(K_dvf*G_dvf(['Dnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i=1:6
    plot(freqs, 180/pi*angle(squeeze(freqresp(K_dvf*G_dvf(['Dnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/dvf_open_loop.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:dvf_open_loop
#+CAPTION: Loop Gain for the Integral Force Feedback ([[./figs/dvf_open_loop.png][png]], [[./figs/dvf_open_loop.pdf][pdf]])
[[file:figs/dvf_open_loop.png]]

*** Diagonal Controller
We create the diagonal controller and we add a minus sign as we have a positive feedback architecture.
#+begin_src matlab
  K_dvf = -K_dvf*eye(6);
#+end_src

We save the controller for further analysis.
#+begin_src matlab
  save('./active_damping/mat/K_dvf.mat', 'K_dvf');
#+end_src

** TODO Identification of the damped plant                         :noexport:
*** Initialize the Simulation
Let's initialize the system prior to identification.
#+begin_src matlab
  initializeReferences();
  initializeGround();
  initializeGranite();
  initializeTy();
  initializeRy();
  initializeRz();
  initializeMicroHexapod();
  initializeAxisc();
  initializeMirror();
  initializeNanoHexapod('actuator', 'piezo');
  initializeSample('mass', 50);
#+end_src

And initialize the controllers.
#+begin_src matlab
  K = tf(zeros(6));
  K_ine = tf(zeros(6));
  save('./mat/controllers.mat', 'K_ine', '-append');
  save('./mat/controllers.mat', 'K', '-append');
  K_iff = tf(zeros(6));
  save('./mat/controllers.mat', 'K_iff', '-append');
  K_dvf = K_dvf;
  save('./mat/controllers.mat', 'K_dvf', '-append');
#+end_src

*** Identification
We identify the system dynamics now that the DVF controller is ON.
#+begin_src matlab
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'sim_nass_active_damping';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fnl'], 1, 'openinput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Dnlm'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Fnlm'); io_i = io_i + 1;

  %% Run the linearization
  G = linearize(mdl, io, options);
  G.InputName  = {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'};
  G.OutputName = {'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6', ...
                  'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6'};
#+end_src

And we save the damped plant for further analysis.
#+begin_src matlab
  save('./active_damping/mat/plants.mat', 'G_dvf', '-append');
#+end_src

*** Sensitivity to disturbances
As shown in figure [[fig:sensitivity_dist_dvf]], DVF control succeed in lowering the sensitivity to disturbances near resonance of the system.

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  subplot(2, 1, 1);
  title('$D_g$ to $D$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_gm('Dx', 'Dgx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / D_{g,x}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_gm('Dy', 'Dgy'), freqs, 'Hz'))), 'DisplayName', '$\left|D_y / D_{g,y}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_gm('Dz', 'Dgz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / D_{g,z}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_gm('Dx', 'Dgx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_gm('Dy', 'Dgy'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_gm('Dz', 'Dgz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]'); xlabel('Frequency [Hz]');
  legend('location', 'southeast');

  subplot(2, 1, 2);
  title('$F_s$ to $D$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_fs('Dx', 'Fsx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / F_{s,x}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_fs('Dy', 'Fsy'), freqs, 'Hz'))), 'DisplayName', '$\left|D_y / F_{s,y}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_fs('Dz', 'Fsz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{s,z}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_fs('Dx', 'Fsx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_fs('Dy', 'Fsy'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_fs('Dz', 'Fsz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_dist_dvf.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_dist_dvf
#+CAPTION: Sensitivity to disturbance once the DVF controller is applied to the system ([[./figs/sensitivity_dist_dvf.png][png]], [[./figs/sensitivity_dist_dvf.pdf][pdf]])
[[file:figs/sensitivity_dist_dvf.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_dist('Dz', 'Frzz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{rz, z}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_dist('Dz', 'Ftyz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{ty, z}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_dist('Dx', 'Ftyx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / F_{ty, x}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_dist('Dz', 'Frzz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_dist('Dz', 'Ftyz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_dist('Dx', 'Ftyx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_dist_stages_dvf.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_dist_stages_dvf
#+CAPTION: Sensitivity to force disturbances in various stages when DVF is applied ([[./figs/sensitivity_dist_stages_dvf.png][png]], [[./figs/sensitivity_dist_stages_dvf.pdf][pdf]])
[[file:figs/sensitivity_dist_stages_dvf.png]]

*** Damped Plant
#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 2, 1);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Dx', 'Fnx'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Dy', 'Fny'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Dz', 'Fnz'), freqs, 'Hz'))));
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), '--');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_cart('Dy', 'Fny'), freqs, 'Hz'))), '--');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), '--');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');

  ax2 = subplot(2, 2, 2);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Rx', 'Mnx'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Ry', 'Mny'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Rz', 'Mnz'), freqs, 'Hz'))));
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_cart('Rx', 'Mnx'), freqs, 'Hz'))), '--');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_cart('Ry', 'Mny'), freqs, 'Hz'))), '--');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_cart('Rz', 'Mnz'), freqs, 'Hz'))), '--');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [rad/(Nm)]'); xlabel('Frequency [Hz]');

  ax3 = subplot(2, 2, 3);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / F_{n,x}\right|$');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Dy', 'Fny'), freqs, 'Hz'))), 'DisplayName', '$\left|D_y / F_{n,y}\right|$');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{n,z}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf.G_cart('Dy', 'Fny'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'northwest');

  ax4 = subplot(2, 2, 4);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Rx', 'Mnx'), freqs, 'Hz'))), 'DisplayName', '$\left|R_x / M_{n,x}\right|$');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Ry', 'Mny'), freqs, 'Hz'))), 'DisplayName', '$\left|R_y / M_{n,y}\right|$');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Rz', 'Mnz'), freqs, 'Hz'))), 'DisplayName', '$\left|R_z / M_{n,z}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf.G_cart('Rx', 'Mnx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf.G_cart('Ry', 'Mny'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf.G_cart('Rz', 'Mnz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'northwest');

  linkaxes([ax1,ax2,ax3,ax4],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/plant_dvf_damped.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:plant_dvf_damped
#+CAPTION: Damped Plant after DVF is applied ([[./figs/plant_dvf_damped.png][png]], [[./figs/plant_dvf_damped.pdf][pdf]])
[[file:figs/plant_dvf_damped.png]]

** Tomography Experiment
*** Initialize the Simulation
We initialize elements for the tomography experiment.
#+begin_src matlab
  prepareTomographyExperiment();
#+end_src

We set the DVF controller.
#+begin_src matlab
  load('./active_damping/mat/K_dvf.mat', 'K_dvf');
  save('./mat/controllers.mat', 'K_dvf', '-append');
#+end_src

*** Simulation
We change the simulation stop time.
#+begin_src matlab
  load('mat/conf_simscape.mat');
  set_param(conf_simscape, 'StopTime', '3');
#+end_src

And we simulate the system.
#+begin_src matlab
  sim('sim_nass_active_damping');
#+end_src

Finally, we save the simulation results for further analysis
#+begin_src matlab
  En_dvf = En;
  Eg_dvf = Eg;
  save('./active_damping/mat/tomo_exp.mat', 'En_dvf', 'Eg_dvf', '-append');
#+end_src

*** Compare with Undamped system
We load the results of tomography experiments.
#+begin_src matlab
  load('./active_damping/mat/tomo_exp.mat', 'En', 'En_dvf');
  t = linspace(0, 3, length(En(:,1)));
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(En(:,1),     En(:,2),     'DisplayName', '$\epsilon_{x,y}$ - OL')
  plot(En_dvf(:,1), En_dvf(:,2), 'DisplayName', '$\epsilon_{x,y}$ - DVF')
  xlabel('X Motion [m]'); ylabel('Y Motion [m]');
  legend();
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_act_damp_dvf_sim_tomo_xy.pdf" :var figsize="small-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_act_damp_dvf_sim_tomo_xy
#+CAPTION: Position Error during tomography experiment - XY Motion ([[./figs/nass_act_damp_dvf_sim_tomo_xy.png][png]], [[./figs/nass_act_damp_dvf_sim_tomo_xy.pdf][pdf]])
[[file:figs/nass_act_damp_dvf_sim_tomo_xy.png]]

#+begin_src matlab :exports none
  figure;
  ax1 = subplot(3, 1, 1);
  hold on;
  plot(t, En(:,1), 'DisplayName', '$\epsilon_{x}$')
  plot(t, En_dvf(:,1), 'DisplayName', '$\epsilon_{x}$ - DVF')
  legend();

  ax2 = subplot(3, 1, 2);
  hold on;
  plot(t, En(:,2), 'DisplayName', '$\epsilon_{y}$')
  plot(t, En_dvf(:,2), 'DisplayName', '$\epsilon_{y}$ - DVF')
  legend();
  ylabel('Position Error [m]');

  ax3 = subplot(3, 1, 3);
  hold on;
  plot(t, En(:,3), 'DisplayName', '$\epsilon_{z}$')
  plot(t, En_dvf(:,3), 'DisplayName', '$\epsilon_{z}$ - DVF')
  legend();
  xlabel('Time [s]');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_act_damp_dvf_sim_tomo_trans.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_act_damp_dvf_sim_tomo_trans
#+CAPTION: Position Error during tomography experiment - Translations ([[./figs/nass_act_damp_dvf_sim_tomo_trans.png][png]], [[./figs/nass_act_damp_dvf_sim_tomo_trans.pdf][pdf]])
[[file:figs/nass_act_damp_dvf_sim_tomo_trans.png]]

#+begin_src matlab :exports none
  figure;
  ax1 = subplot(3, 1, 1);
  hold on;
  plot(t, En(:,4), 'DisplayName', '$\epsilon_{\theta_x}$')
  plot(t, En_dvf(:,4), 'DisplayName', '$\epsilon_{\theta_x}$ - DVF')
  legend();

  ax2 = subplot(3, 1, 2);
  hold on;
  plot(t, En(:,5), 'DisplayName', '$\epsilon_{\theta_y}$')
  plot(t, En_dvf(:,5), 'DisplayName', '$\epsilon_{\theta_y}$ - DVF')
  legend();
  ylabel('Position Error [rad]');

  ax3 = subplot(3, 1, 3);
  hold on;
  plot(t, En(:,6), 'DisplayName', '$\epsilon_{\theta_z}$')
  plot(t, En_dvf(:,6), 'DisplayName', '$\epsilon_{\theta_z}$ - DVF')
  legend();
  xlabel('Time [s]');

  linkaxes([ax1,ax2,ax3],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_act_damp_dvf_sim_tomo_rot.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_act_damp_dvf_sim_tomo_rot
#+CAPTION: Position Error during tomography experiment - Rotations ([[./figs/nass_act_damp_dvf_sim_tomo_rot.png][png]], [[./figs/nass_act_damp_dvf_sim_tomo_rot.pdf][pdf]])
[[file:figs/nass_act_damp_dvf_sim_tomo_rot.png]]

** Conclusion
#+begin_important
Direct Velocity Feedback:
-
#+end_important

* Inertial Control
:PROPERTIES:
:header-args:matlab+: :tangle matlab/ine.m
:header-args:matlab+: :comments none :mkdirp yes
:END:
<<sec:ine>>

** ZIP file containing the data and matlab files                     :ignore:
#+begin_src bash :exports none :results none
  if [ matlab/ine.m -nt data/ine.zip ]; then
    cp matlab/ine.m ine.m;
    zip data/ine \
        mat/plant.mat \
        ine.m
    rm ine.m;
  fi
#+end_src

#+begin_note
  All the files (data and Matlab scripts) are accessible [[file:data/ine.zip][here]].
#+end_note

** Introduction                                                      :ignore:
In Inertial Control, a feedback is applied between the measured *absolute* motion (velocity or acceleration) of the platform to the actuator force input.

** Matlab Init                                              :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
  <<matlab-dir>>
#+end_src

#+begin_src matlab :exports none :results silent :noweb yes
  <<matlab-init>>
#+end_src

#+begin_src matlab :tangle no
  simulinkproject('../');
#+end_src

#+begin_src matlab
  addpath('active_damping/src/');
#+end_src

#+begin_src matlab
  open('active_damping/matlab/sim_nass_active_damping.slx')
#+end_src

** Control Design
*** Plant
Let's load the undamped plant:
#+begin_src matlab
  load('./active_damping/mat/undamped_plants.mat', 'G_ine');
#+end_src

Let's look at the transfer function from actuator forces in the nano-hexapod to the measured velocity of the nano-hexapod platform in the direction of the corresponding actuator for all 6 pairs of actuator/sensor (figure [[fig:ine_plant]]).

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i=1:6
    plot(freqs, abs(squeeze(freqresp(G_ine(['Vnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [$\frac{m/s^2}{N}$]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i=1:6
    plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine(['Vnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/ine_plant.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:ine_plant
#+CAPTION: Transfer function from forces applied in the legs to leg velocity sensor ([[./figs/ine_plant.png][png]], [[./figs/ine_plant.pdf][pdf]])
[[file:figs/ine_plant.png]]

*** Control Design
The controller is defined below and the obtained loop gain is shown in figure [[fig:ine_open_loop_gain]].

#+begin_src matlab
  K_ine = 1e4/(1+s/(2*pi*100));
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  for i=1:6
    plot(freqs, abs(squeeze(freqresp(K_ine*G_ine(['Vnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  for i=1:6
    plot(freqs, 180/pi*angle(squeeze(freqresp(K_ine*G_ine(['Vnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/ine_open_loop_gain.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:ine_open_loop_gain
#+CAPTION: Loop Gain for Inertial Control ([[./figs/ine_open_loop_gain.png][png]], [[./figs/ine_open_loop_gain.pdf][pdf]])
[[file:figs/ine_open_loop_gain.png]]

*** Diagonal Controller
We create the diagonal controller and we add a minus sign as we have a positive feedback architecture.
#+begin_src matlab
  K_ine = -K_ine*eye(6);
#+end_src

We save the controller for further analysis.
#+begin_src matlab
  save('./active_damping/mat/K_ine.mat', 'K_ine');
#+end_src

** TODO Identification of the damped plant                         :noexport:
*** Initialize the Simulation
Let's initialize the system prior to identification.
#+begin_src matlab
  initializeReferences();
  initializeGround();
  initializeGranite();
  initializeTy();
  initializeRy();
  initializeRz();
  initializeMicroHexapod();
  initializeAxisc();
  initializeMirror();
  initializeNanoHexapod('actuator', 'piezo');
  initializeSample('mass', 50);
#+end_src

And initialize the controllers.
#+begin_src matlab
  K = tf(zeros(6));
  save('./mat/controllers.mat', 'K', '-append');
  K_ine = -K_ine*eye(6);
  save('./mat/controllers.mat', 'K_ine', '-append');
  K_iff = tf(zeros(6));
  save('./mat/controllers.mat', 'K_iff', '-append');
  K_dvf = tf(zeros(6));
  save('./mat/controllers.mat', 'K_dvf', '-append');
#+end_src

*** Identification
We identify the system dynamics now that the Inertial controller is ON.
#+begin_src matlab
  G_ine = identifyPlant();
#+end_src

And we save the damped plant for further analysis.
#+begin_src matlab
  save('./active_damping/mat/plants.mat', 'G_ine', '-append');
#+end_src

*** Sensitivity to disturbances

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  subplot(2, 1, 1);
  title('$D_g$ to $D$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_gm('Dx', 'Dgx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / D_{g,x}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_gm('Dy', 'Dgy'), freqs, 'Hz'))), 'DisplayName', '$\left|D_y / D_{g,y}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_gm('Dz', 'Dgz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / D_{g,z}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_ine.G_gm('Dx', 'Dgx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_gm('Dy', 'Dgy'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_gm('Dz', 'Dgz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');

  subplot(2, 1, 2);
  title('$F_s$ to $D$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_fs('Dx', 'Fsx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / F_{s,x}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_fs('Dy', 'Fsy'), freqs, 'Hz'))), 'DisplayName', '$\left|D_y / F_{s,y}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_fs('Dz', 'Fsz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{s,z}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_ine.G_fs('Dx', 'Fsx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_fs('Dy', 'Fsy'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_fs('Dz', 'Fsz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_dist_ine.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_dist_ine
#+CAPTION: Sensitivity to disturbance once the INE controller is applied to the system ([[./figs/sensitivity_dist_ine.png][png]], [[./figs/sensitivity_dist_ine.pdf][pdf]])
[[file:figs/sensitivity_dist_ine.png]]


#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_dist('Dz', 'Frzz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{rz, z}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_dist('Dz', 'Ftyz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{ty, z}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_dist('Dx', 'Ftyx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / F_{ty, x}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_ine.G_dist('Dz', 'Frzz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_dist('Dz', 'Ftyz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_dist('Dx', 'Ftyx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_dist_stages_ine.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_dist_stages_ine
#+CAPTION: Sensitivity to force disturbances in various stages when INE is applied ([[./figs/sensitivity_dist_stages_ine.png][png]], [[./figs/sensitivity_dist_stages_ine.pdf][pdf]])
[[file:figs/sensitivity_dist_stages_ine.png]]

*** Damped Plant
#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 2, 1);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Dx', 'Fnx'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Dy', 'Fny'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Dz', 'Fnz'), freqs, 'Hz'))));
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_ine.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), '--');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_cart('Dy', 'Fny'), freqs, 'Hz'))), '--');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), '--');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');

  ax2 = subplot(2, 2, 2);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Rx', 'Mnx'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Ry', 'Mny'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Rz', 'Mnz'), freqs, 'Hz'))));
  set(gca,'ColorOrderIndex',1);
  plot(freqs, abs(squeeze(freqresp(G_ine.G_cart('Rx', 'Mnx'), freqs, 'Hz'))), '--');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_cart('Ry', 'Mny'), freqs, 'Hz'))), '--');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_cart('Rz', 'Mnz'), freqs, 'Hz'))), '--');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [rad/(Nm)]'); xlabel('Frequency [Hz]');

  ax3 = subplot(2, 2, 3);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / F_{n,x}\right|$');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Dy', 'Fny'), freqs, 'Hz'))), 'DisplayName', '$\left|D_y / F_{n,y}\right|$');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{n,z}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine.G_cart('Dy', 'Fny'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'northwest');

  ax4 = subplot(2, 2, 4);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Rx', 'Mnx'), freqs, 'Hz'))), 'DisplayName', '$\left|R_x / M_{n,x}\right|$');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Ry', 'Mny'), freqs, 'Hz'))), 'DisplayName', '$\left|R_y / M_{n,y}\right|$');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Rz', 'Mnz'), freqs, 'Hz'))), 'DisplayName', '$\left|R_z / M_{n,z}\right|$');
  set(gca,'ColorOrderIndex',1);
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine.G_cart('Rx', 'Mnx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine.G_cart('Ry', 'Mny'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine.G_cart('Rz', 'Mnz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'northwest');

  linkaxes([ax1,ax2,ax3,ax4],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/plant_ine_damped.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:plant_ine_damped
#+CAPTION: Damped Plant after INE is applied ([[./figs/plant_ine_damped.png][png]], [[./figs/plant_ine_damped.pdf][pdf]])
[[file:figs/plant_ine_damped.png]]

** Tomography Experiment
*** Initialize the Simulation
We initialize elements for the tomography experiment.
#+begin_src matlab
  prepareTomographyExperiment();
#+end_src

We set the Inertial controller.
#+begin_src matlab
  load('./active_damping/mat/K_ine.mat', 'K_ine');
  save('./mat/controllers.mat', 'K_ine', '-append');
#+end_src

*** Simulation
We change the simulation stop time.
#+begin_src matlab
  load('mat/conf_simscape.mat');
  set_param(conf_simscape, 'StopTime', '3');
#+end_src

And we simulate the system.
#+begin_src matlab
  sim('sim_nass_active_damping');
#+end_src

Finally, we save the simulation results for further analysis
#+begin_src matlab
  En_ine = En;
  Eg_ine = Eg;
  save('./active_damping/mat/tomo_exp.mat', 'En_ine', 'Eg_ine', '-append');
#+end_src

*** Compare with Undamped system
We load the results of tomography experiments.
#+begin_src matlab
  load('./active_damping/mat/tomo_exp.mat', 'En', 'En_ine');
  t = linspace(0, 3, length(En_ine(:,1)));
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(En(:,1),     En(:,2),     'DisplayName', '$\epsilon_{x,y}$ - OL')
  plot(En_ine(:,1), En_ine(:,2), 'DisplayName', '$\epsilon_{x,y}$ - Inertial')
  xlabel('X Motion [m]'); ylabel('Y Motion [m]');
  legend();
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_act_damp_ine_sim_tomo_xy.pdf" :var figsize="small-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_act_damp_ine_sim_tomo_xy
#+CAPTION: Position Error during tomography experiment - XY Motion ([[./figs/nass_act_damp_ine_sim_tomo_xy.png][png]], [[./figs/nass_act_damp_ine_sim_tomo_xy.pdf][pdf]])
[[file:figs/nass_act_damp_ine_sim_tomo_xy.png]]

#+begin_src matlab :exports none
  figure;
  ax1 = subplot(3, 1, 1);
  hold on;
  plot(t, En(:,1), 'DisplayName', '$\epsilon_{x}$')
  plot(t, En_ine(:,1), 'DisplayName', '$\epsilon_{x}$ - Inertial')
  legend();

  ax2 = subplot(3, 1, 2);
  hold on;
  plot(t, En(:,2), 'DisplayName', '$\epsilon_{y}$')
  plot(t, En_ine(:,2), 'DisplayName', '$\epsilon_{y}$ - Inertial')
  legend();
  ylabel('Position Error [m]');

  ax3 = subplot(3, 1, 3);
  hold on;
  plot(t, En(:,3), 'DisplayName', '$\epsilon_{z}$')
  plot(t, En_ine(:,3), 'DisplayName', '$\epsilon_{z}$ - Inertial')
  legend();
  xlabel('Time [s]');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_act_damp_ine_sim_tomo_trans.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_act_damp_ine_sim_tomo_trans
#+CAPTION: Position Error during tomography experiment - Translations ([[./figs/nass_act_damp_ine_sim_tomo_trans.png][png]], [[./figs/nass_act_damp_ine_sim_tomo_trans.pdf][pdf]])
[[file:figs/nass_act_damp_ine_sim_tomo_trans.png]]

#+begin_src matlab :exports none
  figure;
  ax1 = subplot(3, 1, 1);
  hold on;
  plot(t, En(:,4), 'DisplayName', '$\epsilon_{\theta_x}$')
  plot(t, En_ine(:,4), 'DisplayName', '$\epsilon_{\theta_x}$ - Inertial')
  legend();

  ax2 = subplot(3, 1, 2);
  hold on;
  plot(t, En(:,5), 'DisplayName', '$\epsilon_{\theta_y}$')
  plot(t, En_ine(:,5), 'DisplayName', '$\epsilon_{\theta_y}$ - Inertial')
  legend();
  ylabel('Position Error [rad]');

  ax3 = subplot(3, 1, 3);
  hold on;
  plot(t, En(:,6), 'DisplayName', '$\epsilon_{\theta_z}$')
  plot(t, En_ine(:,6), 'DisplayName', '$\epsilon_{\theta_z}$ - Inertial')
  legend();
  xlabel('Time [s]');

  linkaxes([ax1,ax2,ax3],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/nass_act_damp_ine_sim_tomo_rot.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:nass_act_damp_ine_sim_tomo_rot
#+CAPTION: Position Error during tomography experiment - Rotations ([[./figs/nass_act_damp_ine_sim_tomo_rot.png][png]], [[./figs/nass_act_damp_ine_sim_tomo_rot.pdf][pdf]])
[[file:figs/nass_act_damp_ine_sim_tomo_rot.png]]

** Conclusion
#+begin_important
Inertial Control:
#+end_important

* Comparison
<<sec:comparison>>
** Introduction                                                     :ignore:
** Matlab Init                                             :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
  <<matlab-dir>>
#+end_src

#+begin_src matlab :exports none :results silent :noweb yes
  <<matlab-init>>
#+end_src

#+begin_src matlab
  cd('../');
#+end_src

** Load the plants
#+begin_src matlab
  load('./active_damping/mat/plants.mat', 'G', 'G_iff', 'G_ine', 'G_dvf');
#+end_src

** Sensitivity to Disturbance
#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;
  title('$D_{g,z}$ to $D_z$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_gm(    'Dz', 'Dgz'), freqs, 'Hz'))), 'k-' , 'DisplayName', 'Undamped');
  plot(freqs, abs(squeeze(freqresp(G_iff.G_gm('Dz', 'Dgz'), freqs, 'Hz'))), 'k:' , 'DisplayName', 'IFF');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_gm('Dz', 'Dgz'), freqs, 'Hz'))), 'k--', 'DisplayName', 'INE');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_gm('Dz', 'Dgz'), freqs, 'Hz'))), 'k-.', 'DisplayName', 'DVF');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_comp_ground_motion_z.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_comp_ground_motion_z
#+CAPTION: Sensitivity to ground motion in the Z direction on the Z motion error ([[./figs/sensitivity_comp_ground_motion_z.png][png]], [[./figs/sensitivity_comp_ground_motion_z.pdf][pdf]])
[[file:figs/sensitivity_comp_ground_motion_z.png]]


#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;
  title('$F_{s,z}$ to $D_z$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_fs(    'Dz', 'Fsz'), freqs, 'Hz'))), 'k-' , 'DisplayName', 'Undamped');
  plot(freqs, abs(squeeze(freqresp(G_iff.G_fs('Dz', 'Fsz'), freqs, 'Hz'))), 'k:' , 'DisplayName', 'IFF');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_fs('Dz', 'Fsz'), freqs, 'Hz'))), 'k--', 'DisplayName', 'INE');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_fs('Dz', 'Fsz'), freqs, 'Hz'))), 'k-.', 'DisplayName', 'DVF');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_comp_direct_forces_z.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_comp_direct_forces_z
#+CAPTION: Compliance in the Z direction: Sensitivity of direct forces applied on the sample in the Z direction on the Z motion error ([[./figs/sensitivity_comp_direct_forces_z.png][png]], [[./figs/sensitivity_comp_direct_forces_z.pdf][pdf]])
[[file:figs/sensitivity_comp_direct_forces_z.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;
  title('$F_{rz,z}$ to $D_z$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_dist(    'Dz', 'Frzz'), freqs, 'Hz'))), 'k-' , 'DisplayName', 'Undamped');
  plot(freqs, abs(squeeze(freqresp(G_iff.G_dist('Dz', 'Frzz'), freqs, 'Hz'))), 'k:' , 'DisplayName', 'IFF');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_dist('Dz', 'Frzz'), freqs, 'Hz'))), 'k--', 'DisplayName', 'INE');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_dist('Dz', 'Frzz'), freqs, 'Hz'))), 'k-.', 'DisplayName', 'DVF');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_comp_spindle_z.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_comp_spindle_z
#+CAPTION: Sensitivity to forces applied in the Z direction by the Spindle on the Z motion error ([[./figs/sensitivity_comp_spindle_z.png][png]], [[./figs/sensitivity_comp_spindle_z.pdf][pdf]])
[[file:figs/sensitivity_comp_spindle_z.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;
  title('$F_{ty,z}$ to $D_z$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_dist(    'Dz', 'Ftyz'), freqs, 'Hz'))), 'k-' , 'DisplayName', 'Undamped');
  plot(freqs, abs(squeeze(freqresp(G_iff.G_dist('Dz', 'Ftyz'), freqs, 'Hz'))), 'k:' , 'DisplayName', 'IFF');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_dist('Dz', 'Ftyz'), freqs, 'Hz'))), 'k--', 'DisplayName', 'INE');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_dist('Dz', 'Ftyz'), freqs, 'Hz'))), 'k-.', 'DisplayName', 'DVF');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_comp_ty_z.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_comp_ty_z
#+CAPTION: Sensitivity to forces applied in the Z direction by the Y translation stage on the Z motion error ([[./figs/sensitivity_comp_ty_z.png][png]], [[./figs/sensitivity_comp_ty_z.pdf][pdf]])
[[file:figs/sensitivity_comp_ty_z.png]]


#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;
  title('$F_{ty,x}$ to $D_x$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_dist(    'Dx', 'Ftyx'), freqs, 'Hz'))), 'k-' , 'DisplayName', 'Undamped');
  plot(freqs, abs(squeeze(freqresp(G_iff.G_dist('Dx', 'Ftyx'), freqs, 'Hz'))), 'k:' , 'DisplayName', 'IFF');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_dist('Dx', 'Ftyx'), freqs, 'Hz'))), 'k--', 'DisplayName', 'INE');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_dist('Dx', 'Ftyx'), freqs, 'Hz'))), 'k-.', 'DisplayName', 'DVF');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_comp_ty_x.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_comp_ty_x
#+CAPTION: Sensitivity to forces applied in the X direction by the Y translation stage on the X motion error ([[./figs/sensitivity_comp_ty_x.png][png]], [[./figs/sensitivity_comp_ty_x.pdf][pdf]])
[[file:figs/sensitivity_comp_ty_x.png]]

** Damped Plant
#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  title('$F_{n,z}$ to $D_z$');
  ax1 = subplot(2, 1, 1);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_cart(    'Dz', 'Fnz'), freqs, 'Hz'))), 'k-' , 'DisplayName', 'Undamped');
  plot(freqs, abs(squeeze(freqresp(G_iff.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), 'k:' , 'DisplayName', 'IFF');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), 'k--', 'DisplayName', 'INE');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), 'k-.', 'DisplayName', 'DVF');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);
  legend('location', 'northeast');

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart    ('Dz', 'Fnz'), freqs, 'Hz'))), 'k-');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), 'k:');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), 'k--');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), 'k-.');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/plant_comp_damping_z.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:plant_comp_damping_z
#+CAPTION: Plant for the $z$ direction for different active damping technique used ([[./figs/plant_comp_damping_z.png][png]], [[./figs/plant_comp_damping_z.pdf][pdf]])
[[file:figs/plant_comp_damping_z.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  title('$F_{n,z}$ to $D_z$');
  ax1 = subplot(2, 1, 1);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_cart(    'Dx', 'Fnx'), freqs, 'Hz'))), 'k-' , 'DisplayName', 'Undamped');
  plot(freqs, abs(squeeze(freqresp(G_iff.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), 'k:' , 'DisplayName', 'IFF');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), 'k--', 'DisplayName', 'INE');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), 'k-.', 'DisplayName', 'DVF');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);
  legend('location', 'northeast');

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart    ('Dx', 'Fnx'), freqs, 'Hz'))), 'k-');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), 'k:');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), 'k--');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), 'k-.');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/plant_comp_damping_x.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:plant_comp_damping_x
#+CAPTION: Plant for the $x$ direction for different active damping technique used ([[./figs/plant_comp_damping_x.png][png]], [[./figs/plant_comp_damping_x.pdf][pdf]])
[[file:figs/plant_comp_damping_x.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  title('$F_{n,x}$ to $R_z$');
  ax1 = subplot(2, 1, 1);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_cart(    'Rz', 'Fnx'), freqs, 'Hz'))), 'k-' , 'DisplayName', 'Undamped');
  plot(freqs, abs(squeeze(freqresp(G_iff.G_cart('Rz', 'Fnx'), freqs, 'Hz'))), 'k:' , 'DisplayName', 'IFF');
  plot(freqs, abs(squeeze(freqresp(G_ine.G_cart('Rz', 'Fnx'), freqs, 'Hz'))), 'k--', 'DisplayName', 'INE');
  plot(freqs, abs(squeeze(freqresp(G_dvf.G_cart('Rz', 'Fnx'), freqs, 'Hz'))), 'k-.', 'DisplayName', 'DVF');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);
  legend('location', 'northeast');

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart    ('Ry', 'Fnx'), freqs, 'Hz'))), 'k-');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff.G_cart('Ry', 'Fnx'), freqs, 'Hz'))), 'k:');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_ine.G_cart('Ry', 'Fnx'), freqs, 'Hz'))), 'k--');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf.G_cart('Ry', 'Fnx'), freqs, 'Hz'))), 'k-.');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/plant_comp_damping_coupling.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:plant_comp_damping_coupling
#+CAPTION: Comparison of one off-diagonal plant for different damping technique applied ([[./figs/plant_comp_damping_coupling.png][png]], [[./figs/plant_comp_damping_coupling.pdf][pdf]])
[[file:figs/plant_comp_damping_coupling.png]]

** Tomography Experiment
*** Load the Simulation Data
#+begin_src matlab
  load('./active_damping/mat/tomo_exp.mat', 'En', 'En_iff_hpf', 'En_dvf', 'En_ine');
  En_iff = En_iff_hpf;
  t = linspace(0, 3, length(En(:,1)));
#+end_src

*** Frequency Domain Analysis
Window used for =pwelch= function.
#+begin_src matlab
  n_av = 8;
  han_win = hanning(ceil(length(En(:, 1))/n_av));
#+end_src

#+begin_src matlab :exports none
  Ts = t(2)-t(1); % Sample Time for the Data [s]

  [pxx,     f] = pwelch(En(:,     1), han_win, [], [], 1/Ts);
  [pxx_ine, ~] = pwelch(En_ine(:, 1), han_win, [], [], 1/Ts);
  [pxx_dvf, ~] = pwelch(En_dvf(:, 1), han_win, [], [], 1/Ts);
  [pxx_iff, ~] = pwelch(En_iff(:, 1), han_win, [], [], 1/Ts);

  [pyy,     ~] = pwelch(En(:,     2), han_win, [], [], 1/Ts);
  [pyy_ine, ~] = pwelch(En_ine(:, 2), han_win, [], [], 1/Ts);
  [pyy_dvf, ~] = pwelch(En_dvf(:, 2), han_win, [], [], 1/Ts);
  [pyy_iff, ~] = pwelch(En_iff(:, 2), han_win, [], [], 1/Ts);

  [pzz,     ~] = pwelch(En(:,     3), han_win, [], [], 1/Ts);
  [pzz_ine, ~] = pwelch(En_ine(:, 3), han_win, [], [], 1/Ts);
  [pzz_dvf, ~] = pwelch(En_dvf(:, 3), han_win, [], [], 1/Ts);
  [pzz_iff, ~] = pwelch(En_iff(:, 3), han_win, [], [], 1/Ts);

  [prx,     ~] = pwelch(En(:,     4), han_win, [], [], 1/Ts);
  [prx_ine, ~] = pwelch(En_ine(:, 4), han_win, [], [], 1/Ts);
  [prx_dvf, ~] = pwelch(En_dvf(:, 4), han_win, [], [], 1/Ts);
  [prx_iff, ~] = pwelch(En_iff(:, 4), han_win, [], [], 1/Ts);

  [pry,     ~] = pwelch(En(:,     5), han_win, [], [], 1/Ts);
  [pry_ine, ~] = pwelch(En_ine(:, 5), han_win, [], [], 1/Ts);
  [pry_dvf, ~] = pwelch(En_dvf(:, 5), han_win, [], [], 1/Ts);
  [pry_iff, ~] = pwelch(En_iff(:, 5), han_win, [], [], 1/Ts);

  [prz,     ~] = pwelch(En(:,     6), han_win, [], [], 1/Ts);
  [prz_ine, ~] = pwelch(En_ine(:, 6), han_win, [], [], 1/Ts);
  [prz_dvf, ~] = pwelch(En_dvf(:, 6), han_win, [], [], 1/Ts);
  [prz_iff, ~] = pwelch(En_iff(:, 6), han_win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f, pxx_ine,    'DisplayName', 'Inertial')
  plot(f, pxx_dvf,    'DisplayName', 'DVF')
  plot(f, pxx_iff,    'DisplayName', 'IFF')
  plot(f, pxx, 'k--', 'DisplayName', 'Undamped')
  hold off;
  xlabel('Frequency [Hz]');
  ylabel('Power Spectral Density [$m^2/Hz$]');
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  legend('location', 'northeast');
  xlim([2, 500]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_tomo_exp_comp_psd_trans.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_tomo_exp_comp_psd_trans
#+CAPTION: PSD of the translation errors in the X direction for applied Active Damping techniques ([[./figs/act_damp_tomo_exp_comp_psd_trans.png][png]], [[./figs/act_damp_tomo_exp_comp_psd_trans.pdf][pdf]])
[[file:figs/act_damp_tomo_exp_comp_psd_trans.png]]

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f, prx_ine,    'DisplayName', 'Inertial')
  plot(f, prx_dvf,    'DisplayName', 'DVF')
  plot(f, prx_iff,    'DisplayName', 'IFF')
  plot(f, prx, 'k--', 'DisplayName', 'Undamped')
  hold off;
  xlabel('Frequency [Hz]');
  ylabel('Power Spectral Density [$rad^2/Hz$]');
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  legend('location', 'northeast');
  xlim([2, 500]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_tomo_exp_comp_psd_rot.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_tomo_exp_comp_psd_rot
#+CAPTION: PSD of the rotation errors in the X direction for applied Active Damping techniques ([[./figs/act_damp_tomo_exp_comp_psd_rot.png][png]], [[./figs/act_damp_tomo_exp_comp_psd_rot.pdf][pdf]])
[[file:figs/act_damp_tomo_exp_comp_psd_rot.png]]

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f, flip(-cumtrapz(flip(f), flip(pxx_ine))),    'DisplayName', 'Inertial')
  plot(f, flip(-cumtrapz(flip(f), flip(pxx_dvf))),    'DisplayName', 'DVF')
  plot(f, flip(-cumtrapz(flip(f), flip(pxx_iff))),    'DisplayName', 'IFF')
  plot(f, flip(-cumtrapz(flip(f), flip(pxx))), 'k--', 'DisplayName', 'Undamped')
  hold off;
  xlabel('Frequency [Hz]');
  ylabel('Cumulative Power Spectrum [$m^2$]');
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  legend('location', 'northeast');
  xlim([2, 500]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_tomo_exp_comp_cps_trans.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_tomo_exp_comp_cps_trans
#+CAPTION: CPS of the translation errors in the X direction for applied Active Damping techniques ([[./figs/act_damp_tomo_exp_comp_cps_trans.png][png]], [[./figs/act_damp_tomo_exp_comp_cps_trans.pdf][pdf]])
[[file:figs/act_damp_tomo_exp_comp_cps_trans.png]]

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f, flip(-cumtrapz(flip(f), flip(prx_ine))),    'DisplayName', 'Inertial')
  plot(f, flip(-cumtrapz(flip(f), flip(prx_dvf))),    'DisplayName', 'DVF')
  plot(f, flip(-cumtrapz(flip(f), flip(prx_iff))),    'DisplayName', 'IFF')
  plot(f, flip(-cumtrapz(flip(f), flip(prx))), 'k--', 'DisplayName', 'Undamped')
  hold off;
  xlabel('Frequency [Hz]');
  ylabel('Cumulative Power Spectrum [$rad^2$]');
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  legend('location', 'northeast');
  xlim([2, 500]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/act_damp_tomo_exp_comp_cps_rot.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+NAME: fig:act_damp_tomo_exp_comp_cps_rot
#+CAPTION: CPS of the rotation errors in the X direction for applied Active Damping techniques ([[./figs/act_damp_tomo_exp_comp_cps_rot.png][png]], [[./figs/act_damp_tomo_exp_comp_cps_rot.pdf][pdf]])
[[file:figs/act_damp_tomo_exp_comp_cps_rot.png]]

* Useful Functions
** prepareTomographyExperiment
:PROPERTIES:
:header-args:matlab+: :tangle src/prepareTomographyExperiment.m
:header-args:matlab+: :comments none :mkdirp yes :eval no
:END:
<<sec:prepareTomographyExperiment>>

This Matlab function is accessible [[file:src/prepareTomographyExperiment.m][here]].

*** Function Description
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
  function [] = prepareTomographyExperiment(args)
#+end_src

*** Optional Parameters
:PROPERTIES:
:UNNUMBERED: t
:END:
#+begin_src matlab
  arguments
      args.nass_actuator       char   {mustBeMember(args.nass_actuator,{'piezo', 'lorentz'})} = 'piezo'
      args.sample_mass   (1,1) double {mustBeNumeric, mustBePositive} = 50 % [kg]
      args.Rz_period     (1,1) double {mustBeNumeric, mustBePositive} = 1 % [s]
  end
#+end_src

*** Initialize the Simulation
:PROPERTIES:
:UNNUMBERED: t
:END:
We initialize all the stages with the default parameters.
#+begin_src matlab
  initializeGround();
  initializeGranite();
  initializeTy();
  initializeRy();
  initializeRz();
  initializeMicroHexapod();
  initializeAxisc();
  initializeMirror();
#+end_src

The nano-hexapod is a piezoelectric hexapod and the sample has a mass of 50kg.
#+begin_src matlab
  initializeNanoHexapod('actuator', args.nass_actuator);
  initializeSample('mass', args.sample_mass);
#+end_src

We set the references that corresponds to a tomography experiment.
#+begin_src matlab
  initializeReferences('Rz_type', 'rotating', 'Rz_period', args.Rz_period);
#+end_src

#+begin_src matlab
  initializeDisturbances();
#+end_src

And all the controllers are set to 0.
#+begin_src matlab
  K = tf(zeros(6));
  save('./mat/controllers.mat', 'K', '-append');
  K_ine = tf(zeros(6));
  save('./mat/controllers.mat', 'K_ine', '-append');
  K_iff = tf(zeros(6));
  save('./mat/controllers.mat', 'K_iff', '-append');
  K_dvf = tf(zeros(6));
  save('./mat/controllers.mat', 'K_dvf', '-append');
#+end_src

* TODO Order                                                        :noexport:
** Undamped
*** Identification of the transfer function from disturbance to position error
#+begin_src matlab
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'sim_nass_active_damping';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Disturbances'], 1, 'openinput', [], 'Dwx'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Disturbances'], 1, 'openinput', [], 'Dwy'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Disturbances'], 1, 'openinput', [], 'Dwz'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Compute Error in NASS base'], 2, 'openoutput'); io_i = io_i + 1;

  %% Run the linearization
  G = linearize(mdl, io, options);
  G.InputName  = {'Dwx', 'Dwy', 'Dwz'};
  G.OutputName = {'Edx', 'Edy', 'Edz', 'Erx', 'Ery', 'Erz'};
#+end_src

*** Identification of the plant
#+begin_src matlab
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'sim_nass_active_damping';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fnl'], 1, 'openinput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Compute Error in NASS base'], 2, 'openoutput'); io_i = io_i + 1;

  %% Run the linearization
  G = linearize(mdl, io, options);
  G.InputName  = {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'};
  G.OutputName = {'Edx', 'Edy', 'Edz', 'Erx', 'Ery', 'Erz'};
#+end_src

#+begin_src matlab
  load('mat/stages.mat', 'nano_hexapod');
  G_cart = G*inv(nano_hexapod.J');
  G_cart.InputName = {'Fnx', 'Fny', 'Fnz', 'Mnx', 'Mny', 'Mnz'};
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_cart('Edx', 'Fnx'), freqs, 'Hz'))), 'DisplayName', '$T_{x}$');
  plot(freqs, abs(squeeze(freqresp(G_cart('Edy', 'Fny'), freqs, 'Hz'))), 'DisplayName', '$T_{y}$');
  plot(freqs, abs(squeeze(freqresp(G_cart('Edz', 'Fnz'), freqs, 'Hz'))), 'DisplayName', '$T_{z}$');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);
  legend('location', 'southwest')

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_cart('Edx', 'Fnx'), freqs, 'Hz'))));
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_cart('Edy', 'Fny'), freqs, 'Hz'))));
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_cart('Edz', 'Fnz'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_cart('Erx', 'Mnx'), freqs, 'Hz'))), 'DisplayName', '$R_{x}$');
  plot(freqs, abs(squeeze(freqresp(G_cart('Ery', 'Mny'), freqs, 'Hz'))), 'DisplayName', '$R_{y}$');
  plot(freqs, abs(squeeze(freqresp(G_cart('Erz', 'Mnz'), freqs, 'Hz'))), 'DisplayName', '$R_{z}$');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);
  legend('location', 'southwest')

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_cart('Erx', 'Mnx'), freqs, 'Hz'))));
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_cart('Ery', 'Mny'), freqs, 'Hz'))));
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_cart('Erz', 'Mnz'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

*** TODO test
#+begin_src matlab
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'sim_nass_active_damping';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Micro-Station/Dy'], 1, 'openinput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Compute Error in NASS base'], 2, 'openoutput'); io_i = io_i + 1;

  %% Run the linearization
  G = linearize(mdl, io, options);
  G.InputName  = {'Dy'};
  G.OutputName = {'Edx', 'Edy', 'Edz', 'Erx', 'Ery', 'Erz'};
#+end_src

#+begin_important
Why is the transfer function from Ty displacement to position error is equal to
1 at all frequencies?
Why don't we see any resonance?
#+end_important

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G('Edy', 'Dy(1)'), freqs, 'Hz'))), 'DisplayName', '$T_{x}$');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);
  legend('location', 'southwest')

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G('Edy', 'Dy(1)'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

*** TODO test on hexapod
#+begin_src matlab
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'test_nano_hexapod';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fnl'], 1, 'openinput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/x'],  1, 'openoutput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/y'],  1, 'openoutput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/z'],  1, 'openoutput'); io_i = io_i + 1;

  %% Run the linearization
  G = linearize(mdl, io, options);
  G.InputName  = {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'};
  G.OutputName = {'x', 'y', 'z'};
#+end_src

#+begin_src matlab
  %% Options for Linearized
  options = linearizeOptions;
  options.SampleTime = 0;

  %% Name of the Simulink File
  mdl = 'test_nano_hexapod';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fx'], 1, 'openinput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/x'], 1, 'openoutput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/z'], 1, 'openoutput'); io_i = io_i + 1;

  %% Run the linearization
  G = linearize(mdl, io, options);
  G.InputName  = {'Fx'};
  G.OutputName = {'x', 'y', 'z'};
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G('Edy', 'Dy(1)'), freqs, 'Hz'))), 'DisplayName', '$T_{x}$');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);
  legend('location', 'southwest')

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G('Edy', 'Dy(1)'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_cart('Erx', 'Mnx'), freqs, 'Hz'))), 'DisplayName', '$R_{x}$');
  plot(freqs, abs(squeeze(freqresp(G_cart('Ery', 'Mny'), freqs, 'Hz'))), 'DisplayName', '$R_{y}$');
  plot(freqs, abs(squeeze(freqresp(G_cart('Erz', 'Mnz'), freqs, 'Hz'))), 'DisplayName', '$R_{z}$');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);
  legend('location', 'southwest')

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_cart('Erx', 'Mnx'), freqs, 'Hz'))));
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_cart('Ery', 'Mny'), freqs, 'Hz'))));
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_cart('Erz', 'Mnz'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);

  linkaxes([ax1,ax2],'x');
#+end_src

*** Sensitivity to disturbances
The sensitivity to disturbances are shown on figure [[fig:sensitivity_dist_undamped]].

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  subplot(2, 1, 1);
  title('$D_g$ to $D$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_gm('Dx', 'Dgx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / D_{g,x}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_gm('Dy', 'Dgy'), freqs, 'Hz'))), 'DisplayName', '$\left|D_y / D_{g,y}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_gm('Dz', 'Dgz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / D_{g,z}\right|$');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]'); xlabel('Frequency [Hz]');
  legend('location', 'southeast');

  subplot(2, 1, 2);
  title('$F_s$ to $D$');
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_fs('Dx', 'Fsx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / F_{s,x}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_fs('Dy', 'Fsy'), freqs, 'Hz'))), 'DisplayName', '$\left|D_y / F_{s,y}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_fs('Dz', 'Fsz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{s,z}\right|$');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_dist_undamped.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_dist_undamped
#+CAPTION: Undamped sensitivity to disturbances ([[./figs/sensitivity_dist_undamped.png][png]], [[./figs/sensitivity_dist_undamped.pdf][pdf]])
[[file:figs/sensitivity_dist_undamped.png]]

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_dist('Dz', 'Frzz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{rz, z}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_dist('Dz', 'Ftyz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{ty, z}\right|$');
  plot(freqs, abs(squeeze(freqresp(G.G_dist('Dx', 'Ftyx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / F_{ty, x}\right|$');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_dist_stages.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:sensitivity_dist_stages
#+CAPTION: Sensitivity to force disturbances in various stages ([[./figs/sensitivity_dist_stages.png][png]], [[./figs/sensitivity_dist_stages.pdf][pdf]])
[[file:figs/sensitivity_dist_stages.png]]
*** Undamped Plant
The "plant" (transfer function from forces applied by the nano-hexapod to the measured displacement of the sample with respect to the granite) bode plot is shown on figure [[fig:sensitivity_dist_undamped]].

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;

  ax1 = subplot(2, 1, 1);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Dx', 'Fnx'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Dy', 'Fny'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(G.G_cart('Dz', 'Fnz'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), 'DisplayName', '$D_x / F_x$');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Dy', 'Fny'), freqs, 'Hz'))), 'DisplayName', '$D_y / F_y$');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), 'DisplayName', '$D_z / F_z$');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
  ylim([-180, 180]);
  yticks([-180, -90, 0, 90, 180]);
  legend('location', 'southwest');

  linkaxes([ax1,ax2],'x');
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/plant_undamped.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:plant_undamped
#+CAPTION: Transfer Function from cartesian forces to displacement for the undamped plant ([[./figs/plant_undamped.png][png]], [[./figs/plant_undamped.pdf][pdf]])
[[file:figs/plant_undamped.png]]


** Direct Velocity Feedback
*** One degree-of-freedom example
:PROPERTIES:
:header-args:matlab+: :tangle no
:END:
<<sec:rmc_1dof>>
**** Equations
#+begin_src latex :file rmc_1dof.pdf :post pdf2svg(file=*this*, ext="png") :exports results
  \begin{tikzpicture}
    % Ground
    \draw (-1, 0) -- (1, 0);

    % Ground Displacement
    \draw[dashed] (-1, 0) -- ++(-0.5, 0) coordinate(w);
    \draw[->] (w) -- ++(0, 0.5) node[left]{$w$};

    % Mass
    \draw[fill=white] (-1, 1) rectangle ++(2, 0.8) node[pos=0.5]{$m$};

    % Displacement of the mass
    \draw[dashed] (-1, 1.8) -- ++(-0.5, 0) coordinate(x);
    \draw[->] (x) -- ++(0, 0.5) node[left]{$x$};

    % Spring, Damper, and Actuator
    \draw[spring] (-0.8, 0) -- (-0.8, 1) node[midway, left=0.1]{$k$};
    \draw[damper] (0, 0)    -- (0, 1)    node[midway, left=0.2]{$c$};
    \draw[actuator={0.4}{0.2}] (0.8, 0) -- (0.8, 1)  coordinate[midway, right=0.1](F);

    % Measured deformation
    \draw[dashed] (1, 0) -- ++(2, 0) coordinate(d_bot);
    \draw[dashed] (1, 1) -- ++(2, 0) coordinate(d_top);
    \draw[<->] (d_bot) --coordinate[midway](d) (d_top);

    % Displacements
    \node[block={0.8cm}{0.6cm}, right=0.6 of F] (Krmc) {$K$};
    \draw[->] (Krmc.west) -- (F) node[above right]{$F$};
    \draw[->] (d)node[above left]{$d$} -- (Krmc.east);
  \end{tikzpicture}
#+end_src

#+name: fig:rmc_1dof
#+caption: Relative Motion Control applied to a 1dof system
#+RESULTS:
[[file:figs/rmc_1dof.png]]

The dynamic of the system is:
\begin{equation}
  ms^2x = F_d - kx - csx + kw + csw + F
\end{equation}
In terms of the stage deformation $d = x - w$:
\begin{equation}
  (ms^2 + cs + k) d = -ms^2 w + F_d + F
\end{equation}
The relative motion control law is:
\begin{equation}
  K = -g s
\end{equation}
Thus, the applied force is:
\begin{equation}
  F = -g s d
\end{equation}
And the new dynamics will be:
\begin{equation}
  d = w \frac{-ms^2}{ms^2 + (c + g)s + k} + F_d \frac{1}{ms^2 + (c + g)s + k} + F \frac{1}{ms^2 + (c + g)s + k}
\end{equation}

And thus damping is added.

If critical damping is wanted:
\begin{equation}
  \xi = \frac{1}{2}\frac{c + g}{\sqrt{km}} = \frac{1}{2}
\end{equation}
This corresponds to a gain:
\begin{equation}
  g = \sqrt{km} - c
\end{equation}

**** Matlab Init                                          :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
  <<matlab-dir>>
#+end_src

#+begin_src matlab :exports none :results silent :noweb yes
  <<matlab-init>>
#+end_src

**** Matlab Example
Let define the system parameters.
#+begin_src matlab
  m = 50; % [kg]
  k = 1e6; % [N/m]
  c = 1e3; % [N/(m/s)]
#+end_src

The state space model of the system is defined below.
#+begin_src matlab
  A = [-c/m -k/m;
       1     0];

  B = [1/m 1/m -1;
       0   0    0];

  C = [ 0  1;
       -c -k];

  D = [0 0 0;
       1 0 0];

  sys = ss(A, B, C, D);
  sys.InputName = {'F', 'Fd', 'wddot'};
  sys.OutputName = {'d', 'Fm'};
  sys.StateName = {'ddot', 'd'};
#+end_src

The controller $K_\text{RMC}$ is:
#+begin_src matlab
  Krmc = -(sqrt(k*m)-c)*s;
  Krmc.InputName = {'d'};
  Krmc.OutputName = {'F'};
#+end_src

And the closed loop system is computed below.
#+begin_src matlab
  sys_rmc = feedback(sys, Krmc, 'name', +1);
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(-1, 3, 1000);

  figure;

  subplot(2, 2, 1);
  title('Fd to d')
  hold on;
  plot(freqs, abs(squeeze(freqresp(sys('d', 'Fd'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(sys_rmc('d', 'Fd'), freqs, 'Hz'))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  xlim([freqs(1), freqs(end)]);

  subplot(2, 2, 3);
  title('Fd to x')
  hold on;
  plot(freqs, abs(squeeze(freqresp(sys('d', 'Fd'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(sys_rmc('d', 'Fd'), freqs, 'Hz'))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  xlim([freqs(1), freqs(end)]);

  subplot(2, 2, 2);
  title('w to d')
  hold on;
  plot(freqs, abs(squeeze(freqresp(sys('d', 'wddot')*s^2, freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(sys_rmc('d', 'wddot')*s^2, freqs, 'Hz'))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]'); xlabel('Frequency [Hz]');
  xlim([freqs(1), freqs(end)]);

  subplot(2, 2, 4);
  title('w to x')
  hold on;
  plot(freqs, abs(squeeze(freqresp(1+sys('d', 'wddot')*s^2, freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(1+sys_rmc('d', 'wddot')*s^2, freqs, 'Hz'))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]'); xlabel('Frequency [Hz]');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/rmc_1dof_sensitivitiy.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:rmc_1dof_sensitivitiy
#+CAPTION: Sensitivity to disturbance when RMC is applied on the 1dof system ([[./figs/rmc_1dof_sensitivitiy.png][png]], [[./figs/rmc_1dof_sensitivitiy.pdf][pdf]])
[[file:figs/rmc_1dof_sensitivitiy.png]]
** Inertial Control
*** One degree-of-freedom example
:PROPERTIES:
:header-args:matlab+: :tangle no
:END:
<<sec:ine_1dof>>
**** Equations
#+begin_src latex :file ine_1dof.pdf :post pdf2svg(file=*this*, ext="png") :exports results
  \begin{tikzpicture}
    % Ground
    \draw (-1, 0) -- (1, 0);

    % Ground Displacement
    \draw[dashed] (-1, 0) -- ++(-0.5, 0) coordinate(w);
    \draw[->] (w) -- ++(0, 0.5) node[left]{$w$};

    % Mass
    \draw[fill=white] (-1, 1) rectangle ++(2, 0.8) node[pos=0.5]{$m$};

    % Velocity Sensor
    \node[inertialsensor={0.3}] (velg) at (1, 1.8){};
    \node[above] at (velg.north) {$\dot{x}$};

    % Displacement of the mass
    \draw[dashed] (-1, 1.8) -- ++(-0.5, 0) coordinate(x);
    \draw[->] (x) -- ++(0, 0.5) node[left]{$x$};

    % Spring, Damper, and Actuator
    \draw[spring] (-0.8, 0) -- (-0.8, 1) node[midway, left=0.1]{$k$};
    \draw[damper] (0, 0)    -- (0, 1)    node[midway, left=0.2]{$c$};
    \draw[actuator={0.4}{0.2}] (0.8, 0) -- (0.8, 1)  coordinate[midway, right=0.1](F);

    % Control
    \node[block={0.8cm}{0.6cm}, right=0.6 of F] (Kine) {$K$};
    \draw[->] (Kine.west) -- (F) node[above right]{$F$};
    \draw[<-] (Kine.east) -- ++(0.5, 0) |- (velg.east);
  \end{tikzpicture}
#+end_src

#+name: fig:ine_1dof
#+caption: Direct Velocity Feedback applied to a 1dof system
#+RESULTS:
[[file:figs/ine_1dof.png]]

The dynamic of the system is:
\begin{equation}
  ms^2x = F_d - kx - csx + kw + csw + F
\end{equation}
In terms of the stage deformation $d = x - w$:
\begin{equation}
  (ms^2 + cs + k) d = -ms^2 w + F_d + F
\end{equation}
The direct velocity feedback law shown in figure [[fig:ine_1dof]] is:
\begin{equation}
  K = -g
\end{equation}
Thus, the applied force is:
\begin{equation}
  F = -g \dot{x}
\end{equation}
And the new dynamics will be:
\begin{equation}
  d = w \frac{-ms^2 - gs}{ms^2 + (c + g)s + k} + F_d \frac{1}{ms^2 + (c + g)s + k} + F \frac{1}{ms^2 + (c + g)s + k}
\end{equation}

And thus damping is added.

If critical damping is wanted:
\begin{equation}
  \xi = \frac{1}{2}\frac{c + g}{\sqrt{km}} = \frac{1}{2}
\end{equation}
This corresponds to a gain:
\begin{equation}
  g = \sqrt{km} - c
\end{equation}

**** Matlab Init                                          :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
  <<matlab-dir>>
#+end_src

#+begin_src matlab :exports none :results silent :noweb yes
  <<matlab-init>>
#+end_src

**** Matlab Example
Let define the system parameters.
#+begin_src matlab
  m = 50; % [kg]
  k = 1e6; % [N/m]
  c = 1e3; % [N/(m/s)]
#+end_src

The state space model of the system is defined below.
#+begin_src matlab
  A = [-c/m -k/m;
       1     0];

  B = [1/m 1/m -1;
       0   0    0];

  C = [1 0;
       0 1;
       0 0];

  D = [0 0 0;
       0 0 0;
       0 0 1];

  sys = ss(A, B, C, D);
  sys.InputName = {'F', 'Fd', 'wddot'};
  sys.OutputName = {'ddot', 'd', 'wddot'};
  sys.StateName = {'ddot', 'd'};
#+end_src

Because we need $\dot{x}$ for feedback, we compute it from the outputs
#+begin_src matlab
  G_xdot = [1, 0, 1/s;
            0, 1, 0];
  G_xdot.InputName = {'ddot', 'd', 'wddot'};
  G_xdot.OutputName = {'xdot', 'd'};
#+end_src

Finally, the system is described by =sys= as defined below.
#+begin_src matlab
  sys = series(sys, G_xdot, [1 2 3], [1 2 3]);
#+end_src

The controller $K_\text{INE}$ is:
#+begin_src matlab
  Kine = tf(-(sqrt(k*m)-c));
  Kine.InputName = {'xdot'};
  Kine.OutputName = {'F'};
#+end_src

And the closed loop system is computed below.
#+begin_src matlab
  sys_ine = feedback(sys, Kine, 'name', +1);
#+end_src

The obtained sensitivity to disturbances is shown in figure [[fig:ine_1dof_sensitivitiy]].
#+begin_src matlab :exports none
  freqs = logspace(-1, 3, 1000);

  figure;

  subplot(2, 2, 1);
  title('Fd to d')
  hold on;
  plot(freqs, abs(squeeze(freqresp(sys('d', 'Fd'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(sys_ine('d', 'Fd'), freqs, 'Hz'))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  xlim([freqs(1), freqs(end)]);

  subplot(2, 2, 3);
  title('Fd to x')
  hold on;
  plot(freqs, abs(squeeze(freqresp(sys('xdot', 'Fd')/s, freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(sys_ine('xdot', 'Fd')/s, freqs, 'Hz'))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  xlim([freqs(1), freqs(end)]);

  subplot(2, 2, 2);
  title('w to d')
  hold on;
  plot(freqs, abs(squeeze(freqresp(sys('d', 'wddot')*s^2, freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(sys_ine('d', 'wddot')*s^2, freqs, 'Hz'))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]'); xlabel('Frequency [Hz]');
  xlim([freqs(1), freqs(end)]);

  subplot(2, 2, 4);
  title('w to x')
  hold on;
  plot(freqs, abs(squeeze(freqresp(1+sys('d', 'wddot')*s^2, freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(1+sys_ine('d', 'wddot')*s^2, freqs, 'Hz'))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]'); xlabel('Frequency [Hz]');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/ine_1dof_sensitivitiy.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:ine_1dof_sensitivitiy
#+CAPTION: Sensitivity to disturbance when INE is applied on the 1dof system ([[./figs/ine_1dof_sensitivitiy.png][png]], [[./figs/ine_1dof_sensitivitiy.pdf][pdf]])
[[file:figs/ine_1dof_sensitivitiy.png]]