#+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

#+end_important

* Effect of granite stiffness
<<sec:granite_stiffness>>
** Analytical Analysis
#+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]{$x_{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.4*\massw, 0)   -- (-0.4*\massw, \spaceh) node[midway, left=0.1]{$k$};
      \draw[damper] (0, 0)             -- ( 0, \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.4*\massw, 0) -- (-0.4*\massw, \spaceh) node[midway, left=0.1]{$k^\prime$};
      \draw[damper] (0, 0)           -- ( 0, \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: Figure caption
#+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) (x_w - x)
\end{align}

If we note $d = x^\prime - x$, we obtain:
#+name: eq:plant_ground_transmissibility
\begin{equation}
  \frac{d}{x_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}

** Soft Granite
Let's initialize a soft granite that will act as an isolation stage from ground motion.
#+begin_src matlab
  initializeGranite('K', 5e5*ones(3,1), 'C', 5e3*ones(3,1));
#+end_src

#+begin_src matlab
  Ks = logspace(3,9,7); % [N/m]
#+end_src

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

#+begin_src matlab
  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
#+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', 'southeast');
#+end_src

#+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', 'southeast');
#+end_src


* 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