#+TITLE: Determination of the optimal nano-hexapod's stiffness for reducing the effect of disturbances
:DRAWER:
#+STARTUP: overview

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

#+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 ../matlab/optimal_stiffness.m
#+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/org/}{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 file raw replace
#+PROPERTY: header-args:latex+ :buffer no
#+PROPERTY: header-args:latex+ :eval no-export
#+PROPERTY: header-args:latex+ :exports results
#+PROPERTY: header-args:latex+ :mkdirp yes
#+PROPERTY: header-args:latex+ :output-dir figs
#+PROPERTY: header-args:latex+ :post pdf2svg(file=*this*, ext="png")
:END:

* Introduction                                                        :ignore:
In this document is studied how the stiffness of the nano-hexapod will impact the effect of disturbances on the position error of the sample.

It is divided in the following sections:
- Section [[sec:psd_disturbances]]: the disturbances are listed and their Power Spectral Densities (PSD) are shown
- Section [[sec:effect_disturbances]]: the transfer functions from disturbances to the position error of the sample are computed for a wide range of nano-hexapod stiffnesses
- Section [[sec:granite_stiffness]]:
- Section [[sec:open_loop_budget_error]]: from both the PSD of the disturbances and the transfer function from disturbances to sample's position errors, we compute the resulting PSD and Cumulative Amplitude Spectrum (CAS)
- Section [[sec:closed_loop_budget_error]]: from a simplistic model is computed the required control bandwidth to reduce the position error to acceptable values

* Disturbances
<<sec:psd_disturbances>>
** Introduction                                                      :ignore:
The main disturbances considered here are:
- $D_w$: Ground displacement in the $x$, $y$ and $z$ directions
- $F_{ty}$: Forces applied by the Translation stage in the $x$ and $z$ directions
- $F_{rz}$: Forces applied by the Spindle in the $z$ direction
- $F_d$: Direct forces applied at the center of mass of the Payload

The level of these disturbances has been identified form experiments which are detailed in [[file:disturbances.org][this]] document.

** 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
  load('mat/conf_simulink.mat');

  open('nass_model.slx')
#+end_src
** Plots                                                             :ignore:
The measured Amplitude Spectral Densities (ASD) of these forces are shown in Figures [[fig:opt_stiff_dist_gm]] and [[fig:opt_stiff_dist_fty_frz]].

In this study, the expected frequency content of the direct forces applied to the payload is not considered.

#+begin_src matlab :exports none
  load('./mat/dist_psd.mat', 'dist_f');
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(dist_f.f, sqrt(dist_f.psd_gm));
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('$\Gamma_{D_w}$ $\left[\frac{m}{\sqrt{Hz}}\right]$')
  xlim([1, 500]);
#+end_src

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

#+name: fig:opt_stiff_dist_gm
#+caption: Amplitude Spectral Density of the Ground Displacement ([[./figs/opt_stiff_dist_gm.png][png]], [[./figs/opt_stiff_dist_gm.pdf][pdf]])
[[file:figs/opt_stiff_dist_gm.png]]

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(dist_f.f, sqrt(dist_f.psd_ty), 'DisplayName', '$F_{T_y}$');
  plot(dist_f.f, sqrt(dist_f.psd_rz), 'DisplayName', '$F_{R_z}$');
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('ASD $\left[\frac{F}{\sqrt{Hz}}\right]$')
  legend('Location', 'southwest');
  xlim([1, 500]);
#+end_src

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

#+name: fig:opt_stiff_dist_fty_frz
#+caption: Amplitude Spectral Density of the "parasitic" forces comming from the Translation stage and the spindle ([[./figs/opt_stiff_dist_fty_frz.png][png]], [[./figs/opt_stiff_dist_fty_frz.pdf][pdf]])
[[file:figs/opt_stiff_dist_fty_frz.png]]

* Effect of disturbances on the position error
<<sec:effect_disturbances>>
** Introduction                                                      :ignore:
In this section, we use the Simscape model to identify the transfer function from disturbances to the position error of the sample.
We do that for a wide range of nano-hexapod stiffnesses and we compare the obtained results.

** 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
  load('mat/conf_simulink.mat');

  open('nass_model.slx')
#+end_src
** Initialization
We initialize all the stages with the default parameters.
#+begin_src matlab
  initializeGround();
  initializeGranite();
  initializeTy();
  initializeRy();
  initializeRz();
  initializeMicroHexapod();
  initializeAxisc();
  initializeMirror();
#+end_src

