nass-simscape/org/uncertainty_support.org

#+TITLE: Effect of Uncertainty on the support's dynamics on the isolation platform dynamics
#+SETUPFILE: ./setup/org-setup-file.org

* Introduction                                                        :ignore:
In this document we will consider an *isolation platform* (e.g. the nano-hexapod) on top of a *flexible support* (e.g. the micro-station).

The goal is to study:
- how does the dynamics of the support influence the dynamics of the plant to control
- similarly: how does the uncertainty on the support's dynamics will be transferred to uncertainty on the plant
- what design choice should be made in order to minimize the resulting uncertainty on the plant

Two models are made to study these effects:
- In section [[sec:introductory_example]], simple mass-spring-damper systems are chosen to model both the isolation platform and the flexible support
- In section [[sec:arbitrary_dynamics]], we consider arbitrary support dynamics with multiplicative input uncertainty to study the unmodelled dynamics of the support

* Simple Introductory Example
<<sec:introductory_example>>

** Introduction                                                     :ignore:
Let's consider the system shown in Figure [[fig:2dof_system_stiffness_uncertainty]] consisting of:
- A *support* represented by a mass $m^\prime$, a stiffness $k^\prime$ and a dashpot $c^\prime$
- An *isolation platform* represented by a mass $m$, a stiffness $k$ and a dashpot $c$ and an actuator $F$

The goal is to stabilize $x$ using $F$ in spite of uncertainty on the support mechanical properties.

#+begin_src latex :file 2dof_system_stiffness_uncertainty.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^\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);
      \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]{Support};
    \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$};

      % 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$};
      \draw[actuator] ( 0.4*\massw, 0) -- (0.4*\massw, \spaceh) node[midway, left=0.1](F){$F$};

      % Displacements
      \draw[dashed] (0.5*\massw, \spaceh) -- ++(\dispw, 0) coordinate(dhigh);
      \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]{Isolation\\Platform};
    \end{scope}

    % \draw[<->] (dlow) -- node[midway, right]{$d$} (dhigh);
  \end{tikzpicture}
#+end_src

