diff --git a/org/uncertainty_stiffness.org b/org/uncertainty_stiffness.org new file mode 100644 index 0000000..c9529ef --- /dev/null +++ b/org/uncertainty_stiffness.org @@ -0,0 +1,1344 @@ +#+TITLE: Uncertainty of Resonant Plant +:DRAWER: +#+STARTUP: overview + +#+LANGUAGE: en +#+EMAIL: dehaeze.thomas@gmail.com +#+AUTHOR: Dehaeze Thomas + +#+HTML_LINK_HOME: ./index.html +#+HTML_LINK_UP: ./index.html + +#+HTML_HEAD: +#+HTML_HEAD: +#+HTML_HEAD: +#+HTML_HEAD: +#+HTML_HEAD: +#+HTML_HEAD: +#+HTML_HEAD: + +#+HTML_MATHJAX: align: center tagside: right font: TeX + +#+PROPERTY: header-args:matlab :session *MATLAB* +#+PROPERTY: header-args:matlab+ :comments org +#+PROPERTY: header-args:matlab+ :results none +#+PROPERTY: header-args:matlab+ :exports both +#+PROPERTY: header-args:matlab+ :eval no-export +#+PROPERTY: header-args:matlab+ :output-dir figs +#+PROPERTY: header-args:matlab+ :tangle no +#+PROPERTY: header-args:matlab+ :mkdirp yes + +#+PROPERTY: header-args:shell :eval no-export + +#+PROPERTY: header-args:latex :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/Cloud/thesis/latex/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: +* 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. + +** Matlab Init :noexport:ignore: +#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) + <> +#+end_src + +#+begin_src matlab :exports none :results silent :noweb yes + <> +#+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") + <> +#+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") + <> +#+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") + <> +#+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]] + +* 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_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) + <> +#+end_src + +#+begin_src matlab :exports none :results silent :noweb yes + <> +#+end_src + +** 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)]$ + +#+begin_src matlab + mpi = 1e3; + cpi = 5e4; + kpi = 1e8; +#+end_src + +#+begin_src matlab + mp = ureal('m', mpi, 'Percentage', 10); + cp = ureal('c', cpi, 'Percentage', 20); + kp = ureal('k', kpi, 'Percentage', 30); +#+end_src + +#+begin_src matlab + bodeFig({1/(mpi*s^2 + cpi*s + kpi)}) +#+end_src + +** Initialization of the isolation platform +The mass of the payload is constant: +#+begin_src matlab + m = 100; +#+end_src + +Soft Isolation Platform: +#+begin_src matlab + 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)); +#+end_src + +Mid Isolation Platform +#+begin_src matlab + 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)); +#+end_src + +Stiff Isolation Platform +#+begin_src matlab + 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)); +#+end_src + + +** Comparison +#+begin_src matlab :exports none + f_soft = logspace(-1, 2, 1000); + f_mid = logspace( 0, 3, 1000); + f_stiff = logspace( 1, 4, 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(f_soft, 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 [dB]'); + hold off; + + ax2 = subplot(2,3,4); + hold on; + for i = 1:length(Gs_soft) + plot(f_soft, 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]); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + hold off; + + linkaxes([ax1,ax2],'x'); + xlim([f_soft(1), f_soft(end)]); + + ax1 = subplot(2,3,2); + hold on; + for i = 1:length(Gs_mid) + plot(f_mid, 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',[]); + ylabel('Magnitude [dB]'); + hold off; + + ax2 = subplot(2,3,5); + hold on; + for i = 1:length(Gs_mid) + plot(f_mid, 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]'); ylabel('Phase [deg]'); + hold off; + + linkaxes([ax1,ax2],'x'); + xlim([f_mid(1), f_mid(end)]); + + ax1 = subplot(2,3,3); + hold on; + for i = 1:length(Gs_stiff) + plot(f_stiff, 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',[]); + ylabel('Magnitude [dB]'); + hold off; + + ax2 = subplot(2,3,6); + hold on; + for i = 1:length(Gs_stiff) + plot(f_stiff, 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]); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + hold off; + + linkaxes([ax1,ax2],'x'); + xlim([f_stiff(1), f_stiff(end)]); +#+end_src + +* 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_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 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[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]] + +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} + +** Matlab Init :noexport:ignore: +#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) + <> +#+end_src + +#+begin_src matlab :exports none :results silent :noweb yes + <> +#+end_src + +** 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$. + +#+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. +#+begin_src matlab + r0 = 0.5; + tau = 1/(50*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; + + ax1 = subplot(2,1,1); + hold on; + plot(freqs, abs(squeeze(freqresp(Gp0, freqs, 'Hz'))), 'k--'); + for i = 1:length(Gps) + plot(freqs, abs(squeeze(freqresp(Gps(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]); + end + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + set(gca, 'XTickLabel',[]); + ylabel('Magnitude [dB]'); + hold off; + + ax2 = subplot(2,1,2); + hold on; + plot(freqs, 180/pi*angle(squeeze(freqresp(Gp0, freqs, 'Hz'))), 'k--'); + for i = 1:length(Gps) + plot(freqs, 180/pi*angle(squeeze(freqresp(Gps(:, :, i), freqs, 'Hz'))), '-', 'color', [0, 0, 0, 0.2]); + end + set(gca,'xscale','log'); + yticks(-360:90:180); + ylim([-180 180]); + xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); + hold off; + + linkaxes([ax1,ax2],'x'); + xlim([freqs(1), freqs(end)]); +#+end_src + +#+header: :tangle no :exports results :results none :noweb yes +#+begin_src matlab :var filepath="figs/compliance_support_uncertainty.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png") +<> +#+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]] + +** 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); +#+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); +#+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); +#+end_src + +The obtained transfer functions $x/F$ for each of the three platforms are shown in Figure [[]]. +#+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 [dB]'); + hold off; + title('Soft Isolator'); + + ax2 = 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([-180 180]); + ylabel('Phase [deg]'); + hold off; + + linkaxes([ax1,ax2],'x'); + xlim([freqs(1), freqs(end)]); + + ax1 = 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',[]); + hold off; + title('Medium Stiff Isolator'); + + ax2 = 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([-180 180]); + xlabel('Frequency [Hz]'); + hold off; + + linkaxes([ax1,ax2],'x'); + xlim([freqs(1), freqs(end)]); + + ax1 = 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',[]); + hold off; + title('Stiff Isolator'); + + ax2 = 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([-180 180]); + hold off; + + linkaxes([ax1,ax2],'x'); + xlim([freqs(1), freqs(end)]); +#+end_src + +#+header: :tangle no :exports results :results none :noweb yes +#+begin_src matlab :var filepath="figs/plant_uncertainty_stiffness_isolator.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png") +<> +#+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]] + +** 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_important +\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]] + +** 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$. + +#+begin_src matlab + 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 + 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') +#+end_src + +#+begin_src matlab + 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]'); +#+end_src + +#+begin_src matlab + 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]'); +#+end_src + +#+begin_src matlab + 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 + +Test to verify the formula +#+begin_src matlab + 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'); +#+end_src + +*** 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$. + +#+begin_src matlab + 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 +#+end_src + +#+begin_src matlab + 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') +#+end_src + +#+begin_src matlab + 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'); +#+end_src + +* Old :noexport: +** Initialization +#+begin_src matlab + m1 = ureal('m', 10, 'Percentage', 5 ); + c1 = ureal('c', 100, 'Percentage', 20); + k1 = ureal('k', 1e6, 'Percentage', 30); +#+end_src + +#+begin_src matlab + G1 = 1/(m1*s^2 + c1*s + k1); +#+end_src + +#+begin_src matlab :exports none + Gus = usample(G1, 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 + 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 + 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/lower_plant_uncertainty.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png") + <> +#+end_src + +#+NAME: fig:lower_plant_uncertainty +#+CAPTION: caption ([[./figs/lower_plant_uncertainty.png][png]], [[./figs/lower_plant_uncertainty.pdf][pdf]]) +[[file:figs/lower_plant_uncertainty.png]] + + +** Soft Plant +Let's consider a simple mass spring damper system with the following parameters. +#+begin_src matlab + m1 = ureal('m', 10, 'Percentage', 5); + c1 = ureal('c', 100, 'Percentage', 20); + k1 = ureal('k', 1e6, 'Percentage', 30); + + m2 = 5; % [kg] + c2 = 100; % [N/(m/s)] + k2 = 1e3; % [N/m] +#+end_src + +#+begin_src matlab :exports none + Mm = diag([m2, m1]); + Km = diag([k2, k1+k2]) - diag([k2], -1) - diag([k2], 1); + Cm = diag([c2, c1+c2]) - diag([c2], -1) - diag([c2], 1); +#+end_src + +#+begin_src matlab :exports none + StateName = {... + 'x2', ... + 'x1', ... + 'v2', ... + 'v1' ... + }; + StateUnit = {'m', 'm', 'm/s', 'm/s'}; + + InputName = {... + 'f' ... % [m] + }; + InputUnit = {'N'}; + + OutputName = {... + 'x2' ... % [m] + }; + OutputUnit = {'m'}; +#+end_src + +#+begin_src matlab :exports none + % A Matrix + A = [zeros(size(Mm)) eye(size(Mm)) ; ... + -Mm\Km -Mm\Cm]; + + % B Matrix + B_low = zeros(length(StateName), length(InputName)); + B_low(strcmp(StateName,'v2'), strcmp(InputName,'f')) = 1; + B_low(strcmp(StateName,'v1'), strcmp(InputName,'f')) = -1; + + B = blkdiag(zeros(length(StateName)/2), diag([1/m2, 1/m1]))*B_low; + + % C Matrix + C = zeros(length(OutputName), length(StateName)); + C(strcmp(OutputName,'x2'), strcmp(StateName,'x2')) = 1; + + % D Matrix + D = zeros(length(OutputName), length(InputName)); +#+end_src + +#+begin_src matlab :exports none + sys_soft = ss(A, B, C, D); + sys_soft.StateName = StateName; + sys_soft.StateUnit = StateUnit; + + sys_soft.InputName = InputName; + sys_soft.InputUnit = InputUnit; + + sys_soft.OutputName = OutputName; + sys_soft.OutputUnit = OutputUnit; +#+end_src + +#+begin_src matlab :exports none + figure; + + % Magnitude + ax1 = subplot(2,1,1); + hold on; + + plot(freqs, abs(squeeze(freqresp(getNominal(sys_soft), freqs, 'Hz'))), '-', 'color', [0 0.4470 0.7410]); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + set(gca, 'XTickLabel',[]); + ylabel('Magnitude [dB]'); + hold off; + + % Phase + ax2 = subplot(2,1,2); + hold on; + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(getNominal(sys_soft), freqs, 'Hz')))), '-', 'color', [0 0.4470 0.7410]); + 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/soft_plant_uncertainty.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png") + <> +#+end_src + +#+NAME: fig:soft_plant_uncertainty +#+CAPTION: caption ([[./figs/soft_plant_uncertainty.png][png]], [[./figs/soft_plant_uncertainty.pdf][pdf]]) +[[file:figs/soft_plant_uncertainty.png]] + + +** Stiff Plant +Let's consider a simple mass spring damper system with the following parameters. +#+begin_src matlab + m2 = 5; % [kg] + c2 = 100; % [N/(m/s)] + k2 = 1e6; % [N/m] +#+end_src + +#+begin_src matlab :exports none + Mm = diag([m2, m1]); + Km = diag([k2, k1+k2]) - diag([k2], -1) - diag([k2], 1); + Cm = diag([c2, c1+c2]) - diag([c2], -1) - diag([c2], 1); +#+end_src + +#+begin_src matlab :exports none + StateName = {... + 'x2', ... + 'x1', ... + 'v2', ... + 'v1' ... + }; + StateUnit = {'m', 'm', 'm/s', 'm/s'}; + + InputName = {... + 'f' ... % [m] + }; + InputUnit = {'N'}; + + OutputName = {... + 'x2' ... % [m] + }; + OutputUnit = {'m'}; +#+end_src + +#+begin_src matlab :exports none + % A Matrix + A = [zeros(size(Mm)) eye(size(Mm)) ; ... + -Mm\Km -Mm\Cm]; + + % B Matrix + B_low = zeros(length(StateName), length(InputName)); + B_low(strcmp(StateName,'v2'), strcmp(InputName,'f')) = 1; + B_low(strcmp(StateName,'v1'), strcmp(InputName,'f')) = -1; + + B = blkdiag(zeros(length(StateName)/2), diag([1/m2, 1/m1]))*B_low; + + % C Matrix + C = zeros(length(OutputName), length(StateName)); + C(strcmp(OutputName,'x2'), strcmp(StateName,'x2')) = 1; + + % D Matrix + D = zeros(length(OutputName), length(InputName)); +#+end_src + +#+begin_src matlab :exports none + sys_stiff = ss(A, B, C, D); + sys_stiff.StateName = StateName; + sys_stiff.StateUnit = StateUnit; + + sys_stiff.InputName = InputName; + sys_stiff.InputUnit = InputUnit; + + sys_stiff.OutputName = OutputName; + sys_stiff.OutputUnit = OutputUnit; +#+end_src + +#+begin_src matlab :exports none + figure; + + % Magnitude + ax1 = subplot(2,1,1); + hold on; + + plot(freqs, abs(squeeze(freqresp(getNominal(sys_stiff), freqs, 'Hz'))), '-', 'color', [0.8500 0.3250 0.0980]); + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + set(gca, 'XTickLabel',[]); + ylabel('Magnitude [dB]'); + hold off; + + % Phase + ax2 = subplot(2,1,2); + hold on; + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(getNominal(sys_stiff), freqs, 'Hz')))), '-', 'color', [0.8500 0.3250 0.0980]); + 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/stiff_plant_uncertainty.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png") + <> +#+end_src + +#+NAME: fig:stiff_plant_uncertainty +#+CAPTION: caption ([[./figs/stiff_plant_uncertainty.png][png]], [[./figs/stiff_plant_uncertainty.pdf][pdf]]) +[[file:figs/stiff_plant_uncertainty.png]] + + +** Comparison +#+begin_src matlab + G_soft = usample(sys_soft('x2', 'f'), 20); + G_stiff = usample(sys_stiff('x2', 'f'), 20); +#+end_src + +#+begin_src matlab :exports none + figure; + + % Magnitude + ax1 = subplot(2,1,1); + hold on; + + plot(freqs, abs(squeeze(freqresp(getNominal(sys_soft), freqs, 'Hz'))), '-', 'color', [0 0.4470 0.7410]); + for i = 1:length(G_soft) + plot(freqs, abs(squeeze(freqresp(G_soft(:, :, i), freqs, 'Hz'))), '-', 'color', [0 0.4470 0.7410, 0.2]); + end + plot(freqs, abs(squeeze(freqresp(getNominal(sys_stiff), freqs, 'Hz'))), '-', 'color', [0.8500 0.3250 0.0980]); + for i = 1:length(G_stiff) + plot(freqs, abs(squeeze(freqresp(G_stiff(:, :, i), freqs, 'Hz'))), '-', 'color', [0.8500 0.3250 0.0980, 0.2]); + end + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + set(gca, 'XTickLabel',[]); + ylabel('Magnitude [dB]'); + hold off; + + % Phase + ax2 = subplot(2,1,2); + hold on; + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(getNominal(sys_soft), freqs, 'Hz')))), '-', 'color', [0 0.4470 0.7410]); + for i = 1:length(G_soft) + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_soft(:, :, i), freqs, 'Hz')))), '-', 'color', [0 0.4470 0.7410, 0.2]); + end + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(getNominal(sys_stiff), freqs, 'Hz')))), '-', 'color', [0.8500 0.3250 0.0980]); + for i = 1:length(G_stiff) + plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_stiff(:, :, i), freqs, 'Hz')))), '-', 'color', [0.8500 0.3250 0.0980, 0.2]); + end + 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/compaire_stiff_soft_plants.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png") + <> +#+end_src + +#+NAME: fig:compaire_stiff_soft_plants +#+CAPTION: caption ([[./figs/compaire_stiff_soft_plants.png][png]], [[./figs/compaire_stiff_soft_plants.pdf][pdf]]) +[[file:figs/compaire_stiff_soft_plants.png]] + + +#+begin_src matlab :exports none + figure; + hold on; + for i = 1:length(G_soft) + plot(freqs, abs(squeeze(freqresp(G_soft(:, :, i), freqs, 'Hz')-freqresp(getNominal(sys_soft), freqs, 'Hz')))./abs(squeeze(freqresp(getNominal(sys_soft), freqs, 'Hz'))), '-', 'color', [0 0.4470 0.7410, 0.2]); + end + for i = 1:length(G_stiff) + plot(freqs, abs(squeeze(freqresp(G_stiff(:, :, i), freqs, 'Hz')-freqresp(getNominal(sys_stiff), freqs, 'Hz')))./abs(squeeze(freqresp(getNominal(sys_stiff), freqs, 'Hz'))), '-', 'color', [0.8500 0.3250 0.0980, 0.2]); + end + 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/compare_stiff_soft_uncertainty.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png") + <> +#+end_src + +#+NAME: fig:compare_stiff_soft_uncertainty +#+CAPTION: caption ([[./figs/compare_stiff_soft_uncertainty.png][png]], [[./figs/compare_stiff_soft_uncertainty.pdf][pdf]]) +[[file:figs/compare_stiff_soft_uncertainty.png]]