We use a sample mass of 10kg.
#+begin_src matlab
  initializeSample('mass', 10);
#+end_src

We include gravity, and we use no controller.
#+begin_src matlab
  initializeSimscapeConfiguration('gravity', true);
  initializeController();
  initializeDisturbances('enable', false);
  initializeLoggingConfiguration('log', 'none');
#+end_src

** Identification
The considered inputs are:
- =Dwx=: Ground displacement in the $x$ direction
- =Dwy=: Ground displacement in the $y$ direction
- =Dwz=: Ground displacement in the $z$ direction
- =Fty_x=: Forces applied by the Translation stage in the $x$ direction
- =Fty_z=: Forces applied by the Translation stage in the $z$ direction
- =Frz_z=: Forces applied by the Spindle in the $z$ direction
- =Fd=: Direct forces applied at the center of mass of the Payload

The outputs are =Ex=, =Ey=, =Ez=, =Erx=, =Ery=, =Erz= which are the 3 positions and 3 orientations errors of the sample.

We initialize the set of the nano-hexapod stiffnesses, and for each of them, we identify the dynamics from defined inputs to defined outputs.
#+begin_src matlab
  Ks = logspace(3,9,7); % [N/m]
#+end_src

#+begin_src matlab :exports none
  %% Name of the Simulink File
  mdl = 'nass_model';

  %% Micro-Hexapod
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Disturbances'], 1, 'openinput', [], 'Dwx');   io_i = io_i + 1; % X Ground motion
  io(io_i) = linio([mdl, '/Disturbances'], 1, 'openinput', [], 'Dwy');   io_i = io_i + 1; % Y Ground motion
  io(io_i) = linio([mdl, '/Disturbances'], 1, 'openinput', [], 'Dwz');   io_i = io_i + 1; % Z Ground motion
  io(io_i) = linio([mdl, '/Disturbances'], 1, 'openinput', [], 'Fty_x'); io_i = io_i + 1; % Parasitic force Ty - X
  io(io_i) = linio([mdl, '/Disturbances'], 1, 'openinput', [], 'Fty_z'); io_i = io_i + 1; % Parasitic force Ty - Z
  io(io_i) = linio([mdl, '/Disturbances'], 1, 'openinput', [], 'Frz_z'); io_i = io_i + 1; % Parasitic force Rz - Z
  io(io_i) = linio([mdl, '/Disturbances'], 1, 'openinput', [], 'Fd');    io_i = io_i + 1; % Direct forces

  io(io_i) = linio([mdl, '/Tracking Error'], 1, 'openoutput', [], 'En');   io_i = io_i + 1; % Position Error
#+end_src

#+begin_src matlab :exports none
  Gd = {zeros(length(Ks), 1)};

  for i = 1:length(Ks)
      initializeNanoHexapod('k', Ks(i));

      % Run the linearization
      G = linearize(mdl, io);
      G.InputName  = {'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_z', 'Fdx', 'Fdy', 'Fdz', 'Mdx', 'Mdy', 'Mdz'};
      G.OutputName = {'Ex', 'Ey', 'Ez', 'Erx', 'Ery', 'Erz'};
      Gd(i) = {minreal(G)};
  end
#+end_src

