nass-simscape/org/uncertainty_stiffness.org

Uncertainty of Resonant Plant

• Introduction
• Active Damping and Uncertainty
  • Introduction
  • Plant
  • Active Damping
  • Uncertainty of the Actively Damped system
• Uncertainty and Stiffness
  • Introduction
  • Equations of motion
  • Initialization of the support
  • Initialization of the isolation platform
  • Comparison
• Generalization to arbitrary dynamics
  • Introduction
  • Verification of the obtained formula
  • Compliance of the Support
  • Effect of the Isolation platform Stiffness.
  • Equivalent Inverse Multiplicative Uncertainty
  • Reduce the Uncertainty on the plant
    • Introduction
    • Effect of the stiffness $k$
    • Effect of the damping $c$

Introduction   ignore

Active Damping and Uncertainty

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.

Plant

Let's consider a simple mass spring damper system with the following parameters.

  m = 10; % [kg]
  c = 100; % [N/(m/s)]
  k = 1e6; % [N/m]

1 degree-of-freedom mass-spring-damper system

The nominal system is thus:

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

Now, let's consider some uncertainty in the parameters:

  um = ureal('m', m, 'Percentage', 5);
  uc = ureal('c', c, 'Percentage', 20);
  uk = ureal('k', k, 'Percentage', 30);

And the uncertain model is defined as:

  Gu = 1/(um*s^2 + uc*s + uk);

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

  <<plt-matlab>>

Uncertain System - Bode Plot (png, pdf)

Active Damping

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

Direct Velocity Feedback

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.

  alpha = 1e3;

  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);

  <<plt-matlab>>

Bode plot of the actively damped uncertain system (png, pdf)

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.

  <<plt-matlab>>

Multiplicative Uncertainty for the Open Loop system and the Actively damped system (png, pdf)

Uncertainty and Stiffness

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 $d$ using $F$ in spite of uncertainty on the support mechanical properties.

  \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}

Two degrees-of-freedom system

Equations of motion

If we write the equation of motion of the system in Figure ref: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:

\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

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)]$

  mpi = 1e3;
  cpi = 5e4;
  kpi = 1e8;

  mp = ureal('m', mpi, 'Percentage', 10);
  cp = ureal('c', cpi, 'Percentage', 20);
  kp = ureal('k', kpi, 'Percentage', 30);

  bodeFig({1/(mpi*s^2 + cpi*s + kpi)})

Initialization of the isolation platform

The mass of the payload is constant:

  m = 100;

Soft Isolation Platform:

  k = m*(2*pi*5)^2;
  c = 0.1*sqrt(m*k);

  G_soft = (mp*s^2 + cp*s + kp)/(m*s^2*(c*s + k) + (m*s^2 + c*s + k)*(mp*s^2 + cp*s + kp));

Mid Isolation Platform

  k = m*(2*pi*50)^2;
  c = 0.1*sqrt(m*k);

  G_mid = (mp*s^2 + cp*s + kp)/(m*s^2*(c*s + k) + (m*s^2 + c*s + k)*(mp*s^2 + cp*s + kp));

Stiff Isolation Platform

  k = m*(2*pi*500)^2;
  c = 0.1*sqrt(m*k);

  G_stiff = (mp*s^2 + cp*s + kp)/(m*s^2*(c*s + k) + (m*s^2 + c*s + k)*(mp*s^2 + cp*s + kp));

Comparison

Generalization to 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{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}

General support

Now let's consider the system consisting of a mass-spring-system on top of a general support as shown in Figure fig:general_support_with_isolator.

  \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[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}

Mass-Spring-Damper system on top of a general support

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:

\begin{equation} \frac{x}{F} = \frac{1}{ms^2 + cs + k + ms^2(cs + k)G^\prime(s)} \end{equation}

Verification of the obtained formula

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

\begin{equation} \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)} \end{equation}

Which is the same transfer function that was obtained.

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 them $\omega^\prime = \sqrt{\frac{m^\prime}{k^\prime}} \approx 50\ Hz$.

  m0 = 1e3;
  c0 = 5e4;
  k0 = 1e8;

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

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{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}

Input Multiplicative Uncertainty

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.

  r0 = 0.5;
  tau = 1/(50*2*pi);
  rinf = 10;

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

We then generate a complex $\Delta$.

  DeltaI = ucomplex('A',0);

We generate the uncertain plant $G^\prime(s)$.

  Gp = Gp0*(1+wI*DeltaI);

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

<<plt-matlab>>

Uncertainty of the support's compliance (png, pdf)

Effect of the Isolation platform Stiffness.

Let's first fix the mass of the payload to be isolated:

  m = 100;

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:

  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);

Mid Isolation Platform

  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);

Stiff Isolation Platform

  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);

The obtained transfer functions $x/F$ for each of the three platforms are shown in Figure /tdehaeze/nass-simscape/src/commit/7478d7c968c7d71f5ce488b84929072de57f67d5/org/.

<<plt-matlab>>

Obtained plant for the three isolators (png, pdf)

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

We can rewrite 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)$

  \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}

Inverse Multiplicative Uncertainty

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.

We can do that by reducing: \[ 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)} \]

Effect of the stiffness $k$

Let's fix $\xi = 0.1$, $m = 100\ [kg]$ and see the evolution of $|w_{iI}(j\omega)|$ with $k$.

  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

  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)|$');
  colorbar;
  set(gca,'ColorScale','log')

  figure;
  contour(freqs, Ks, 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('Platform Stiffness $k$ [N/m]');

  Ws = sqrt(Ks./m)/2/pi;

  figure;
  contour(freqs, Ws, 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('Resonance Frequency $\omega_0$ [Hz]');

  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}$');

Test to verify the formula

  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 = %g$', k))
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Uncertainty');
  hold off;
  legend('location', 'southeast');

Effect of the 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$.

  freqs = logspace(0, 3, 1000);

  m = 100;
  k = 1e3;

  xi = logspace(-2, 1, 100)

  wiI_c = zeros(length(freqs), length(xi));

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

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

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

  figure;
  surf(freqs, xi, wiI_c', 'FaceColor', 'interp', 'EdgeColor', 'none')
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); set(gca, 'ZScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Damping Ratio'); zlabel('$|w_{iI}(j\omega)|$');
  colorbar;
  set(gca,'ColorScale','log')

  figure;
  contour(freqs, xi, wiI_c', [0.01, 0.1, 0.5, 1], 'LineWidth', 2, 'ShowText', 'on')
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Damping Ratio');