tdehaeze/nass-simscape
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
7ca627ff91
nass-simscape/disturbances/index.org

764 lines
30 KiB
Org Mode
Raw Blame History

#+TITLE: Identification of the disturbances
: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: <link rel="stylesheet" type="text/css" href="../css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/readtheorg.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/zenburn.css"/>
#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.stickytableheaders.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/readtheorg.js"></script>
#+HTML_MATHJAX: align: center tagside: right font: TeX
#+PROPERTY: header-args:matlab :session *MATLAB*
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :results none
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :output-dir figs
#+PROPERTY: header-args:matlab+ :tangle matlab/modal_frf_coh.m
#+PROPERTY: header-args:matlab+ :mkdirp yes
#+PROPERTY: header-args:shell :eval no-export
#+PROPERTY: header-args:latex :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/Cloud/thesis/latex/}{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 raw replace :buffer no
#+PROPERTY: header-args:latex+ :eval no-export
#+PROPERTY: header-args:latex+ :exports both
#+PROPERTY: header-args:latex+ :mkdirp yes
#+PROPERTY: header-args:latex+ :output-dir figs
:END:
* Introduction :ignore:
The goal here is to extract the Power Spectral Density of the sources of perturbation.
The sources of perturbations are (schematically shown in figure [[fig:uniaxial-model-micro-station]]):
- $D_w$: Ground Motion
- Parasitic forces applied in the system when scanning with the Translation Stage and the Spindle ($F_{rz}$ and $F_{ty}$).
These forces can be due to imperfect guiding for instance.
Because we cannot measure directly the perturbation forces, we have the measure the effect of those perturbations on the system (in terms of velocity for instance using geophones, $D$ on figure [[fig:uniaxial-model-micro-station]]) and then, using a model, compute the forces that induced such velocity.
#+begin_src latex :file uniaxial-model-micro-station.pdf :post pdf2svg(file=*this*, ext="png") :exports results
\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.4} % Width of the dashed line for the displacement
\def\disph{0.3} % Height of the arrow for the displacements
\def\bracs{0.05} % Brace spacing vertically
\def\brach{-12pt} % Brace shift horizontaly
\def\fsensh{0.2} % Height of the force sensor
\def\velsize{0.2} % Size of the velocity sensor
% ====================
% ====================
% Ground
% ====================
\draw (-0.5*\massw, 0) -- (0.5*\massw, 0);
\draw[dashed] (0.5*\massw, 0) -- ++(1, 0);
\draw[->, dasehd] (0.5*\massw+0.8, 0) -- ++(0, 0.5) node[right]{$D_{w}$};
% ====================
% ====================
% Marble
\begin{scope}[shift={(0, 0)}]
% Mass
\draw[fill=white] (-0.5*\massw, \spaceh) rectangle (0.5*\massw, \spaceh+\massh) node[pos=0.5]{$m_{m}$};
% Spring, Damper, and Actuator
\draw[spring] (-0.4*\massw, 0) -- (-0.4*\massw, \spaceh) node[midway, left=0.1]{$k_{m}$};
\draw[damper] (0, 0) -- ( 0, \spaceh) node[midway, left=0.2]{$c_{m}$};
% 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]{Marble};
% Displacements
\draw[dashed] (0.5*\massw, \spaceh+\massh) -- ++(2*\dispw, 0) coordinate(xm) -- ++(2.2*\dispw, 0) coordinate(dbot);
\end{scope}
% ====================
% ====================
% Ty
\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_{t}$};
% Spring, Damper, and Actuator
\draw[spring] (-0.4*\massw, 0) -- (-0.4*\massw, \spaceh) node[midway, left=0.1]{$k_{t}$};
\draw[damper] (0, 0) -- ( 0, \spaceh) node[midway, left=0.2]{$c_{t}$};
\draw[actuator={0.45}{0.2}] ( 0.4*\massw, 0) -- ( 0.4*\massw, \spaceh) node[midway, right=0.1](ft){$F_{ty}$};
% 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]{Ty};
\end{scope}
% ====================
% ====================
% Rz
\begin{scope}[shift={(0, 2*(\spaceh+\massh))}]
% Mass
\draw[fill=white] (-0.5*\massw, \spaceh) rectangle (0.5*\massw, \spaceh+\massh) node[pos=0.5]{$m_{z}$};
% Spring, Damper, and Actuator
\draw[spring] (-0.4*\massw, 0) -- (-0.4*\massw, \spaceh) node[midway, left=0.1]{$k_{z}$};
\draw[damper] (0, 0) -- ( 0, \spaceh) node[midway, left=0.2]{$c_{z}$};
\draw[actuator] ( 0.4*\massw, 0) -- ( 0.4*\massw, \spaceh) node[midway, right=0.1](F){$F_{rz}$};
% 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]{Rz};
\end{scope}
% ====================
% ====================
% Hexapod
\begin{scope}[shift={(0, 3*(\spaceh+\massh))}]
% Mass
\draw[fill=white] (-0.5*\massw, \spaceh) rectangle (0.5*\massw, \spaceh+\massh) node[pos=0.5]{$m_{h}$};
% Spring, Damper, and Actuator
\draw[spring] (-0.4*\massw, 0) -- (-0.4*\massw, \spaceh) node[midway, left=0.1]{$k_{h}$};
\draw[damper] (0, 0) -- ( 0, \spaceh) node[midway, left=0.2]{$c_{h}$};
% Displacements
\draw[dashed] (0.5*\massw, \spaceh+\massh) -- ++(2*\dispw, 0) coordinate(xs) -- ++(2.2*\dispw, 0) coordinate(dtop);
% 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]{Hexapod};
\end{scope}
% ====================
\draw[<->] ($(dbot)+(-0.3,0)$) --node[midway, right]{$D$} ($(dtop)+(-0.3,0)$);
\end{tikzpicture}
#+end_src
#+name: fig:uniaxial-model-micro-station
#+caption: Schematic of the Micro Station and the sources of disturbance
#+RESULTS:
[[file:figs/uniaxial-model-micro-station.png]]
This file is divided in the following sections:
- Section [[sec:identification]]: transfer functions from the disturbance forces to the relative velocity of the hexapod with respect to the granite are computed using the Simscape Model representing the experimental setup
- Section [[sec:sensitivity_disturbances]]: the bode plot of those transfer functions are shown
- Section [[sec:psd_dist]]: the measured PSD of the effect of the disturbances are shown
- Section [[sec:psd_force_dist]]: from the model and the measured PSD, the PSD of the disturbance forces are computed
- Section [[sec:noise_budget]]: with the computed PSD, the noise budget of the system is done
* Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
<<matlab-dir>>
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+end_src
#+begin_src matlab :tangle no
simulinkproject('../');
#+end_src
#+begin_src matlab
open 'simscape/sim_micro_station_disturbances.slx'
#+end_src
* Identification
<<sec:identification>>
The transfer functions from the disturbance forces to the relative velocity of the hexapod with respect to the granite are computed using the Simscape Model representing the experimental setup with the code below.
#+begin_src matlab
%% Options for Linearized
options = linearizeOptions;
options.SampleTime = 0;
%% Name of the Simulink File
mdl = 'sim_micro_station_disturbances';
#+end_src
#+begin_src matlab
%% Micro-Hexapod
% Input/Output definition
io(1) = linio([mdl, '/Dw'], 1, 'input'); % Ground Motion
io(2) = linio([mdl, '/Fty'], 1, 'input'); % Parasitic force Ty
io(3) = linio([mdl, '/Frz'], 1, 'input'); % Parasitic force Rz
io(4) = linio([mdl, '/Dgm'], 1, 'output'); % Absolute motion - Granite
io(5) = linio([mdl, '/Dhm'], 1, 'output'); % Absolute Motion - Hexapod
io(6) = linio([mdl, '/Vm'], 1, 'output'); % Relative Velocity hexapod/granite
#+end_src
#+begin_src matlab
% Run the linearization
G = linearize(mdl, io, 0);
% Input/Output names
G.InputName = {'Dw', 'Fty', 'Frz'};
G.OutputName = {'Dgm', 'Dhm', 'Vm'};
#+end_src
* Sensitivity to Disturbances
<<sec:sensitivity_disturbances>>
#+begin_src matlab :exports none
freqs = logspace(0, 3, 1000);
figure;
title('$D_w$ to $D$');
hold on;
plot(freqs, abs(squeeze(freqresp(G('Vm', 'Dw')/s, freqs, 'Hz'))), 'k-');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude [m/m]'); xlabel('Frequency [Hz]');
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_dist_gm.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:sensitivity_dist_gm
#+CAPTION: Sensitivity to Ground Motion ([[./figs/sensitivity_dist_gm.png][png]], [[./figs/sensitivity_dist_gm.pdf][pdf]])
[[file:figs/sensitivity_dist_gm.png]]
#+begin_src matlab :exports none
freqs = logspace(0, 3, 1000);
figure;
title('$F_{ty}$ to $D$');
hold on;
plot(freqs, abs(squeeze(freqresp(G('Vm', 'Fty')/s, freqs, 'Hz'))), 'k-');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_dist_fty.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:sensitivity_dist_fty
#+CAPTION: Sensitivity to vertical forces applied by the Ty stage ([[./figs/sensitivity_dist_fty.png][png]], [[./figs/sensitivity_dist_fty.pdf][pdf]])
[[file:figs/sensitivity_dist_fty.png]]
#+begin_src matlab :exports none
freqs = logspace(0, 3, 1000);
figure;
title('$F_{rz}$ to $D$');
hold on;
plot(freqs, abs(squeeze(freqresp(G('Vm', 'Frz')/s, freqs, 'Hz'))), 'k-');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/sensitivity_dist_frz.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:sensitivity_dist_frz
#+CAPTION: Sensitivity to vertical forces applied by the Rz stage ([[./figs/sensitivity_dist_frz.png][png]], [[./figs/sensitivity_dist_frz.pdf][pdf]])
[[file:figs/sensitivity_dist_frz.png]]
* Power Spectral Density of the effect of the disturbances
<<sec:psd_dist>>
The PSD of the relative velocity between the hexapod and the marble in $[(m/s)^2/Hz]$ are loaded for the following sources of disturbance:
- Slip Ring Rotation
- Scan of the translation stage (effect in the vertical direction and in the horizontal direction)
Also, the Ground Motion is measured.
#+begin_src matlab
gm = load('./disturbances/mat/psd_gm.mat', 'f', 'psd_gm', 'psd_gv');
rz = load('./disturbances/mat/pxsp_r.mat', 'f', 'pxsp_r');
tyz = load('./disturbances/mat/pxz_ty_r.mat', 'f', 'pxz_ty_r');
tyx = load('./disturbances/mat/pxe_ty_r.mat', 'f', 'pxe_ty_r');
#+end_src
#+begin_src matlab :exports none
gm.f = gm.f(2:end);
rz.f = rz.f(2:end);
tyz.f = tyz.f(2:end);
tyx.f = tyx.f(2:end);
#+end_src
#+begin_src matlab :exports none
gm.psd_gm = gm.psd_gm(2:end); % PSD of Ground Motion [m^2/Hz]
gm.psd_gv = gm.psd_gv(2:end); % PSD of Ground Velocity [(m/s)^2/Hz]
rz.pxsp_r = rz.pxsp_r(2:end); % PSD of Relative Velocity [(m/s)^2/Hz]
tyz.pxz_ty_r = tyz.pxz_ty_r(2:end); % PSD of Relative Velocity [(m/s)^2/Hz]
tyx.pxe_ty_r = tyx.pxe_ty_r(2:end); % PSD of Relative Velocity [(m/s)^2/Hz]
#+end_src
We now compute the relative velocity between the hexapod and the granite due to ground motion.
#+begin_src matlab
gm.psd_rv = gm.psd_gm.*abs(squeeze(freqresp(G('Vm', 'Dw'), gm.f, 'Hz'))).^2;
#+end_src
The Power Spectral Density of the relative motion/velocity of the hexapod with respect to the granite are shown in figures [[fig:dist_effect_relative_velocity]] and [[fig:dist_effect_relative_motion]].
The Cumulative Amplitude Spectrum of the relative motion is shown in figure [[fig:dist_effect_relative_motion_cas]].
#+begin_src matlab :exports none
figure;
hold on;
plot(gm.f, sqrt(gm.psd_rv), 'DisplayName', 'Ground Motion');
plot(tyz.f, sqrt(tyz.pxz_ty_r), 'DisplayName', 'Ty');
plot(rz.f, sqrt(rz.pxsp_r), 'DisplayName', 'Rz');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD of the measured velocity $\left[\frac{m/s}{\sqrt{Hz}}\right]$')
legend('Location', 'southwest');
xlim([2, 500]);
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/dist_effect_relative_velocity.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:dist_effect_relative_velocity
#+CAPTION: Amplitude Spectral Density of the relative velocity of the hexapod with respect to the granite due to different sources of perturbation ([[./figs/dist_effect_relative_velocity.png][png]], [[./figs/dist_effect_relative_velocity.pdf][pdf]])
[[file:figs/dist_effect_relative_velocity.png]]
#+begin_src matlab :exports none
figure;
hold on;
plot(gm.f, sqrt(gm.psd_rv)./(2*pi*gm.f), 'DisplayName', 'Ground Motion');
plot(tyz.f, sqrt(tyz.pxz_ty_r)./(2*pi*tyz.f), 'DisplayName', 'Ty');
plot(rz.f, sqrt(rz.pxsp_r)./(2*pi*rz.f), 'DisplayName', 'Rz');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD of the displacement $\left[\frac{m}{\sqrt{Hz}}\right]$')
legend('Location', 'southwest');
xlim([2, 500]);
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/dist_effect_relative_motion.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:dist_effect_relative_motion
#+CAPTION: Amplitude Spectral Density of the relative displacement of the hexapod with respect to the granite due to different sources of perturbation ([[./figs/dist_effect_relative_motion.png][png]], [[./figs/dist_effect_relative_motion.pdf][pdf]])
[[file:figs/dist_effect_relative_motion.png]]
#+begin_src matlab :exports none
figure;
hold on;
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(gm.psd_rv./((2*pi*gm.f).^2))))), 'DisplayName', 'Ground Motion');
plot(tyz.f, flip(sqrt(-cumtrapz(flip(tyz.f), flip(tyz.pxz_ty_r./((2*pi*tyz.f).^2))))), 'DisplayName', 'Ty');
plot(rz.f, flip(sqrt(-cumtrapz(flip(rz.f), flip(rz.pxsp_r./((2*pi*rz.f).^2))))), 'DisplayName', 'Rz');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('CAS of the relative displacement $[m]$')
legend('Location', 'southwest');
xlim([2, 500]); ylim([1e-11, 1e-6]);
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/dist_effect_relative_motion_cas.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:dist_effect_relative_motion_cas
#+CAPTION: Cumulative Amplitude Spectrum of the relative motion due to different sources of perturbation ([[./figs/dist_effect_relative_motion_cas.png][png]], [[./figs/dist_effect_relative_motion_cas.pdf][pdf]])
[[file:figs/dist_effect_relative_motion_cas.png]]
* Compute the Power Spectral Density of the disturbance force
<<sec:psd_force_dist>>
Now, from the extracted transfer functions from the disturbance force to the relative motion of the hexapod with respect to the granite (section [[sec:sensitivity_disturbances]]) and from the measured PSD of the relative motion (section [[sec:psd_dist]]), we can compute the PSD of the disturbance force.
#+begin_src matlab
rz.psd_f = rz.pxsp_r./abs(squeeze(freqresp(G('Vm', 'Frz'), rz.f, 'Hz'))).^2;
tyz.psd_f = tyz.pxz_ty_r./abs(squeeze(freqresp(G('Vm', 'Fty'), tyz.f, 'Hz'))).^2;
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
set(gca,'ColorOrderIndex',2);
plot(tyz.f, sqrt(tyz.psd_f), 'DisplayName', 'F - Ty');
plot(rz.f, sqrt(rz.psd_f), 'DisplayName', 'F - Rz');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD of the disturbance force $\left[\frac{F}{\sqrt{Hz}}\right]$')
legend('Location', 'southwest');
xlim([2, 500]);
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/dist_force_psd.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:dist_force_psd
#+CAPTION: Amplitude Spectral Density of the disturbance force ([[./figs/dist_force_psd.png][png]], [[./figs/dist_force_psd.pdf][pdf]])
[[file:figs/dist_force_psd.png]]
* Noise Budget
<<sec:noise_budget>>
Now, from the compute spectral density of the disturbance sources, we can compute the resulting relative motion of the Hexapod with respect to the granite using the model.
We should verify that this is coherent with the measurements.
#+begin_src matlab :exports none
% Power Spectral Density of the relative Displacement
psd_gm_d = gm.psd_gm.*abs(squeeze(freqresp(G('Vm', 'Dw')/s, gm.f, 'Hz'))).^2;
psd_ty_d = tyz.psd_f.*abs(squeeze(freqresp(G('Vm', 'Fty')/s, tyz.f, 'Hz'))).^2;
psd_rz_d = rz.psd_f.*abs(squeeze(freqresp(G('Vm', 'Frz')/s, rz.f, 'Hz'))).^2;
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
plot(gm.f, sqrt(psd_gm_d), 'DisplayName', 'Ground Motion');
plot(tyz.f, sqrt(psd_ty_d), 'DisplayName', 'Ty');
plot(rz.f, sqrt(psd_rz_d), 'DisplayName', 'Rz');
plot(rz.f, sqrt(psd_gm_d + psd_ty_d + psd_rz_d), 'k--', 'DisplayName', 'tot');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD of the relative motion $\left[\frac{m}{\sqrt{Hz}}\right]$')
legend('Location', 'southwest');
xlim([2, 500]);
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/psd_effect_dist_verif.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:psd_effect_dist_verif
#+CAPTION: Computed Effect of the disturbances on the relative displacement hexapod/granite ([[./figs/psd_effect_dist_verif.png][png]], [[./figs/psd_effect_dist_verif.pdf][pdf]])
[[file:figs/psd_effect_dist_verif.png]]
#+begin_src matlab :exports none
figure;
hold on;
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(psd_gm_d)))), 'DisplayName', 'Ground Motion');
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(psd_ty_d)))), 'DisplayName', 'Ty');
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(psd_rz_d)))), 'DisplayName', 'Rz');
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(psd_gm_d + psd_ty_d + psd_rz_d)))), 'k-', 'DisplayName', 'tot');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('Cumulative Amplitude Spectrum [m]')
legend('location', 'northeast');
xlim([2, 500]); ylim([1e-11, 1e-6]);
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/cas_computed_relative_displacement.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:cas_computed_relative_displacement
#+CAPTION: CAS of the total Relative Displacement due to all considered sources of perturbation ([[./figs/cas_computed_relative_displacement.png][png]], [[./figs/cas_computed_relative_displacement.pdf][pdf]])
[[file:figs/cas_computed_relative_displacement.png]]
* Approximation
We approximate the PSD of the disturbance with the following transfer functions.
#+begin_src matlab
G_ty = 0.1*(s+634.3)*(s+283.7)/((s+2*pi)*(s+2*pi));
G_rz = 0.5*(s+418.8)*(s+36.51)*(s^2 + 110.9*s + 3.375e04)/((s+0.7324)*(s+0.546)*(s^2 + 0.6462*s + 2.391e04));
G_gm = 0.002*(s^2 + 3.169*s + 27.74)/(s*(s+32.73)*(s+8.829)*(s+7.983)^2);
#+end_src
We compute the effect of these approximate disturbances on $D$.
#+begin_src matlab :exports none
% Power Spectral Density of the relative Displacement
psd_gm_s = abs(squeeze(freqresp(G_gm*G('Vm', 'Dw')/s, gm.f, 'Hz'))).^2;
psd_ty_s = abs(squeeze(freqresp(G_ty*G('Vm', 'Fty')/s, gm.f, 'Hz'))).^2;
psd_rz_s = abs(squeeze(freqresp(G_rz*G('Vm', 'Frz')/s, gm.f, 'Hz'))).^2;
#+end_src
#+begin_src matlab :exports none
figure;
ax1 = subplot(1, 2, 1);
hold on;
set(gca,'ColorOrderIndex',2);
plot(gm.f, sqrt(psd_ty_d), 'DisplayName', 'F - Ty');
plot(gm.f, sqrt(psd_rz_d), 'DisplayName', 'F - Rz');
set(gca,'ColorOrderIndex',2);
plot(gm.f, sqrt(psd_ty_s), '--', 'HandleVisibility', 'off');
plot(gm.f, sqrt(psd_rz_s), '--', 'HandleVisibility', 'off');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD of the disturbance force $\left[\frac{F}{\sqrt{Hz}}\right]$')
legend('Location', 'southwest');
xlim([2, 500]);
ax2 = subplot(1, 2, 2);
hold on;
plot(gm.f, sqrt(psd_gm_d), 'DisplayName', 'D - Gm');
set(gca,'ColorOrderIndex',1);
plot(gm.f, sqrt(psd_gm_s), '--', 'HandleVisibility', 'off');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD of the ground displacement $\left[\frac{m}{\sqrt{Hz}}\right]$')
legend('Location', 'southwest');
xlim([2, 500]);
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/estimate_spectral_density_disturbances.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:estimate_spectral_density_disturbances
#+CAPTION: Estimated spectral density of the disturbances ([[./figs/estimate_spectral_density_disturbances.png][png]], [[./figs/estimate_spectral_density_disturbances.pdf][pdf]])
[[file:figs/estimate_spectral_density_disturbances.png]]
#+begin_src matlab :exports none
figure;
hold on;
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(psd_gm_d)))), 'DisplayName', 'Gm');
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(psd_ty_d)))), 'DisplayName', 'Ty');
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(psd_rz_d)))), 'DisplayName', 'Rz');
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(psd_gm_d + psd_ty_d + psd_rz_d)))), 'k-', 'DisplayName', 'tot');
set(gca,'ColorOrderIndex',1);
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(psd_gm_s)))), '--', 'HandleVisibility', 'off');
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(psd_ty_s)))), '--', 'HandleVisibility', 'off');
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(psd_rz_s)))), '--', 'HandleVisibility', 'off');
plot(gm.f, flip(sqrt(-cumtrapz(flip(gm.f), flip(psd_gm_s + psd_ty_s + psd_rz_s)))), 'k--', 'HandleVisibility', 'off');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('Cumulative Amplitude Spectrum [m]')
legend('location', 'northeast');
xlim([2, 500]); ylim([1e-11, 1e-6]);
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/comp_estimation_cas_disturbances.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:comp_estimation_cas_disturbances
#+CAPTION: Comparison of the CAS of the disturbances with the approximate ones ([[./figs/comp_estimation_cas_disturbances.png][png]], [[./figs/comp_estimation_cas_disturbances.pdf][pdf]])
[[file:figs/comp_estimation_cas_disturbances.png]]
* Save
The PSD of the disturbance force are now saved for further noise budgeting when control is applied (the mat file is accessible [[file:mat/dist_psd.mat][here]]).
#+begin_src matlab
dist_f = struct();
dist_f.f = gm.f; % Frequency Vector [Hz]
dist_f.psd_gm = gm.psd_gm; % Power Spectral Density of the Ground Motion [m^2/Hz]
dist_f.psd_ty = tyz.psd_f; % Power Spectral Density of the force induced by the Ty stage in the Z direction [N^2/Hz]
dist_f.psd_rz = rz.psd_f; % Power Spectral Density of the force induced by the Rz stage in the Z direction [N^2/Hz]
dist_f.G_gm = G_ty;
dist_f.G_ty = G_rz;
dist_f.G_rz = G_gm;
save('./disturbances/mat/dist_psd.mat', 'dist_f');
#+end_src
* Function that initialize the disturbances
:PROPERTIES:
:header-args:matlab+: :tangle ../src/initDisturbances.m
:header-args:matlab+: :comments none :mkdirp yes
:header-args:matlab+: :eval no :results none
:END:
<<sec:initDisturbances>>
This Matlab function is accessible [[file:src/initDisturbances.m][here]].
#+begin_src matlab
function [] = initDisturbances(opts_param)
% initDisturbances - Initialize the disturbances
%
% Syntax: [] = initDisturbances(opts_param)
%
% Inputs:
% - opts_param -
#+end_src
** Default values for the Options
#+begin_src matlab
%% Default values for opts
opts = struct();
%% Populate opts with input parameters
if exist('opts_param','var')
for opt = fieldnames(opts_param)'
opts.(opt{1}) = opts_param.(opt{1});
end
end
#+end_src
** Load Data
#+begin_src matlab
load('./disturbances/mat/dist_psd.mat', 'dist_f');
#+end_src
We remove the first frequency point that usually is very large.
#+begin_src matlab :exports none
dist_f.f = dist_f.f(2:end);
dist_f.psd_gm = dist_f.psd_gm(2:end);
dist_f.psd_ty = dist_f.psd_ty(2:end);
dist_f.psd_rz = dist_f.psd_rz(2:end);
#+end_src
** Parameters
We define some parameters that will be used in the algorithm.
#+begin_src matlab
Fs = 2*dist_f.f(end); % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
N = 2*length(dist_f.f); % Number of Samples match the one of the wanted PSD
T0 = N/Fs; % Signal Duration [s]
df = 1/T0; % Frequency resolution of the DFT [Hz]
% Also equal to (dist_f.f(2)-dist_f.f(1))
t = linspace(0, T0, N+1)'; % Time Vector [s]
Ts = 1/Fs; % Sampling Time [s]
#+end_src
** Ground Motion
#+begin_src matlab
phi = dist_f.psd_gm;
C = zeros(N/2,1);
for i = 1:N/2
C(i) = sqrt(phi(i)*df);
end
theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
Cx = [0 ; C.*complex(cos(theta),sin(theta))];
Cx = [Cx; flipud(conj(Cx(2:end)))];;
u = N/sqrt(2)*ifft(Cx); % Ground Motion - x direction [m]
% Dwx = struct('time', t, 'signals', struct('values', u));
Dwx = u;
theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
Cx = [0 ; C.*complex(cos(theta),sin(theta))];
Cx = [Cx; flipud(conj(Cx(2:end)))];;
u = N/sqrt(2)*ifft(Cx); % Ground Motion - y direction [m]
Dwy = u;
theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
Cx = [0 ; C.*complex(cos(theta),sin(theta))];
Cx = [Cx; flipud(conj(Cx(2:end)))];;
u = N/sqrt(2)*ifft(Cx); % Ground Motion - z direction [m]
Dwz = u;
#+end_src
** Translation Stage - X direction
#+begin_src matlab
phi = dist_f.psd_ty; % TODO - we take here the vertical direction which is wrong but approximate
C = zeros(N/2,1);
for i = 1:N/2
C(i) = sqrt(phi(i)*df);
end
theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
Cx = [0 ; C.*complex(cos(theta),sin(theta))];
Cx = [Cx; flipud(conj(Cx(2:end)))];;
u = N/sqrt(2)*ifft(Cx); % Disturbance Force Ty x [N]
Fty_x = u;
#+end_src
** Translation Stage - Z direction
#+begin_src matlab
phi = dist_f.psd_ty;
C = zeros(N/2,1);
for i = 1:N/2
C(i) = sqrt(phi(i)*df);
end
theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
Cx = [0 ; C.*complex(cos(theta),sin(theta))];
Cx = [Cx; flipud(conj(Cx(2:end)))];;
u = N/sqrt(2)*ifft(Cx); % Disturbance Force Ty z [N]
Fty_z = u;
#+end_src
** Spindle - Z direction
#+begin_src matlab
phi = dist_f.psd_rz;
C = zeros(N/2,1);
for i = 1:N/2
C(i) = sqrt(phi(i)*df);
end
theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
Cx = [0 ; C.*complex(cos(theta),sin(theta))];
Cx = [Cx; flipud(conj(Cx(2:end)))];;
u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N]
Frz_z = u;
#+end_src
** Direct Forces
#+begin_src matlab
u = zeros(length(t), 6);
Fd = u;
#+end_src
** Set initial value to zero
#+begin_src matlab
Dwx = Dwx - Dwx(1);
Dwy = Dwy - Dwy(1);
Dwz = Dwz - Dwz(1);
Fty_x = Fty_x - Fty_x(1);
Fty_z = Fty_z - Fty_z(1);
Frz_z = Frz_z - Frz_z(1);
#+end_src
** Save
#+begin_src matlab
save('./mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_z', 'Fd', 'Ts', 't');
#+end_src
* Display Obtained Disturbances
#+begin_src matlab
initDisturbances();
load('./mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_z', 'Fd', 'Ts', 't');
#+end_src
#+begin_src matlab
figure;
hold on;
plot(t, Dwx, 'DisplayName', 'Dw x');
plot(t, Dwy, 'DisplayName', 'Dw y');
plot(t, Dwz, 'DisplayName', 'Dw z');
hold off;
xlabel('Time [s]'); ylabel('Amplitude [m]');
legend('location', 'north east');
#+end_src
#+begin_src matlab
figure;
hold on;
plot(t, Fty_x, 'DisplayName', 'Ty x');
plot(t, Fty_z, 'DisplayName', 'Ty z');
plot(t, Frz_z, 'DisplayName', 'Rz z');
hold off;
xlabel('Time [s]'); ylabel('Force [N]');
legend('location', 'north east');
#+end_src
Reference in New Issue View Git Blame Copy Permalink