** Sensitivity to Stages vibration (Filtering)
The sensitivity the stage vibrations are displayed:
- Figure [[fig:opt_stiff_sensitivity_Frz]]: sensitivity to vertical spindle vibrations
- Figure [[fig:opt_stiff_sensitivity_Fty_z]]: sensitivity to vertical translation stage vibrations
- Figure [[fig:opt_stiff_sensitivity_Fty_x]]: sensitivity to horizontal (x) translation stage vibrations

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

  figure;
  hold on;
  for i = 1:length(Ks)
      plot(freqs, abs(squeeze(freqresp(Gd{i}('Ez', 'Frz_z'), freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Effect of $F_{rz}$ on $E_z$ [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'southwest');
#+end_src

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

#+name: fig:opt_stiff_sensitivity_Frz
#+caption: Sensitivity to Spindle vertical motion error ($F_{rz}$) to the vertical error position of the sample ($E_z$) ([[./figs/opt_stiff_sensitivity_Frz.png][png]], [[./figs/opt_stiff_sensitivity_Frz.pdf][pdf]])
[[file:figs/opt_stiff_sensitivity_Frz.png]]

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

  figure;
  hold on;
  for i = 1:length(Ks)
      plot(freqs, abs(squeeze(freqresp(Gd{i}('Ez', 'Fty_z'), freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Effect of $F_{ty}$ on $E_z$ [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'southwest');
#+end_src

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

#+name: fig:opt_stiff_sensitivity_Fty_z
#+caption: Sensitivity to Translation stage vertical motion error ($F_{ty,z}$) to the vertical error position of the sample ($E_z$) ([[./figs/opt_stiff_sensitivity_Fty_z.png][png]], [[./figs/opt_stiff_sensitivity_Fty_z.pdf][pdf]])
[[file:figs/opt_stiff_sensitivity_Fty_z.png]]

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

  figure;
  hold on;
  for i = 1:length(Ks)
      plot(freqs, abs(squeeze(freqresp(Gd{i}('Ex', 'Fty_x'), freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Effect of $F_{ty}$ on $E_x$ [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/opt_stiff_sensitivity_Fty_x.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+name: fig:opt_stiff_sensitivity_Fty_x
#+caption: Sensitivity to Translation stage $x$ motion error ($F_{ty,x}$) to the error position of the sample in the $x$ direction ($E_x$) ([[./figs/opt_stiff_sensitivity_Fty_x.png][png]], [[./figs/opt_stiff_sensitivity_Fty_x.pdf][pdf]])
[[file:figs/opt_stiff_sensitivity_Fty_x.png]]

** Effect of Ground motion (Transmissibility).
The effect of Ground motion on the position error of the sample is shown in Figure [[fig:opt_stiff_sensitivity_Dw]].

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

  figure;

  ax1 = subplot(1, 2, 1);
  hold on;
  for i = 1:length(Ks)
      plot(freqs, abs(squeeze(freqresp(Gd{i}('Ey', 'Dwy'), freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('$E_y/D_{wy}$ [m/m]'); xlabel('Frequency [Hz]');

  ax2 = subplot(1, 2, 2);
  hold on;
  for i = 1:length(Ks)
      plot(freqs, abs(squeeze(freqresp(Gd{i}('Ez', 'Dwz'), freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('$E_z/D_{wz}$ [m/m]'); xlabel('Frequency [Hz]');
#+end_src

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

#+name: fig:opt_stiff_sensitivity_Dw
#+caption: Sensitivity to Ground motion ($D_{w}$) to the position error of the sample ($E_y$ and $E_z$) ([[./figs/opt_stiff_sensitivity_Dw.png][png]], [[./figs/opt_stiff_sensitivity_Dw.pdf][pdf]])
[[file:figs/opt_stiff_sensitivity_Dw.png]]

** Direct Forces (Compliance).
The effect of direct forces/torques applied on the sample (cable forces for instance) on the position error of the sample is shown in Figure [[fig:opt_stiff_sensitivity_Fd]].

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

  figure;

  ax1 = subplot(1, 2, 1);
  hold on;
  for i = 1:length(Ks)
      plot(freqs, abs(squeeze(freqresp(Gd{i}('Ery', 'Mdy'), freqs, 'Hz'))), '-');
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('$E_{ry}/M_{d,y}\ \left[\frac{rad}{N m}\right]$'); xlabel('Frequency [Hz]');

  ax2 = subplot(1, 2, 2);
  hold on;
  for i = 1:length(Ks)
      plot(freqs, abs(squeeze(freqresp(Gd{i}('Ez', 'Fdz'), freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('$E_{z}/F_{d,z}$ [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'northeast');

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

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

#+name: fig:opt_stiff_sensitivity_Fd
#+caption: Sensitivity to Direct forces and torques applied to the sample ($F_d$, $M_d$) to the position error of the sample ([[./figs/opt_stiff_sensitivity_Fd.png][png]], [[./figs/opt_stiff_sensitivity_Fd.pdf][pdf]])
[[file:figs/opt_stiff_sensitivity_Fd.png]]

** Save                                                            :noexport:
#+begin_src matlab
  save('./mat/opt_stiffness_disturbances.mat', 'Ks', 'Gd')
#+end_src

** Conclusion
#+begin_important
  Reducing the nano-hexapod stiffness generally lowers the sensitivity to stages vibration but increases the sensitivity to ground motion and direct forces.

  In order to conclude on the optimal stiffness that will yield the smallest sample vibration, one has to include the level of disturbances. This is done in Section [[sec:open_loop_budget_error]].
#+end_important

* Effect of granite stiffness
<<sec:granite_stiffness>>
** Introduction                                                      :ignore:
In this section, we wish to see if a soft granite suspension could help in reducing the effect of disturbances on the position error of the sample.

** Analytical Analysis
*** Simple mass-spring-damper model
Let's consider the system shown in Figure [[fig:2dof_system_granite_stiffness]] consisting of two stacked mass-spring-damper systems.
The bottom one represents the granite, and the top one all the positioning stages.
We want the smallest stage "deformation" $d = x^\prime - x$ due to ground motion $w$.

#+begin_src latex :file 2dof_system_granite_stiffness.pdf
  \begin{tikzpicture}
    % ====================
    % Parameters
    % ====================
    \def\massw{2.2}  % Width of the masses
    \def\massh{0.8}  % Height of the masses
    \def\spaceh{1.2} % Height of the springs/dampers
    \def\dispw{0.3}  % Width of the dashed line for the displacement
    \def\disph{0.5}  % Height of the arrow for the displacements
    \def\bracs{0.05} % Brace spacing vertically
    \def\brach{-10pt} % Brace shift horizontaly
    % ====================


    % ====================
    % Ground
    % ====================
    \draw (-0.5*\massw, 0) -- (0.5*\massw, 0);
    \draw[dashed] (0.5*\massw, 0) -- ++(\dispw, 0) coordinate(dlow);
    \draw[->] (0.5*\massw+0.5*\dispw, 0) -- ++(0, \disph) node[right]{$w$};

    % ====================
    % Micro Station
    % ====================
    \begin{scope}[shift={(0, 0)}]
      % Mass
      \draw[fill=white] (-0.5*\massw, \spaceh) rectangle (0.5*\massw, \spaceh+\massh) node[pos=0.5]{$m$};

      % Spring, Damper, and Actuator
      \draw[spring] (-0.3*\massw, 0) -- (-0.3*\massw, \spaceh) node[midway, left=0.1]{$k$};
      \draw[damper] ( 0.3*\massw, 0) -- ( 0.3*\massw, \spaceh) node[midway, left=0.2]{$c$};

      % Displacements
      \draw[dashed] (0.5*\massw, \spaceh) -- ++(\dispw, 0);
      \draw[->] (0.5*\massw+0.5*\dispw, \spaceh) -- ++(0, \disph) node[right]{$x$};

      % Legend
      \draw[decorate, decoration={brace, amplitude=8pt}, xshift=\brach] %
      (-0.5*\massw, \bracs) -- (-0.5*\massw, \spaceh+\massh-\bracs) %
      node[midway,rotate=90,anchor=south,yshift=10pt,align=center]{Granite};
    \end{scope}

    % ====================
    % Nano Station
    % ====================
    \begin{scope}[shift={(0, \spaceh+\massh)}]
      % Mass
      \draw[fill=white] (-0.5*\massw, \spaceh) rectangle (0.5*\massw, \spaceh+\massh) node[pos=0.5]{$m^\prime$};

      % Spring, Damper, and Actuator
      \draw[spring] (-0.3*\massw, 0) -- (-0.3*\massw, \spaceh) node[midway, left=0.1]{$k^\prime$};
      \draw[damper] ( 0.3*\massw, 0) -- ( 0.3*\massw, \spaceh) node[midway, left=0.2]{$c^\prime$};

      % Displacements
      \draw[dashed] (0.5*\massw, \spaceh) -- ++(\dispw, 0) coordinate(dhigh);
      \draw[->] (0.5*\massw+0.5*\dispw, \spaceh) -- ++(0, \disph) node[right]{$x^\prime$};

      % Legend
      \draw[decorate, decoration={brace, amplitude=8pt}, xshift=\brach] %
      (-0.5*\massw, \bracs) -- (-0.5*\massw, \spaceh+\massh-\bracs) %
      node[midway,rotate=90,anchor=south,yshift=10pt,align=center]{Positioning\\Stages};
    \end{scope}
  \end{tikzpicture}
#+end_src

#+name: fig:2dof_system_granite_stiffness
#+caption: Mass Spring Damper system consisting of a granite and a positioning stage
#+RESULTS:
[[file:figs/2dof_system_granite_stiffness.png]]

If we write the equation of motion of the system in Figure [[fig:2dof_system_granite_stiffness]], we obtain:
\begin{align}
  m^\prime s^2 x^\prime &= (c^\prime s + k^\prime) (x - x^\prime) \\
  ms^2 x &= (c^\prime s + k^\prime) (x^\prime - x) + (cs + k) (w - x)
\end{align}

If we note $d = x^\prime - x$, we obtain:
\begin{equation}
  \frac{d}{w} = \frac{-m^\prime s^2 (cs + k)}{ (m^\prime s^2 + c^\prime s + k^\prime) (ms^2 + cs + k) + m^\prime s^2(c^\prime s + k^\prime)}
\end{equation}

*** General Case
Let's now considering a general positioning stage defined by:
- $G^\prime(s) = \frac{F}{x}$: its mechanical "impedance"
- $D^\prime(s) = \frac{d}{x}$: its "deformation" transfer function

#+begin_src latex :file general_system_granite_stiffness.pdf
  \begin{tikzpicture}
    \def\massw{2.2}  % Width of the masses
    \def\massh{0.8}  % Height of the masses
    \def\spaceh{1.2} % Height of the springs/dampers
    \def\dispw{0.3}  % Width of the dashed line for the displacement
    \def\disph{0.5}  % Height of the arrow for the displacements
    \def\bracs{0.05} % Brace spacing vertically
    \def\brach{-10pt} % Brace shift horizontaly

    % Mass
    \draw[fill=white] (-0.5*\massw, \spaceh) rectangle node[left=6pt]{$m$} (0.5*\massw, \spaceh+\massh);

    % Spring, Damper, and Actuator
    \draw[spring] (-0.3*\massw, 0) -- (-0.3*\massw, \spaceh) node[midway, left=0.1]{$k$};
    \draw[damper] ( 0.3*\massw, 0) -- ( 0.3*\massw, \spaceh) node[midway, left=0.2]{$c$};

    % Ground
    \draw (-0.5*\massw, 0) -- (0.5*\massw, 0);
    % Groud Motion
    \draw[dashed] (0.5*\massw, 0) -- ++(\dispw, 0);
    \draw[->] (0.5*\massw+0.5*\dispw, 0) -- ++(0, \disph) node[right]{$w$};

    % Displacements
    \draw[dashed] (0.5*\massw, \spaceh+\massh) -- ++(2*\dispw, 0) coordinate(dhigh);
    \draw[->] (0.5*\massw+1.5*\dispw, \spaceh+\massh) -- ++(0, \disph) node[right]{$x$};

    % Legend
    \draw[decorate, decoration={brace, amplitude=8pt}, xshift=\brach] %
      (-0.5*\massw, \bracs) -- (-0.5*\massw, \spaceh+\massh-\bracs) %
      node[midway,rotate=90,anchor=south,yshift=10pt,align=center]{Granite};

    \begin{scope}[shift={(0, \spaceh+\massh)}]
      \node[piezo={2.2}{1.5}{6}, anchor=south] (piezo) at (0, 0){};
      \draw[->] (0,0)node[branch]{} -- ++(0, -0.6)node[above right]{$F$}

      \draw[<->] (1.1+0.5*\dispw,0) -- node[midway, right]{$d$} ++(0,1.5);

      \draw[decorate, decoration={brace, amplitude=8pt}, xshift=\brach] %
        ($(piezo.south west) + (-10pt, 0)$) -- ($(piezo.north west) + (-10pt, 0)$) %
        node[midway,rotate=90,anchor=south,yshift=10pt,align=center]{Positioning\\Stages};
    \end{scope}
  \end{tikzpicture}
#+end_src

#+name: fig:general_system_granite_stiffness
#+caption: Mass Spring Damper representing the granite and a general representation of positioning stages
#+RESULTS:
[[file:figs/general_system_granite_stiffness.png]]

The equation of motion are:
\begin{align}
  ms^2 x &= (cs + k) (x - w) - F \\
  F &= G^\prime(s) x \\
  d &= D^\prime(s) x
\end{align}
where:
- $F$ is the force applied by the position stages on the granite mass

#+begin_important
We can express $d$ as a function of $w$:
\begin{equation}
  \frac{d}{w} = \frac{D^\prime(s) (cs + k)}{ms^2 + cs + k + G^\prime(s)}
\end{equation}

This is the transfer function that we would like to minimize.
#+end_important

Let's verify this formula for a simple mass/spring/damper positioning stage.
In that case, we have:
\begin{align*}
  D^\prime(s) &= \frac{d}{x} = \frac{- m^\prime s^2}{m^\prime s^2 + c^\prime s + k^\prime} \\
  G^\prime(s) &= \frac{F}{x} = \frac{m^\prime s^2(c^\prime s + k)}{m^\prime s^2 + c^\prime s + k^\prime}
\end{align*}

And finally:
\begin{equation}
  \frac{d}{w} = \frac{-m^\prime s^2 (cs + k)}{ (m^\prime s^2 + c^\prime s + k^\prime) (ms^2 + cs + k) + m^\prime s^2(c^\prime s + k^\prime)}
\end{equation}
which is the same as the previously derived equation.

** Soft Granite
Let's initialize a soft granite and see how the sensitivity to disturbances will change.
#+begin_src matlab
  initializeGranite('K', 5e5*ones(3,1), 'C', 5e3*ones(3,1));
#+end_src

#+begin_src matlab :exports none
  Gdr = {zeros(length(Ks), 1)};

  for i = 1:length(Ks)
      initializeNanoHexapod('k', Ks(i));

      G = linearize(mdl, io);
      G.InputName  = {'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_z', 'Fdx', 'Fdy', 'Fdz', 'Mdx', 'Mdy', 'Mdz'};
      G.OutputName = {'Ex', 'Ey', 'Ez', 'Erx', 'Ery', 'Erz'};
      Gdr(i) = {minreal(G)};
  end
#+end_src

** Effect of the Granite transfer function
From Figure [[fig:opt_stiff_soft_granite_Dw]], we can see that having a "soft" granite suspension greatly lowers the sensitivity to ground motion.
The sensitivity is indeed lowered starting from the resonance of the granite on its soft suspension (few Hz here).

From Figures [[fig:opt_stiff_soft_granite_Frz]] and [[fig:opt_stiff_soft_granite_Fd]], we see that the change of granite suspension does not change a lot the sensitivity to both direct forces and stage vibrations.

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

  figure;
  hold on;
  for i = 1:length(Ks)
      set(gca,'ColorOrderIndex',i);
      plot(freqs, abs(squeeze(freqresp(Gd{i}('Ez', 'Dwz'), freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
      set(gca,'ColorOrderIndex',i);
      plot(freqs, abs(squeeze(freqresp(Gdr{i}('Ez', 'Dwz'), freqs, 'Hz'))), '--', ...
           'HandleVisibility', 'off');
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]'); xlabel('Frequency [Hz]');
  legend('location', 'southwest');
#+end_src

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

#+name: fig:opt_stiff_soft_granite_Dw
#+caption: Change of sensibility to Ground motion when using a stiff Granite (solid curves) and a soft Granite (dashed curves) ([[./figs/opt_stiff_soft_granite_Dw.png][png]], [[./figs/opt_stiff_soft_granite_Dw.pdf][pdf]])
[[file:figs/opt_stiff_soft_granite_Dw.png]]

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

  figure;
  hold on;
  for i = 1:length(Ks)
      set(gca,'ColorOrderIndex',i);
      plot(freqs, abs(squeeze(freqresp(Gd{i}('Ez', 'Frz_z'), freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
      set(gca,'ColorOrderIndex',i);
      plot(freqs, abs(squeeze(freqresp(Gdr{i}('Ez', 'Frz_z'), freqs, 'Hz'))), '--', ...
           'HandleVisibility', 'off');
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
  legend('location', 'southwest');
#+end_src

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

#+name: fig:opt_stiff_soft_granite_Frz
#+caption: Change of sensibility to Spindle vibrations when using a stiff Granite (solid curves) and a soft Granite (dashed curves) ([[./figs/opt_stiff_soft_granite_Frz.png][png]], [[./figs/opt_stiff_soft_granite_Frz.pdf][pdf]])
[[file:figs/opt_stiff_soft_granite_Frz.png]]

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

  figure;
  hold on;
  for i = 1:length(Ks)
      set(gca,'ColorOrderIndex',i);
      plot(freqs, abs(squeeze(freqresp(Gd{i}( 'Ez', 'Fdz'), freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
      set(gca,'ColorOrderIndex',i);
      plot(freqs, abs(squeeze(freqresp(Gdr{i}('Ez', 'Fdz'), freqs, 'Hz'))), '--', ...
           'HandleVisibility', 'off');
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('$E_{z}/F_{d,z}$ [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/opt_stiff_soft_granite_Fd.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src

#+name: fig:opt_stiff_soft_granite_Fd
#+caption: Change of sensibility to direct forces when using a stiff Granite (solid curves) and a soft Granite (dashed curves) ([[./figs/opt_stiff_soft_granite_Fd.png][png]], [[./figs/opt_stiff_soft_granite_Fd.pdf][pdf]])
[[file:figs/opt_stiff_soft_granite_Fd.png]]

** Conclusion
#+begin_important
  Having a soft granite suspension could greatly improve the sensitivity the ground motion and thus the level of sample vibration if it is found that ground motion is the limiting factor.
#+end_important

* Open Loop Budget Error
<<sec:open_loop_budget_error>>
** 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 :tangle no
  simulinkproject('../');
#+end_src

#+begin_src matlab
  load('mat/conf_simulink.mat');

  open('nass_model.slx')
#+end_src

** Load of the identified disturbances and transfer functions
#+begin_src matlab
  load('./mat/dist_psd.mat', 'dist_f');
  load('./mat/opt_stiffness_disturbances.mat', 'Gd')
#+end_src

** Equations

** Results
Effect of all disturbances

#+begin_src matlab
  freqs = dist_f.f;

  figure;
  hold on;
  for i = 1:length(Ks)
    plot(freqs, sqrt(dist_f.psd_rz).*abs(squeeze(freqresp(Gd{i}('Ez', 'Frz_z'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('ASD $\left[\frac{m}{\sqrt{Hz}}\right]$')
  legend('Location', 'southwest');
  xlim([2, 500]);
#+end_src

** Cumulative Amplitude Spectrum
#+begin_src matlab
  freqs = dist_f.f;

  figure;
  hold on;
  for i = 1:length(Ks)
    plot(freqs, sqrt(flip(-cumtrapz(flip(freqs), flip(dist_f.psd_ty.*abs(squeeze(freqresp(Gd{i}('Ez', 'Fty_z'), freqs, 'Hz'))).^2)))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
  end
  plot([freqs(1) freqs(end)], [10e-9 10e-9], 'k--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('CAS $[m]$')
  legend('Location', 'southwest');
  xlim([2, 500]); ylim([1e-10 1e-6]);
#+end_src

#+begin_src matlab
  freqs = dist_f.f;

  figure;
  hold on;
  for i = 1:length(Ks)
    plot(freqs, sqrt(flip(-cumtrapz(flip(freqs), flip(dist_f.psd_rz.*abs(squeeze(freqresp(Gd{i}('Ez', 'Frz_z'), freqs, 'Hz'))).^2)))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
  end
  plot([freqs(1) freqs(end)], [10e-9 10e-9], 'k--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('CAS $[m]$')
  legend('Location', 'southwest');
  xlim([2, 500]); ylim([1e-10 1e-6]);
#+end_src

Ground motion
#+begin_src matlab
  freqs = dist_f.f;

  figure;
  hold on;
  for i = 1:length(Ks)
    plot(freqs, sqrt(flip(-cumtrapz(flip(freqs), flip(dist_f.psd_gm.*abs(squeeze(freqresp(Gd{i}('Ez', 'Dwz'), freqs, 'Hz'))).^2)))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
  end
  plot([freqs(1) freqs(end)], [10e-9 10e-9], 'k--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('CAS $E_y$ $[m]$')
  legend('Location', 'northeast');
  xlim([2, 500]); ylim([1e-10 1e-6]);
#+end_src

#+begin_src matlab
  freqs = dist_f.f;

  figure;
  hold on;
  for i = 1:length(Ks)
    plot(freqs, sqrt(flip(-cumtrapz(flip(freqs), flip(dist_f.psd_gm.*abs(squeeze(freqresp(Gd{i}('Ex', 'Dwx'), freqs, 'Hz'))).^2)))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
  end
  plot([freqs(1) freqs(end)], [10e-9 10e-9], 'k--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'lin');
  xlabel('Frequency [Hz]'); ylabel('CAS $E_y$ $[m]$')
  legend('Location', 'northeast');
  xlim([2, 500]);
#+end_src

#+begin_src matlab
  freqs = dist_f.f;

  figure;
  hold on;
  for i = 1:length(Ks)
    plot(freqs, sqrt(flip(-cumtrapz(flip(freqs), flip(dist_f.psd_gm.*abs(squeeze(freqresp(Gd{i}('Ey', 'Dwy'), freqs, 'Hz'))).^2)))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
  end
  plot([freqs(1) freqs(end)], [10e-9 10e-9], 'k--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'lin');
  xlabel('Frequency [Hz]'); ylabel('CAS $E_y$ $[m]$')
  legend('Location', 'northeast');
  xlim([2, 500]);
#+end_src

Sum of all perturbations
#+begin_src matlab
  psd_tot = zeros(length(freqs), length(Ks));

  for i = 1:length(Ks)
      psd_tot(:,i) = dist_f.psd_gm.*abs(squeeze(freqresp(Gd{i}('Ez', 'Dwz'  ), freqs, 'Hz'))).^2 + ...
          dist_f.psd_ty.*abs(squeeze(freqresp(Gd{i}('Ez', 'Fty_z'), freqs, 'Hz'))).^2 + ...
          dist_f.psd_rz.*abs(squeeze(freqresp(Gd{i}('Ez', 'Frz_z'), freqs, 'Hz'))).^2;
  end
#+end_src

#+begin_src matlab
  freqs = dist_f.f;

  figure;
  hold on;
  for i = 1:length(Ks)
    plot(freqs, sqrt(flip(-cumtrapz(flip(freqs), flip(psd_tot(:,i))))), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)));
  end
  plot([freqs(1) freqs(end)], [10e-9 10e-9], 'k--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('CAS $E_z$ $[m]$')
  legend('Location', 'northeast');
  xlim([1, 500]); ylim([1e-10 1e-6]);
#+end_src

* Closed Loop Budget Error
<<sec:closed_loop_budget_error>>
** Introduction                                                      :ignore:
** Reduction thanks to feedback - Required bandwidth
#+begin_src matlab
  wc = 1*2*pi; % [rad/s]
  xic = 0.5;

  S = (s/wc)/(1 + s/wc);

  bodeFig({S}, logspace(-1,2,1000))
#+end_src

#+begin_src matlab
  wc = [1, 5, 10, 20, 50, 100, 200];

  S1 = {zeros(length(wc), 1)};
  S2 = {zeros(length(wc), 1)};

  for j = 1:length(wc)
      L = (2*pi*wc(j))/s; % Simple integrator
      S1{j} = 1/(1 + L);
      L = ((2*pi*wc(j))/s)^2*(1 + s/(2*pi*wc(j)/2))/(1 + s/(2*pi*wc(j)*2));
      S2{j} = 1/(1 + L);
  end
#+end_src

#+begin_src matlab
  freqs = dist_f.f;

  figure;
  hold on;
  i = 6;
  for j = 1:length(wc)
      set(gca,'ColorOrderIndex',j);
      plot(freqs, sqrt(flip(-cumtrapz(flip(freqs), flip(abs(squeeze(freqresp(S1{j}, freqs, 'Hz'))).^2.*psd_tot(:,i))))), '-', ...
           'DisplayName', sprintf('$\\omega_c = %.0f$ [Hz]', wc(j)));
  end
  plot(freqs, sqrt(flip(-cumtrapz(flip(freqs), flip(psd_tot(:,i))))), 'k-', ...
       'DisplayName', 'Open-Loop');
  plot([freqs(1) freqs(end)], [10e-9 10e-9], 'k--', 'HandleVisibility', 'off');
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('CAS $E_y$ $[m]$')
  legend('Location', 'northeast');
  xlim([0.5, 500]); ylim([1e-10 1e-6]);
#+end_src

#+begin_src matlab
  wc = logspace(0, 3, 100);

  Dz1_rms = zeros(length(Ks), length(wc));
  Dz2_rms = zeros(length(Ks), length(wc));
  for i = 1:length(Ks)
      for j = 1:length(wc)
          L = (2*pi*wc(j))/s;
          Dz1_rms(i, j) = sqrt(trapz(freqs, abs(squeeze(freqresp(1/(1 + L), freqs, 'Hz'))).^2.*psd_tot(:,i)));

          L = ((2*pi*wc(j))/s)^2*(1 + s/(2*pi*wc(j)/2))/(1 + s/(2*pi*wc(j)*2));
          Dz2_rms(i, j) = sqrt(trapz(freqs, abs(squeeze(freqresp(1/(1 + L), freqs, 'Hz'))).^2.*psd_tot(:,i)));
      end
  end
#+end_src

#+begin_src matlab
  freqs = dist_f.f;

  figure;
  hold on;
  for i = 1:length(Ks)
    set(gca,'ColorOrderIndex',i);
    plot(wc, Dz1_rms(i, :), '-', ...
           'DisplayName', sprintf('$k = %.0g$ [N/m]', Ks(i)))

    set(gca,'ColorOrderIndex',i);
    plot(wc, Dz2_rms(i, :), '--', ...
           'HandleVisibility', 'off')
  end
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  xlabel('Control Bandwidth [Hz]'); ylabel('$E_z\ [m, rms]$')
  legend('Location', 'southwest');
  xlim([1, 500]);
#+end_src
* Conclusion