tdehaeze/nass-metrology-test-bench
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
61067522d0
nass-metrology-test-bench/index.org

572 lines
17 KiB
Org Mode
Raw Blame History

#+TITLE: Cercalo Test Bench
:DRAWER:
#+STARTUP: overview
#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="./css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="./css/readtheorg.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="./css/zenburn.css"/>
#+HTML_HEAD: <script type="text/javascript" src="./js/jquery.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="./js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="./js/jquery.stickytableheaders.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="./js/readtheorg.js"></script>
#+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/frf_processing.m
#+PROPERTY: header-args:matlab+ :mkdirp yes
:END:
* Introduction
* Identification
:PROPERTIES:
:header-args:matlab+: :tangle matlab/plant_identification.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
<<sec:plant_identification>>
** ZIP file containing the data and matlab files :ignore:
#+begin_src bash :exports none :results none
if [ matlab/plant_identification.m -nt data/plant_identification.zip ]; then
cp matlab/plant_identification.m plant_identification.m;
zip data/plant_identification \
mat/data_ux.mat \
mat/data_uy.mat \
plant_identification.m
rm plant_identification.m;
fi
#+end_src
#+begin_note
All the files (data and Matlab scripts) are accessible [[file:data/plant_identification.zip][here]].
#+end_note
** 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
** Excitation Data
#+begin_src matlab
fs = 1e4;
Ts = 1/fs;
#+end_src
We generate white noise with the "random number" simulink block, and we filter that noise.
#+begin_src matlab
Gi = (1)/(1+s/2/pi/100);
#+end_src
#+begin_src matlab :results output replace
c2d(Gi, Ts, 'tustin')
#+end_src
#+RESULTS:
#+begin_example
c2d(Gi, Ts, 'tustin')
ans =
0.030459 (z+1)
--------------
(z-0.9391)
Sample time: 0.0001 seconds
Discrete-time zero/pole/gain model.
#+end_example
** Huddle Test
We load the data taken during the Huddle Test.
#+begin_src matlab
load('mat/data_huddle_test.mat', 't', 'xh', 'xv', 'cuh', 'cuv');
#+end_src
The variables are:
- $x_h$ :: Normalized position of the beam in the horizontal direction as measured by the 4 quadrant diode
- $x_v$ :: Normalized position of the beam in the vertical direction as measured by the 4 quadrant diode
- $cu_h$ :: Voltage across the inductance used for the horizontal positioning of the Cercalo
- $vu_v$ :: Voltage across the inductance used for the vertical positioning of the Cercalo
We remove the first second of data where everything is settling down.
#+begin_src matlab
xh(t<1) = [];
xv(t<1) = [];
cuh(t<1) = [];
cuv(t<1) = [];
t(t<1) = [];
t = t - t(1);
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
plot(t, xh, 'DisplayName', '$x_h$');
plot(t, xv, 'DisplayName', '$x_v$');
hold off;
xlabel('Time [s]');
ylabel('Amplitude [V]');
xlim([t(1), t(end)]);
legend();
#+end_src
We compute the Power Spectral Density of the horizontal and vertical positions of the beam as measured by the 4 quadrant diode.
#+begin_src matlab
[psd_xh, f] = pwelch(xh, hanning(ceil(1*fs)), [], [], fs);
[psd_xv, ~] = pwelch(xv, hanning(ceil(1*fs)), [], [], fs);
#+end_src
#+begin_src matlab :results none
figure;
hold on;
plot(f, sqrt(psd_xh), 'DisplayName', '$\Gamma_{x_h}$');
plot(f, sqrt(psd_xv), 'DisplayName', '$\Gamma_{x_v}$');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD $\left[\frac{V}{\sqrt{Hz}}\right]$')
legend('Location', 'southwest');
xlim([1, 1000]);
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
plot(t, cuh, 'DisplayName', '$cu_h$');
plot(t, cuv, 'DisplayName', '$cu_v$');
hold off;
xlabel('Time [s]');
ylabel('Amplitude [V]');
xlim([t(1), t(end)]);
legend();
#+end_src
We compute the Power Spectral Density of the voltage across the inductance used for horizontal and vertical positioning of the Cercalo.
#+begin_src matlab
[psd_cuh, f] = pwelch(cuh, hanning(ceil(1*fs)), [], [], fs);
[psd_cuv, ~] = pwelch(cuv, hanning(ceil(1*fs)), [], [], fs);
#+end_src
#+begin_src matlab :results none
figure;
hold on;
plot(f, sqrt(psd_cuh), 'DisplayName', '$\Gamma_{cu_h}$');
plot(f, sqrt(psd_cuv), 'DisplayName', '$\Gamma_{cu_v}$');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD $\left[\frac{V}{\sqrt{Hz}}\right]$')
legend('Location', 'southwest');
xlim([1, 1000]);
#+end_src
** Input / Output data
The identification data is loaded
#+begin_src matlab
uh = load('mat/data_uh.mat', 't', 'xh', 'xv', 'uh', 'cuh', 'cuv');
uv = load('mat/data_uv.mat', 't', 'xh', 'xv', 'uv', 'cuh', 'cuv');
#+end_src
We remove the first seconds where the Cercalo is turned on.
#+begin_src matlab
i0x = 20*fs;
i0y = 10*fs;
uh.t = uh.t( i0x:end) - uh.t(i0x);
uh.uh = uh.uh(i0x:end);
uh.xh = uh.xh(i0x:end);
uh.xv = uh.xv(i0x:end);
uv.t = uv.t( i0y:end) - uv.t(i0x);
uv.uv = uv.uv(i0y:end);
uv.xh = uv.xh(i0y:end);
uv.xv = uv.xv(i0y:end);
#+end_src
#+begin_src matlab
uh.uh = uh.uh-mean(uh.uh);
uh.xh = uh.xh-mean(uh.xh);
uh.xv = uh.xv-mean(uh.xv);
uv.uv = uv.uv-mean(uv.uv);
uv.xh = uv.xh-mean(uv.xh);
uv.xv = uv.xv-mean(uv.xv);
#+end_src
#+begin_src matlab :exports none
figure;
ax1 = subplot(1, 2, 1);
plot(uh.t, uh.uh);
xlabel('Time [s]');
ylabel('Amplitude [V]');
legend({'$u_x$'});
ax2 = subplot(1, 2, 2);
hold on;
plot(uh.t, uh.xh, 'DisplayName', '$y_x$');
plot(uh.t, uh.xv, 'DisplayName', '$y_y$');
hold off;
xlabel('Time [s]');
ylabel('Amplitude [V]');
legend()
linkaxes([ax1,ax2],'x');
xlim([uh.t(1), uh.t(end)])
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/identification_uh.pdf" :var figsize="full-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:identification_uh
#+CAPTION: Identification signals when exciting the $x$ axis ([[./figs/identification_uh.png][png]], [[./figs/identification_uh.pdf][pdf]])
[[file:figs/identification_uh.png]]
#+begin_src matlab :exports none
figure;
ax1 = subplot(1, 2, 1);
plot(uv.t, uv.uv);
xlabel('Time [s]');
ylabel('Amplitude [V]');
legend({'$u_y$'});
ax2 = subplot(1, 2, 2);
hold on;
plot(uv.t, uv.xh, 'DisplayName', '$y_x$');
plot(uv.t, uv.xv, 'DisplayName', '$y_y$');
hold off;
xlabel('Time [s]');
ylabel('Amplitude [V]');
legend()
linkaxes([ax1,ax2],'x');
xlim([uv.t(1), uv.t(end)])
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/identification_uv.pdf" :var figsize="full-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:identification_uv
#+CAPTION: Identification signals when exciting the $y$ axis ([[./figs/identification_uv.png][png]], [[./figs/identification_uv.pdf][pdf]])
[[file:figs/identification_uv.png]]
** Estimation of the Frequency Response Function Matrix
We compute an estimate of the transfer functions.
#+begin_src matlab
[tf_uh_xh, f] = tfestimate(uh.uh, uh.xh, hanning(ceil(1*fs)), [], [], fs);
[tf_uh_xv, ~] = tfestimate(uh.uh, uh.xv, hanning(ceil(1*fs)), [], [], fs);
[tf_uv_xh, ~] = tfestimate(uv.uv, uv.xh, hanning(ceil(1*fs)), [], [], fs);
[tf_uv_xv, ~] = tfestimate(uv.uv, uv.xv, hanning(ceil(1*fs)), [], [], fs);
#+end_src
#+begin_src matlab :exports none
figure;
ax11 = subplot(2, 2, 1);
hold on;
plot(f, abs(tf_uh_xh))
title('Frequency Response Function $\frac{y_x}{u_x}$')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude')
hold off;
ax12 = subplot(2, 2, 2);
hold on;
plot(f, abs(tf_uv_xh))
title('Frequency Response Function $\frac{y_x}{u_y}$')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
hold off;
ax21 = subplot(2, 2, 3);
hold on;
plot(f, abs(tf_uh_xv))
title('Frequency Response Function $\frac{y_y}{u_x}$')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude')
xlabel('Frequency [Hz]')
hold off;
ax22 = subplot(2, 2, 4);
hold on;
plot(f, abs(tf_uv_xv))
title('Frequency Response Function $\frac{y_y}{u_y}$')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
xlabel('Frequency [Hz]')
hold off;
linkaxes([ax11,ax12,ax21,ax22],'x');
xlim([10, 1000]);
linkaxes([ax11,ax12,ax21,ax22],'y');
ylim([1e-2, 1e3])
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/frequency_response_matrix.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:frequency_response_matrix
#+CAPTION: Frequency Response Matrix ([[./figs/frequency_response_matrix.png][png]], [[./figs/frequency_response_matrix.pdf][pdf]])
[[file:figs/frequency_response_matrix.png]]
** Coherence
#+begin_src matlab
[coh_uh_xh, f] = mscohere(uh.uh, uh.xh, hanning(ceil(1*fs)), [], [], fs);
[coh_uh_xv, ~] = mscohere(uh.uh, uh.xv, hanning(ceil(1*fs)), [], [], fs);
[coh_uv_xh, ~] = mscohere(uv.uv, uv.xh, hanning(ceil(1*fs)), [], [], fs);
[coh_uv_xv, ~] = mscohere(uv.uv, uv.xv, hanning(ceil(1*fs)), [], [], fs);
#+end_src
#+begin_src matlab :exports none
figure;
ax11 = subplot(2, 2, 1);
hold on;
plot(f, coh_uh_xh)
set(gca, 'Xscale', 'log');
title('Coherence $\frac{y_x}{u_x}$')
ylabel('Coherence')
hold off;
ax12 = subplot(2, 2, 2);
hold on;
plot(f, coh_uv_xh)
set(gca, 'Xscale', 'log');
title('Coherence $\frac{y_x}{u_y}$')
hold off;
ax21 = subplot(2, 2, 3);
hold on;
plot(f, coh_uh_xv)
set(gca, 'Xscale', 'log');
title('Coherence $\frac{y_y}{u_x}$')
ylabel('Coherence')
xlabel('Frequency [Hz]')
hold off;
ax22 = subplot(2, 2, 4);
hold on;
plot(f, coh_uv_xv)
set(gca, 'Xscale', 'log');
title('Coherence $\frac{y_y}{u_y}$')
xlabel('Frequency [Hz]')
hold off;
linkaxes([ax11,ax12,ax21,ax22],'x');
xlim([10, 1000]);
linkaxes([ax11,ax12,ax21,ax22],'y');
ylim([0, 1])
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/identification_coherence.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:identification_coherence
#+CAPTION: Coherence ([[./figs/identification_coherence.png][png]], [[./figs/identification_coherence.pdf][pdf]])
[[file:figs/identification_coherence.png]]
** Extraction of a transfer function matrix
First we define the initial guess for the resonance frequencies and the weights associated.
#+begin_src matlab
freqs_res = [410, 250]; % [Hz]
freqs_res_weights = [10, 10]; % [Hz]
#+end_src
From the number of resonance frequency we want to fit, we define the order =N= of the system we want to obtain.
#+begin_src matlab
N = 2*length(freqs_res);
#+end_src
We then make an initial guess on the complex values of the poles.
#+begin_src matlab
xi = 0.001; % Approximate modal damping
poles = [2*pi*freqs_res*(xi + 1i), 2*pi*freqs_res*(xi - 1i)];
#+end_src
We then define the weight that will be used for the fitting.
Basically, we want more weight around the resonance and at low frequency (below the first resonance).
Also, we want more importance where we have a better coherence.
#+begin_src matlab
weight = ones(1, length(f));
% weight = G_coh';
% alpha = 0.1;
% for freq_i = 1:length(freqs_res)
% weight(f>(1-alpha)*freqs_res(freq_i) & omega<(1 + alpha)*2*pi*freqs_res(freq_i)) = freqs_res_weights(freq_i);
% end
#+end_src
Ignore data above some frequency.
#+begin_src matlab
weight(f>1000) = 0;
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
plot(f, weight);
plot(freqs_res, ones(size(freqs_res)), 'rx');
hold off;
xlabel('Frequency [Hz]');
xlabel('Weight Amplitude');
set(gca, 'xscale', 'log');
xlim([f(1), f(end)]);
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/weights.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:weights
#+CAPTION: Weights amplitude ([[./figs/weights.png][png]], [[./figs/weights.pdf][pdf]])
[[file:figs/weights.png]]
When we set some options for =vfit3=.
#+begin_src matlab
opts = struct();
opts.stable = 1; % Enforce stable poles
opts.asymp = 1; % Force D matrix to be null
opts.relax = 1; % Use vector fitting with relaxed non-triviality constraint
opts.skip_pole = 0; % Do NOT skip pole identification
opts.skip_res = 0; % Do NOT skip identification of residues (C,D,E)
opts.cmplx_ss = 0; % Create real state space model with block diagonal A
opts.spy1 = 0; % No plotting for first stage of vector fitting
opts.spy2 = 0; % Create magnitude plot for fitting of f(s)
#+end_src
We define the number of iteration.
#+begin_src matlab
Niter = 5;
#+end_src
An we run the =vectfit3= algorithm.
#+begin_src matlab
for iter = 1:Niter
[SER_uh_xh, poles, ~, fit_uh_xh] = vectfit3(tf_uh_xh.', 1i*2*pi*f, poles, weight, opts);
end
for iter = 1:Niter
[SER_uv_xh, poles, ~, fit_uv_xh] = vectfit3(tf_uv_xh.', 1i*2*pi*f, poles, weight, opts);
end
for iter = 1:Niter
[SER_uh_xv, poles, ~, fit_uh_xv] = vectfit3(tf_uh_xv.', 1i*2*pi*f, poles, weight, opts);
end
for iter = 1:Niter
[SER_uv_xv, poles, ~, fit_uv_xv] = vectfit3(tf_uv_xv.', 1i*2*pi*f, poles, weight, opts);
end
#+end_src
#+begin_src matlab :exports none
figure;
ax11 = subplot(2, 2, 1);
hold on;
plot(f, abs(tf_uh_xh))
plot(f, abs(fit_uh_xh))
title('Frequency Response Function $\frac{y_x}{u_x}$')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude')
hold off;
ax12 = subplot(2, 2, 2);
hold on;
plot(f, abs(tf_uv_xh))
plot(f, abs(fit_uv_xh))
title('Frequency Response Function $\frac{y_x}{u_y}$')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
hold off;
ax21 = subplot(2, 2, 3);
hold on;
plot(f, abs(tf_uh_xv))
plot(f, abs(fit_uh_xv))
title('Frequency Response Function $\frac{y_y}{u_x}$')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude')
xlabel('Frequency [Hz]')
hold off;
ax22 = subplot(2, 2, 4);
hold on;
plot(f, abs(tf_uv_xv))
plot(f, abs(fit_uv_xv))
title('Frequency Response Function $\frac{y_y}{u_y}$')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
xlabel('Frequency [Hz]')
hold off;
linkaxes([ax11,ax12,ax21,ax22],'x');
xlim([10, 1000]);
linkaxes([ax11,ax12,ax21,ax22],'y');
ylim([1e-2, 1e3])
#+end_src
#+HEADER: :tangle no :exports results :results none :noweb yes
#+begin_src matlab :var filepath="figs/identification_matrix_fit.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:identification_matrix_fit
#+CAPTION: Transfer Function Extraction of the FRF matrix ([[./figs/identification_matrix_fit.png][png]], [[./figs/identification_matrix_fit.pdf][pdf]])
[[file:figs/identification_matrix_fit.png]]
And finally, we create the identified state space model:
#+begin_src matlab
G_uh_xh = minreal(ss(full(SER_uh_xh.A),SER_uh_xh.B,SER_uh_xh.C,SER_uh_xh.D));
G_uv_xh = minreal(ss(full(SER_uv_xh.A),SER_uv_xh.B,SER_uv_xh.C,SER_uv_xh.D));
G_uh_xv = minreal(ss(full(SER_uh_xv.A),SER_uh_xv.B,SER_uh_xv.C,SER_uh_xv.D));
G_uv_xv = minreal(ss(full(SER_uv_xv.A),SER_uv_xv.B,SER_uv_xv.C,SER_uv_xv.D));
G = [G_uh_xh, G_uv_xh;
G_uh_xv, G_uv_xv];
#+end_src
#+begin_src matlab
save('mat/plant.mat', 'G');
#+end_src
* Sensor Noise
* Plant Analysis
** Rotation Matrix
#+begin_src matlab
G0 = freqresp(G, 0);
#+end_src
* Control Objective
The maximum expected stroke is $y_\text{max} = 3mm \approx 5e^{-2} rad$ at $1Hz$.
The maximum wanted error is $e_\text{max} = 10 \mu rad$.
Thus, we require the sensitivity function at $\omega_0 = 1\text{ Hz}$:
\begin{align*}
|S(j\omega_0)| &< \left| \frac{e_\text{max}}{y_\text{max}} \right| \\
&< 2 \cdot 10^{-4}
\end{align*}
In terms of loop gain, this is equivalent to:
\[ |L(j\omega_0)| > 5 \cdot 10^{3} \]
* Plant Scaling
- measured noise
- expected perturbations
- maximum input usage
- maximum wanted error
* Control Design
Reference in New Issue View Git Blame Copy Permalink