#+name: fig:2dof_system_stiffness_uncertainty
#+caption: Two degrees-of-freedom system
#+RESULTS:
[[file:figs/2dof_system_stiffness_uncertainty.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

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

After eliminating $x^\prime$, we obtain:
#+name: eq:plant_simple_system
\begin{equation}
  \frac{x}{F} = \frac{m^\prime s^2 + c^\prime s + k^\prime}{ms^2(cs + k) + (ms^2 + cs + k)(m^\prime s^2 + c^\prime s + k^\prime)}
\end{equation}

** Initialization of the support dynamics
Let the support have:
- a nominal mass of $m^\prime = 1000\ [kg]$
- a nominal stiffness of $k^\prime = 10^8\ [N/m]$
- a nominal damping of $c^\prime = 10^5\ [N/(m/s)]$

#+begin_src matlab
  mpi = 1e3;
  cpi = 5e4;
  kpi = 1e8;
#+end_src

Let's also consider some uncertainty in those parameters:
#+begin_src matlab
  mp = ureal('m', mpi, 'Percentage', 30);
  cp = ureal('c', cpi, 'Percentage', 30);
  kp = ureal('k', kpi, 'Percentage', 30);
#+end_src

The compliance of the support without the isolation platform is $\frac{1}{m^\prime s^2 + c^\prime s + k^\prime}$ and its bode plot is shown in Figure [[fig:nominal_support_compliance_dynamics]].

One can see that support has a resonance frequency of $\omega_0^\prime = 50\ Hz$.

#+begin_src matlab :exports none
  Gps = usample(1/(mp*s^2 + cp*s + kp), 10);

  freqs = logspace(1, 3, 1000);

  figure;

  hold on;
  for i = 1:length(Gps)
    plot(freqs, abs(squeeze(freqresp(Gps(:,:,i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  plot(freqs, abs(squeeze(freqresp(1/(mpi*s^2 + cpi*s + kpi), freqs, 'Hz'))), 'k-');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]');
  ylabel('Magnitude [dB]');
  xlim([freqs(1), freqs(end)]);
#+end_src

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

#+name: fig:nominal_support_compliance_dynamics
#+caption: Nominal compliance of the support ([[./figs/nominal_support_compliance_dynamics.png][png]], [[./figs/nominal_support_compliance_dynamics.pdf][pdf]])
[[file:figs/nominal_support_compliance_dynamics.png]]

** Initialization of the isolation platform
Let's first fix the mass of the payload to be isolated:
#+begin_src matlab
  m = 100;
#+end_src

And we generate three isolation platforms:
- A soft one with $\omega_0 = 0.1 \omega_0^\prime = 5\ Hz$
- A medium stiff one with $\omega_0 = \omega_0^\prime = 50\ Hz$
- A stiff one with $\omega_0 = 10 \omega_0^\prime = 500\ Hz$

#+begin_src matlab :exports none
  k_soft = m*(2*pi*5)^2;
  c_soft = 0.1*sqrt(m*k_soft);
  G_soft = (mp*s^2 + cp*s + kp)/(m*s^2*(c_soft*s + k_soft) + (m*s^2 + c_soft*s + k_soft)*(mp*s^2 + cp*s + kp));

  k_mid = m*(2*pi*50)^2;
  c_mid = 0.1*sqrt(m*k_mid);
  G_mid = (mp*s^2 + cp*s + kp)/(m*s^2*(c_mid*s + k_mid) + (m*s^2 + c_mid*s + k_mid)*(mp*s^2 + cp*s + kp));

  k_stiff = m*(2*pi*500)^2;
  c_stiff = 0.1*sqrt(m*k_stiff);
  G_stiff = (mp*s^2 + cp*s + kp)/(m*s^2*(c_stiff*s + k_stiff) + (m*s^2 + c_stiff*s + k_stiff)*(mp*s^2 + cp*s + kp));
#+end_src

** Comparison
The obtained dynamics from $F$ to $x$ for the three isolation platform are shown in Figure [[fig:plant_dynamics_uncertainty_stiff_mid_soft]].

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

  Gs_soft  = usample(G_soft,  10);
  Gs_mid   = usample(G_mid,   10);
  Gs_stiff = usample(G_stiff, 10);

  figure;

  ax1 = subplot(2,3,1);
  hold on;
  for i = 1:length(Gs_soft)
    plot(freqs, abs(squeeze(freqresp(Gs_soft(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  ylabel('Magnitude [m/N]');
  title('$\omega_0 \ll \omega_0^\prime$');
  hold off;

  ax4 = subplot(2,3,4);
  hold on;
  for i = 1:length(Gs_soft)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gs_soft(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  set(gca,'xscale','log');
  yticks(-360:90:180);
  ylim([-270 90]);
  ylabel('Phase [deg]');
  hold off;

  ax2 = subplot(2,3,2);
  hold on;
  for i = 1:length(Gs_mid)
    plot(freqs, abs(squeeze(freqresp(Gs_mid(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  title('$\omega_0 \approx \omega_0^\prime$');
  hold off;

  ax5 = subplot(2,3,5);
  hold on;
  for i = 1:length(Gs_mid)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gs_mid(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  set(gca,'xscale','log');
  yticks(-360:90:180);
  ylim([-270 90]);
  xlabel('Frequency [Hz]');
  hold off;

  ax3 = subplot(2,3,3);
  hold on;
  for i = 1:length(Gs_stiff)
    plot(freqs, abs(squeeze(freqresp(Gs_stiff(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  title('$\omega_0 \gg \omega_0^\prime$');
  hold off;

  ax6 = subplot(2,3,6);
  hold on;
  for i = 1:length(Gs_stiff)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gs_stiff(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  set(gca,'xscale','log');
  yticks(-360:90:180);
  ylim([-270 90]);
  hold off;

  linkaxes([ax1,ax2,ax3,ax4,ax5,ax6],'x');
  xlim([freqs(1), freqs(end)]);

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

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

#+name: fig:plant_dynamics_uncertainty_stiff_mid_soft
#+caption: Obtained plant for the three isolation platforms considered ([[./figs/plant_dynamics_uncertainty_stiff_mid_soft.png][png]], [[./figs/plant_dynamics_uncertainty_stiff_mid_soft.pdf][pdf]])
[[file:figs/plant_dynamics_uncertainty_stiff_mid_soft.png]]

** Conclusion
#+begin_important
  The soft platform dynamics does not seems to depend on the dynamics of the support nor to be affect by the dynamic uncertainty of the support.
#+end_important

* Generalization to arbitrary dynamics
<<sec:arbitrary_dynamics>>
** Introduction
Let's now consider a general support described by its *compliance* $G^\prime(s) = \frac{x^\prime}{F^\prime}$ as shown in Figure [[fig:general_support_compliance]].

#+begin_src latex :file general_support_compliance.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

    \node[piezo={2.2}{3}{10}] (piezo) at (0, 0){};
    \draw[] ($(piezo.north)+(-1.2, 0)$) -- ++(2.4, 0);
    \draw[] ($(piezo.south)+(-1.2, 0)$) -- ++(2.4, 0);
    \draw[dashed] (piezo.north east) -- ++(\dispw, 0) coordinate(dhigh);
    \draw[->] ($(piezo.north east)+(0.5*\dispw, 0)$) -- ++(0, \disph) node[right]{$x^\prime$};
    \draw[->] (piezo.north) node[branch]{} -- ++(0, 1) node[below right]{$F^\prime$};

    \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]{General Support};
  \end{tikzpicture}
#+end_src

#+name: fig:general_support_compliance
#+caption: General support
#+RESULTS:
[[file:figs/general_support_compliance.png]]

Now let's consider the system consisting of a mass-spring-system (the isolation platform) on top of a general support as shown in Figure [[fig:general_support_with_isolator]].
#+begin_src latex :file general_support_with_isolator.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

    \node[piezo={2.2}{3}{10}] (piezo) at (0, 0){};
    \draw[] ($(piezo.north)+(-1.2, 0)$) -- ++(2.4, 0);
    \draw[] ($(piezo.south)+(-1.2, 0)$) -- ++(2.4, 0);
    \draw[dashed] (piezo.north east) -- ++(\dispw, 0) coordinate(dhigh);
    \draw[->] ($(piezo.north east)+(0.5*\dispw, 0)$) -- ++(0, \disph) node[right]{$x^\prime$};
    \draw[->] (piezo.north) node[branch]{} -- ++(0, -1) node[above right]{$-F^\prime$};

    \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]{General Support};

    \begin{scope}[shift={(piezo.north)}]
      % 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$};
      \draw[actuator] ( 0.4*\massw, 0) -- (0.4*\massw, \spaceh) node[midway, left=0.1](F){$F$};

      % Displacements
      \draw[dashed] (0.5*\massw, \spaceh) -- ++(\dispw, 0) coordinate(dhigh);
      \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]{Isolation\\Platform};
    \end{scope}
  \end{tikzpicture}
#+end_src

#+name: fig:general_support_with_isolator
#+caption: Mass-Spring-Damper system on top of a general support
#+RESULTS:
[[file:figs/general_support_with_isolator.png]]

** Equations of motion
We have to following equations of motion:
\begin{align}
  ms^2 x &= F + (cs + k) (x^\prime - x) \\
  F^\prime &= -F +  (cs + k)(x - x^\prime) \\
  \frac{x^\prime}{F^\prime} &= G^\prime(s)
\end{align}

And by eliminating $F^\prime$ and $x^\prime$, we find the plant dynamics $G(s) = \frac{x}{F}$.

#+begin_important
\begin{equation}
  \frac{x}{F} = \frac{1}{ms^2 + cs + k + ms^2(cs + k)G^\prime(s)} \label{eq:plant_dynamics_general_support}
\end{equation}
#+end_important

In order to verify that the formula is correct, let's take the same mass-spring-damper system used in the system shown in Figure [[fig:2dof_system_stiffness_uncertainty]]:
\[ \frac{x^\prime}{F^\prime} = \frac{1}{m^\prime s^2 + c^\prime s + k^\prime} \]

And we obtain
\[ \frac{x}{F} = \frac{m^\prime s^2 + c^\prime s + k^\prime}{(ms^2 + cs + k)(m^\prime s^2 + c^\prime s + k^\prime) + ms^2(cs + k)} \]
Which is the same transfer function that was obtained in section [[sec:introductory_example]] (Eq. eqref:eq:plant_simple_system).

** 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

** Compliance of the Support
We model the support by a mass-spring-damper model with some uncertainty.

The nominal compliance of the support is corresponding to the compliance of a mass-spring-damper system with a mass of $1000\ kg$ and a stiffness of $10^8\ [N/m]$.
The main resonance of the support is then $\omega^\prime = \sqrt{\frac{m^\prime}{k^\prime}} \approx 50\ Hz$.

#+begin_src matlab
  m0 = 1e3;
  c0 = 5e4;
  k0 = 1e8;

  Gp0 = 1/(m0*s^2 + c0*s + k0);
#+end_src

Let's represent the uncertainty on the compliance of the support by a multiplicative uncertainty (Figure [[fig:input_uncertainty_set]]):
\[ G^\prime(s) = G_0^\prime(s)(1 + w_I^\prime(s)\Delta_I(s)) \quad |\Delta_I(j\omega)| < 1\ \forall \omega \]

This could represent *unmodelled dynamics* or unknown parameters of the support.

#+begin_src latex :file input_uncertainty_set.pdf
  \begin{tikzpicture}
    \tikzset{block/.default={0.8cm}{0.8cm}}
    \tikzset{addb/.append style={scale=0.7}}
    \tikzset{node distance=0.6}
    % Blocs
    \node[block] (G) {$G_0^\prime$};

    \node[addb, left= of G] (addi) {};
    \node[block, above left=0.3 and 0.3 of addi] (deltai) {$\Delta_I$};
    \node[block, left= of deltai] (wi) {$w_I$};
    \node[branch] (branch) at ($(wi.west|-addi)+(-0.4, 0)$) {};

    % Connections and labels
    \draw[->] (branch.center) |- (wi.west);
    \draw[->] ($(branch)+(-0.6, 0)$) -- (addi.west);
    \draw[->] (wi.east) -- (deltai.west);
    \draw[->] (deltai.east) -| (addi.north);
    \draw[->] (addi.east) -- (G.west);
    \draw[->] (G.east) -- ++(0.6, 0);

    \begin{scope}[on background layer]
      \node[fit={(branch|-wi.north) (G.south east)}, inner sep=6pt, draw, dashed, fill=black!20!white] (Gp) {};
        \node[below left] at (Gp.north east) {$G^\prime$};
    \end{scope}
  \end{tikzpicture}
#+end_src

#+name: fig:input_uncertainty_set
#+caption: Input Multiplicative Uncertainty
#+RESULTS:
[[file:figs/input_uncertainty_set.png]]

We choose a simple uncertainty weight:
\[ w_I(s) = \frac{\tau s + r_0}{(\tau/r_\infty) s + 1} \]
where $r_0$ is the relative uncertainty at steady-state, $1/\tau$ is the frequency at which the relative uncertainty reaches $100\ \%$, and $r_\infty$ is the magnitude of the weight at high frequency.

The parameters are defined below.
#+begin_src matlab
  r0 = 0.5;
  tau = 1/(100*2*pi);
  rinf = 10;

  wI = (tau*s + r0)/((tau/rinf)*s + 1);
#+end_src

We then generate a complex $\Delta$.
#+begin_src matlab
  DeltaI = ucomplex('A',0);
#+end_src

We generate the uncertain plant $G^\prime(s)$.
#+begin_src matlab
  Gp = Gp0*(1+wI*DeltaI);
#+end_src

A set of uncertainty support's compliance transfer functions is shown in Figure [[fig:compliance_support_uncertainty]].

#+begin_src matlab :exports none
  Gps = usample(Gp, 50);

  freqs  = logspace(-1, 4, 1000);

  figure;

  hold on;
  for i = 1:length(Gps)
    plot(freqs, abs(squeeze(freqresp(Gps(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  plot(freqs, abs(squeeze(freqresp(Gp0, freqs, 'Hz'))), 'r-');
  plot(freqs, abs(squeeze(freqresp(Gp0, freqs, 'Hz'))).*(1 + abs(squeeze(freqresp(wI, freqs, 'Hz')))), 'r--');
  plot(freqs, abs(squeeze(freqresp(Gp0, freqs, 'Hz'))).*max(0, (1 - abs(squeeze(freqresp(wI, freqs, 'Hz'))))), 'r--');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]');
  ylabel('Magnitude [dB]');
  hold off;
#+end_src

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

#+name: fig:compliance_support_uncertainty
#+caption: Uncertainty of the support's compliance ([[./figs/compliance_support_uncertainty.png][png]], [[./figs/compliance_support_uncertainty.pdf][pdf]])
[[file:figs/compliance_support_uncertainty.png]]

** Equivalent Inverse Multiplicative Uncertainty
Let's express the uncertainty of the plant $x/F$ as a function of the parameters as well as of the uncertainty on the platform's compliance:
\begin{align*}
  \frac{x}{F} &= \frac{1}{ms^2 + cs + k + ms^2(cs + k)G_0^\prime(s)(1 + w_I(s)\Delta(s))}\\
              &= \frac{1}{ms^2 + cs + k + ms^2(cs + k)G_0^\prime(s) + ms^2(cs + k)G_0^\prime(s) w_I(s)\Delta(s)}\\
              &= \frac{1}{ms^2 + cs + k + ms^2(cs + k)G_0^\prime(s)} \cdot \frac{1}{1 + \frac{ms^2(cs + k)G_0^\prime(s) w_I(s)}{ms^2 + cs + k + ms^2(cs + k)G_0^\prime(s)} \Delta(s)}\\
\end{align*}

#+begin_important
We can the plant dynamics that as an inverse multiplicative uncertainty (Figure [[fig:inverse_uncertainty_set]]):
\begin{equation}
  \frac{x}{F} = G_0(s) (1 + w_{iI}(s) \Delta(s))^{-1}
\end{equation}
with:
- $G_0(s) = \frac{1}{ms^2 + cs + k + ms^2(cs + k)G_0^\prime(s)}$
- $w_{iI}(s) = \frac{ms^2(cs + k)G_0^\prime(s) w_I(s)}{ms^2 + cs + k + ms^2(cs + k)G_0^\prime(s)} = G_0(s) ms^2(cs + k)G_0^\prime(s) w_I(s)$
#+end_important

#+begin_src latex :file inverse_uncertainty_set.pdf
  \begin{tikzpicture}
    \tikzset{block/.default={0.8cm}{0.8cm}}
    \tikzset{addb/.append style={scale=0.7}}
    \tikzset{node distance=0.6}

    % Blocs
    \node[block] (G) {$G$};

    \node[branch, left=0.4 of G] (branch) {};
    \node[block, above left=0.6 and 0.3 of branch] (deltai) {$\Delta_{iI}$};
    \node[block, left= of deltai] (wi) {$w_{iI}$};
    \node[addb] (addu) at ($(wi.west|-G)+(-0.2, 0)$) {};

    % Connections and labels
    \draw[->] (addu.east) -- (G.west);
    \draw[<-] (addu.north) |- (wi.west);
    \draw[<-] (wi.east) node[above right]{$u_\Delta$} -- (deltai.west);
    \draw[<-] (deltai.east) node[above right]{$y_\Delta$} -| (branch.center);
    \draw[->] (G.east) -- ++(0.8, 0);
    \draw[<-] (addu.west) -- ++(-0.8, 0);

      \begin{scope}[on background layer]
        \node[fit={(wi.north-|addu.west) (G.south east)}, inner sep=6pt, draw, dashed, fill=black!20!white] (Gp) {};
        \node[below left] at (Gp.north east) {$G_p$};
      \end{scope}
  \end{tikzpicture}
#+end_src

#+name: fig:inverse_uncertainty_set
#+caption: Inverse Multiplicative Uncertainty
#+RESULTS:
[[file:figs/inverse_uncertainty_set.png]]

** Effect of the Isolation platform Stiffness
Let's first fix the mass of the payload to be isolated:
#+begin_src matlab
  m = 100;
#+end_src

And we generate three isolation platforms:
- A soft one with $\omega_0 = 5\ Hz$
- A medium stiff one with $\omega_0 = 50\ Hz$
- A stiff one with $\omega_0 = 500\ Hz$

Soft Isolation Platform:
#+begin_src matlab
  k_soft = m*(2*pi*5)^2;
  c_soft = 0.1*sqrt(m*k_soft);

  G_soft = 1/(m*s^2 + c_soft*s + k_soft + m*s^2*(c_soft*s + k_soft)*Gp);
  G0_soft = 1/(m*s^2 + c_soft*s + k_soft + m*s^2*(c_soft*s + k_soft)*Gp0);
  wiI_soft = Gp0*m*s^2*(c_soft*s + k_soft)*G0_soft*wI;
#+end_src

Mid Isolation Platform
#+begin_src matlab
  k_mid = m*(2*pi*50)^2;
  c_mid = 0.1*sqrt(m*k_mid);

  G_mid = 1/(m*s^2 + c_mid*s + k_mid + m*s^2*(c_mid*s + k_mid)*Gp);
  G0_mid = 1/(m*s^2 + c_mid*s + k_mid + m*s^2*(c_mid*s + k_mid)*Gp0);
  wiI_mid = Gp0*m*s^2*(c_mid*s + k_mid)*G0_mid*wI;
#+end_src

Stiff Isolation Platform
#+begin_src matlab
  k_stiff = m*(2*pi*500)^2;
  c_stiff = 0.1*sqrt(m*k_stiff);

  G_stiff = 1/(m*s^2 + c_stiff*s + k_stiff + m*s^2*(c_stiff*s + k_stiff)*Gp);
  G0_stiff = 1/(m*s^2 + c_stiff*s + k_stiff + m*s^2*(c_stiff*s + k_stiff)*Gp0);
  wiI_stiff = Gp0*m*s^2*(c_stiff*s + k_stiff)*G0_stiff*wI;
#+end_src

The obtained transfer functions $x/F$ for each of the three platforms are shown in Figure [[fig:plant_uncertainty_stiffness_isolator]].

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

  Gs_soft  = usample(G_soft, 10);
  Gs_mid   = usample(G_mid, 10);
  Gs_stiff = usample(G_stiff, 10);

  figure;

  ax1 = subplot(2,3,1);
  hold on;
  for i = 1:length(Gs_soft)
    plot(freqs, abs(squeeze(freqresp(Gs_soft(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  plot(freqs, abs(squeeze(freqresp(G0_soft, freqs, 'Hz'))), 'r-');
  plot(freqs, abs(squeeze(freqresp(G0_soft, freqs, 'Hz')))./(1 + abs(squeeze(freqresp(wiI_soft, freqs, 'Hz')))), 'r--');
  plot(freqs, abs(squeeze(freqresp(G0_soft, freqs, 'Hz')))./max(0, (1 - abs(squeeze(freqresp(wiI_soft, freqs, 'Hz'))))), 'r--');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  ylabel('Magnitude [dB]');
  hold off;
  title('$\omega_0 \ll \omega_0^\prime$');

  ax4 = subplot(2,3,4);
  hold on;
  for i = 1:length(Gs_soft)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gs_soft(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G0_soft, freqs, 'Hz'))), 'r-');
  Dphi = 180/pi*asin(abs(squeeze(freqresp(wiI_soft, freqs, 'Hz'))));
  Dphi(find(abs(squeeze(freqresp(wiI_soft, freqs, 'Hz'))) > 1, 1):end) = 360;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G0_soft, freqs, 'Hz')))+Dphi, 'r--');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G0_soft, freqs, 'Hz')))-Dphi, 'r--');
  set(gca,'xscale','log');
  yticks(-360:90:180);
  ylim([-180 180]);
  ylabel('Phase [deg]');
  hold off;

  ax2 = subplot(2,3,2);
  hold on;
  for i = 1:length(Gs_mid)
    plot(freqs, abs(squeeze(freqresp(Gs_mid(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  plot(freqs, abs(squeeze(freqresp(G0_mid, freqs, 'Hz'))), 'r-');
  plot(freqs, abs(squeeze(freqresp(G0_mid, freqs, 'Hz')))./(1 + abs(squeeze(freqresp(wiI_mid, freqs, 'Hz')))), 'r--');
  plot(freqs, abs(squeeze(freqresp(G0_mid, freqs, 'Hz')))./max(0, (1 - abs(squeeze(freqresp(wiI_mid, freqs, 'Hz'))))), 'r--');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  hold off;
  title('$\omega_0 \approx \omega_0^\prime$');

  ax5 = subplot(2,3,5);
  hold on;
  for i = 1:length(Gs_mid)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gs_mid(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G0_mid, freqs, 'Hz'))), 'r-');
  Dphi = 180/pi*asin(abs(squeeze(freqresp(wiI_mid, freqs, 'Hz'))));
  Dphi(find(abs(squeeze(freqresp(wiI_mid, freqs, 'Hz'))) > 1, 1):end) = 360;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G0_mid, freqs, 'Hz')))+Dphi, 'r--');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G0_mid, freqs, 'Hz')))-Dphi, 'r--');
  set(gca,'xscale','log');
  yticks(-360:90:180);
  ylim([-180 180]);
  xlabel('Frequency [Hz]');
  hold off;

  ax3 = subplot(2,3,3);
  hold on;
  for i = 1:length(Gs_stiff)
    plot(freqs, abs(squeeze(freqresp(Gs_stiff(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  plot(freqs, abs(squeeze(freqresp(G0_stiff, freqs, 'Hz'))), 'r-');
  plot(freqs, abs(squeeze(freqresp(G0_stiff, freqs, 'Hz')))./(1 + abs(squeeze(freqresp(wiI_stiff, freqs, 'Hz')))), 'r--');
  plot(freqs, abs(squeeze(freqresp(G0_stiff, freqs, 'Hz')))./max(0, (1 - abs(squeeze(freqresp(wiI_stiff, freqs, 'Hz'))))), 'r--');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  hold off;
  title('$\omega_0 \gg \omega_0^\prime$');

  ax6 = subplot(2,3,6);
  hold on;
  for i = 1:length(Gs_stiff)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gs_stiff(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G0_stiff, freqs, 'Hz'))), 'r-');
  Dphi = 180/pi*asin(abs(squeeze(freqresp(wiI_stiff, freqs, 'Hz'))));
  Dphi(find(abs(squeeze(freqresp(wiI_stiff, freqs, 'Hz'))) > 1, 1):end) = 360;
  plot(freqs, 180/pi*angle(squeeze(freqresp(G0_stiff, freqs, 'Hz')))+Dphi, 'r--');
  plot(freqs, 180/pi*angle(squeeze(freqresp(G0_stiff, freqs, 'Hz')))-Dphi, 'r--');
  set(gca,'xscale','log');
  yticks(-360:90:180);
  ylim([-180 180]);
  hold off;

  linkaxes([ax1,ax2,ax3,ax4,ax5,ax6],'x');
  xlim([freqs(1), freqs(end)]);

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

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

#+name: fig:plant_uncertainty_stiffness_isolator
#+caption: Obtained plant for the three isolators ([[./figs/plant_uncertainty_stiffness_isolator.png][png]], [[./figs/plant_uncertainty_stiffness_isolator.pdf][pdf]])
[[file:figs/plant_uncertainty_stiffness_isolator.png]]

The obtain result is very similar to the one obtain in section [[sec:introductory_example]], except for the stiff isolation that experience lot's of uncertainty at high frequency.
This is due to the fact that with the current model, at high frequency, the support's compliance uncertainty is much higher than the previous model.

** Reduce the Uncertainty on the plant
*** Introduction                                                    :ignore:
Now that we know the expression of the uncertainty on the plant, we can wonder what parameters of the isolation platform would lower the plant uncertainty, or at least bring the uncertainty to reasonable level.

The uncertainty of the plant is described by an inverse multiplicative uncertainty with the following weight:
\[ w_{iI}(s) = \frac{ms^2(cs + k)G_0^\prime(s) w_I(s)}{ms^2 + cs + k + ms^2(cs + k)G_0^\prime(s)} \]

Let's study separately the effect of the platform's mass, damping and stiffness.

*** Effect of the platform's stiffness $k$
Let's fix $\xi = \frac{c}{2\sqrt{km}} = 0.1$, $m = 100\ [kg]$ and see the evolution of $|w_{iI}(j\omega)|$ with $k$.

This is first shown for few values of the stiffness $k$ in figure [[fig:inverse_multiplicative_uncertainty_norm_few_k]]

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

  figure;
  hold on;
  for k = logspace(3,9,7)
    c = 0.1*sqrt(m*k);
    G0 = 1/(m*s^2 + c*s + k + m*s^2*(c*s + k)*Gp0);
    plot(freqs, abs(squeeze(freqresp(G0*m*s^2*(c*s + k)*Gp0*wI, freqs, 'Hz'))), 'DisplayName', sprintf('$k = %.0e$', k))
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Uncertainty');
  hold off;
  legend('location', 'southeast');
#+end_src

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

#+name: fig:inverse_multiplicative_uncertainty_norm_few_k
#+caption: caption ([[./figs/inverse_multiplicative_uncertainty_norm_few_k.png][png]], [[./figs/inverse_multiplicative_uncertainty_norm_few_k.pdf][pdf]])
[[file:figs/inverse_multiplicative_uncertainty_norm_few_k.png]]

The norm of the uncertainty weight $|w_iI(j\omega)|$ is displayed as a function of $\omega$ and $k$ in Figure [[fig:inverse_multiplicative_uncertainty_norm_k]].

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

  m = 100;

  Ks = logspace(3, 9, 100);

  wiI_k = zeros(length(freqs), length(Ks));

  for i = 1:length(Ks)
    k = Ks(i);
    c = 0.1*sqrt(m*k);

    G0 = 1/(m*s^2 + c*s + k + m*s^2*(c*s + k)*Gp0);

    wiI_k(:, i) = abs(squeeze(freqresp(G0*m*s^2*(c*s + k)*Gp0*wI, freqs, 'Hz')));
  end
#+end_src

#+begin_src matlab :exports none
  figure;
  surf(freqs, Ks, wiI_k', 'FaceColor', 'interp', 'EdgeColor', 'none')
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); set(gca, 'ZScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Platform Stiffness [N/m]'); zlabel('$|w_{iI}(j\omega)|$');
  xlim([freqs(1), freqs(end)]); ylim([Ks(1), Ks(end)]);
  view([0 0 1]);
  colorbar('location', 'west');
  set(gca,'ColorScale','log')
  caxis([1e-3, 1]);
#+end_src

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

#+name: fig:inverse_multiplicative_uncertainty_norm_k
#+caption: Evolution of the norm of the uncertainty weight $|w_{iI}(j\omega)|$ as a function of the platform's stiffness $k$ ([[./figs/inverse_multiplicative_uncertainty_norm_k.png][png]], [[./figs/inverse_multiplicative_uncertainty_norm_k.pdf][pdf]])
[[file:figs/inverse_multiplicative_uncertainty_norm_k.png]]

Instead of plotting as a function of the platform's stiffness, we can plot as a function of $\omega_0/\omega_0^\prime$ where:
- $\omega_0$ is the resonance of the platform alone
- $\omega_0^\prime$ is the resonance of the support alone

The obtain plot is shown in Figure [[fig:inverse_multiplicative_uncertainty_k_normalized_frequency]].
In that case, we can see that with a platform's resonance frequency 10 times lower than the resonance of the support, we get less than $1\%$ uncertainty.

#+begin_src matlab :exports none
  Ws = sqrt(Ks./m);
  Wn = Ws./sqrt(k0/m0); % Normalized Frequency

  figure;
  contour(freqs, Wn, wiI_k', [0.001, 0.01, 0.1, 1], 'LineWidth', 2, 'ShowText', 'on')
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Normalized Frequency $\frac{\omega_0}{\omega_0^\prime}$');
#+end_src

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

#+name: fig:inverse_multiplicative_uncertainty_k_normalized_frequency
#+caption: Evolution of the norm of the uncertainty weight $|w_{iI}(j\omega)|$ as a function of the frequency ratio $\omega_0/\omega_0^\prime$ ([[./figs/inverse_multiplicative_uncertainty_k_normalized_frequency.png][png]], [[./figs/inverse_multiplicative_uncertainty_k_normalized_frequency.pdf][pdf]])
[[file:figs/inverse_multiplicative_uncertainty_k_normalized_frequency.png]]

*** Effect of the platform's damping $c$
Let's fix $k = 10^7\ [N/m]$, $m = 100\ [kg]$ and see the evolution of $|w_{iI}(j\omega)|$ with the isolator damping $c$ (Figure [[fig:inverse_multiplicative_uncertainty_norm_c]]).

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

  m = 100;
  xi = logspace(-2, 1, 100);

  wiI_c_soft  = zeros(length(freqs), length(xi));
  wiI_c_mid   = zeros(length(freqs), length(xi));
  wiI_c_stiff = zeros(length(freqs), length(xi));
  for i = 1:length(xi)
    k = 1e3;
    c = 2*xi(i)*sqrt(m*k);
    G0 = 1/(m*s^2 + c*s + k + m*s^2*(c*s + k)*Gp0);
    wiI_c_soft(:, i) = abs(squeeze(freqresp(G0*m*s^2*(c*s + k)*Gp0*wI, freqs, 'Hz')));

    k = 1e7;
    c = 2*xi(i)*sqrt(m*k);
    G0 = 1/(m*s^2 + c*s + k + m*s^2*(c*s + k)*Gp0);
    wiI_c_mid(:, i) = abs(squeeze(freqresp(G0*m*s^2*(c*s + k)*Gp0*wI, freqs, 'Hz')));

    k = 1e9;
    c = 2*xi(i)*sqrt(m*k);
    G0 = 1/(m*s^2 + c*s + k + m*s^2*(c*s + k)*Gp0);
    wiI_c_stiff(:, i) = abs(squeeze(freqresp(G0*m*s^2*(c*s + k)*Gp0*wI, freqs, 'Hz')));
  end
#+end_src

#+begin_src matlab :exports none
  figure;

  ax1 = subplot(1, 3, 1);
  surf(freqs, xi, wiI_c_soft', 'FaceColor', 'interp', 'EdgeColor', 'none')
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Damping Ratio');
  title('$\omega_0 \ll \omega_0^\prime$');
  view([0 0 1]);
  set(gca,'ColorScale','log')
  colorbar('location', 'west');
  caxis([1e-3, 1]);

  ax2 = subplot(1, 3, 2);
  surf(freqs, xi, wiI_c_mid', 'FaceColor', 'interp', 'EdgeColor', 'none')
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]');
  title('$\omega_0 \approx \omega_0^\prime$');
  view([0 0 1]);
  set(gca,'ColorScale','log')
  caxis([1e-3, 1]);

  ax3 = subplot(1, 3, 3);
  surf(freqs, xi, wiI_c_stiff', 'FaceColor', 'interp', 'EdgeColor', 'none')
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  title('$\omega_0 \gg \omega_0^\prime$');
  view([0 0 1]);
  set(gca,'ColorScale','log')
  caxis([1e-3 1e0]);
#+end_src

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

#+name: fig:inverse_multiplicative_uncertainty_norm_c
#+caption: Evolution of the norm of the uncertainty weight $|w_{iI}(j\omega)|$ as a function of the platform's damping ratio $\xi$ ([[./figs/inverse_multiplicative_uncertainty_norm_c.png][png]], [[./figs/inverse_multiplicative_uncertainty_norm_c.pdf][pdf]])
[[file:figs/inverse_multiplicative_uncertainty_norm_c.png]]


*** Effect of the platform's mass $m$
Let's fix $k = 10^7\ [N/m]$, $\xi = \frac{c}{2\sqrt{km}} = 0.1$ and see the evolution of $|w_{iI}(j\omega)|$ with the payload mass $m$ (Figure [[fig:inverse_multiplicative_uncertainty_norm_m]]).

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

  k = 1e7;

  Ms = logspace(0, 2, 100);

  wiI_m = zeros(length(freqs), length(Ms));

  for i = 1:length(Ms)
    m = Ms(i);
    c = 2*0.01*sqrt(m*k);

    G0 = 1/(m*s^2 + c*s + k + m*s^2*(c*s + k)*Gp0);

    wiI_m(:, i) = abs(squeeze(freqresp(G0*m*s^2*(c*s + k)*Gp0*wI, freqs, 'Hz')));
  end
#+end_src

#+begin_src matlab :exports none
  figure;
  surf(freqs, Ms, wiI_m', 'FaceColor', 'interp', 'EdgeColor', 'none')
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); set(gca, 'ZScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Payload Mass [kg]'); zlabel('$|w_{iI}(j\omega)|$');
  view([0 0 1]);
  colorbar('location', 'west');
  set(gca,'ColorScale','log')
  caxis([1e-3, 1]);
#+end_src

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

#+name: fig:inverse_multiplicative_uncertainty_norm_m
#+caption: Evolution of the norm of the uncertainty weight $|w_{iI}(j\omega)|$ as a function of the payload mass $m$ ([[./figs/inverse_multiplicative_uncertainty_norm_m.png][png]], [[./figs/inverse_multiplicative_uncertainty_norm_m.pdf][pdf]])
[[file:figs/inverse_multiplicative_uncertainty_norm_m.png]]

** Conclusion
#+begin_important
If the goal is to have an acceptable ($<10\%$) uncertainty on the plant until the highest frequency, two design choice for the isolation platform are possible:
- a very soft isolation platform $\omega_0 \ll \omega_0^\prime$
- a very stiff isolation platform $\omega_0 \gg \omega_0^\prime$

If a very soft isolation platform is used, the uncertainty due to the support's compliance is filtered out and never reaches problematic values.

If a very stiff isolation platform is used, the uncertainty will be high around $\omega_0^\prime$ and may reach unacceptable value.
It will then be high around $\omega_0$ and probably be higher than one.
Thus, if a stiff isolation platform is used, the recommendation is to have the largest possible resonance frequency, as the control bandwidth will be limited by the first resonance of the isolation platform (if not already limited by the resonance of the support).
#+end_important

* Numerical Analysis for the NASS                                   :noexport:
<<sec:nass>>
** 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

** Micro-Station compliance

* Active Damping and Uncertainty                                    :noexport:
** Introduction                                                     :ignore:
The goal of this part is to study how plant uncertainty evolves when active damping is applied.

It seems that as feedback generally reduces uncertainty on the plant, active damping should reduce the uncertainty of the plant.

This may be one of the reason HAC/LAC control architecture is quite effective.

This is also stated in cite:skogestad07_multiv_feedb_contr that cascade control using extra measurements are used to provide local disturbance rejection, linearization or to reduce the effect of measurement noise.

** 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
  freqs = logspace(0, 3, 1000);
#+end_src

** Plant
Let's consider a simple mass spring damper system with the following parameters.
#+begin_src matlab
  m = 10; % [kg]
  c = 100; % [N/(m/s)]
  k = 1e6; % [N/m]
#+end_src

#+name: fig:mech_sys_1dof
#+caption: 1 degree-of-freedom mass-spring-damper system
[[file:figs/mech_sys_1dof.png]]

The nominal system is thus:
#+begin_src matlab
  G = 1/(m*s^2 + c*s + k);
#+end_src

Now, let's consider some uncertainty in the parameters:
#+begin_src matlab
  um = ureal('m', m, 'Percentage', 5);
  uc = ureal('c', c, 'Percentage', 20);
  uk = ureal('k', k, 'Percentage', 30);
#+end_src

And the uncertain model is defined as:
#+begin_src matlab
  Gu = 1/(um*s^2 + uc*s + uk);
#+end_src

The body plot of some of the elements of the uncertain model is shown in Fig. [[fig:uncertain_system]].

#+begin_src matlab :exports none
  Gus = usample(Gu, 50);

  figure;

  % Magnitude
  ax1 = subplot(2,1,1);
  hold on;
  for i = 1:length(Gus)
    plot(freqs, abs(squeeze(freqresp(Gus(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  plot(freqs, abs(squeeze(freqresp(G, freqs, 'Hz'))), 'k-');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  ylabel('Magnitude [dB]');
  hold off;

  % Phase
  ax2 = subplot(2,1,2);
  hold on;
  for i = 1:length(Gus)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gus(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G, freqs, 'Hz'))), 'k-');
  set(gca,'xscale','log');
  yticks(-360:90:180);
  ylim([-270 90]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;

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

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

#+NAME: fig:uncertain_system
#+CAPTION: Uncertain System - Bode Plot ([[./figs/uncertain_system.png][png]], [[./figs/uncertain_system.pdf][pdf]])
[[file:figs/uncertain_system.png]]

** Active Damping
We know apply Direct Velocity Feedback to the system:
\[ F = -\alpha \dot{x} + f \]

#+name: fig:mech_sys_1dof_inertial_contr
#+caption: Direct Velocity Feedback
[[file:figs/mech_sys_1dof_inertial_contr.png]]

Once the loop is closed, the new dynamics of the system is:
\[ \frac{x}{F} = \frac{1}{m s^2 + (c + \alpha) s + k} \]

We compute the new closed-loop uncertain system. Its bode plot is shown in Fig. [[fig:dvf_uncertain_system]].
#+begin_src matlab
  alpha = 1e3;
#+end_src

#+begin_src matlab
  Gu_dvf = 1/(um*s^2 + (uc + alpha)*s + uk);
  G_dvf = 1/(m*s^2 + (c + alpha)*s + k);

  Gus_dvf = usample(Gu_dvf, 50);
#+end_src

#+begin_src matlab :exports none
  figure;

  % Magnitude
  ax1 = subplot(2,1,1);
  hold on;
  for i = 1:length(Gus_dvf)
    plot(freqs, abs(squeeze(freqresp(Gus_dvf(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  plot(freqs, abs(squeeze(freqresp(G_dvf, freqs, 'Hz'))), 'k-');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  ylabel('Magnitude [dB]');
  hold off;

  % Phase
  ax2 = subplot(2,1,2);
  hold on;
  for i = 1:length(Gus_dvf)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Gus_dvf(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]);
  end
  plot(freqs, 180/pi*angle(squeeze(freqresp(G_dvf, freqs, 'Hz'))), 'k-');
  set(gca,'xscale','log');
  yticks(-360:90:180);
  ylim([-270 90]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;

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

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

#+NAME: fig:dvf_uncertain_system
#+CAPTION: Bode plot of the actively damped uncertain system ([[./figs/dvf_uncertain_system.png][png]], [[./figs/dvf_uncertain_system.pdf][pdf]])
[[file:figs/dvf_uncertain_system.png]]

** Uncertainty of the Actively Damped system
As shown in Fig. [[fig:uncertainty_open_loop_dvf]], the uncertainty of the system near the resonance frequency is lowered when active damping is applied.

However, uncertainty is not reduced at all frequencies.

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp((Gus(:, :, 1) - G)*inv(G), freqs, 'Hz'))), '-', 'DisplayName', 'Open Loop', 'color', [0, 0, 0, 0.2]);
  for i = 2:length(Gus)
    plot(freqs, abs(squeeze(freqresp((Gus(:, :, i) - G)*inv(G), freqs, 'Hz'))), '-', 'HandleVisibility', 'off', 'color', [0, 0, 0, 0.2]);
  end

  plot(freqs, abs(squeeze(freqresp((Gus_dvf(:, :, 1) - G)*inv(G), freqs, 'Hz'))), '-', 'DisplayName', 'DVF', 'color', [0, 0, 1, 0.2]);
  for i = 2:length(Gus_dvf)
    plot(freqs, abs(squeeze(freqresp((Gus_dvf(:, :, i) - G)*inv(G), freqs, 'Hz'))), '-', 'HandleVisibility', 'off', 'color', [0, 0, 1, 0.2]);
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]');
  hold off;
  legend('location', 'northeast');
#+end_src

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

#+NAME: fig:uncertainty_open_loop_dvf
#+CAPTION: Multiplicative Uncertainty for the Open Loop system and the Actively damped system ([[./figs/uncertainty_open_loop_dvf.png][png]], [[./figs/uncertainty_open_loop_dvf.pdf][pdf]])
[[file:figs/uncertainty_open_loop_dvf.png]]