nass-metrology-test-bench/index.org

Cercalo Test Bench

• Introduction
• Identification
  • ZIP file containing the data and matlab files
  • Excitation Data
  • Huddle Test
  • Input / Output data
  • Estimation of the Frequency Response Function Matrix
  • Coherence
  • Extraction of a transfer function matrix
• Sensor Noise
• Plant Analysis
  • Rotation Matrix
• Control Objective
• Plant Scaling
• Control Design

Introduction

Identification

<<sec:plant_identification>>

ZIP file containing the data and matlab files   ignore

All the files (data and Matlab scripts) are accessible here.

Excitation Data

  fs = 1e4;
  Ts = 1/fs;

We generate white noise with the "random number" simulink block, and we filter that noise.

  Gi = (1)/(1+s/2/pi/100);

  c2d(Gi, Ts, 'tustin')

c2d(Gi, Ts, 'tustin')

ans =

  0.030459 (z+1)
  --------------
    (z-0.9391)

Sample time: 0.0001 seconds
Discrete-time zero/pole/gain model.

Huddle Test

We load the data taken during the Huddle Test.

  load('mat/data_huddle_test.mat', 't', 'xh', 'xv', 'cuh', 'cuv');

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.

  xh(t<1) = [];
  xv(t<1) = [];
  cuh(t<1) = [];
  cuv(t<1) = [];
  t(t<1) = [];
  t = t - t(1);

We compute the Power Spectral Density of the horizontal and vertical positions of the beam as measured by the 4 quadrant diode.

  [psd_xh, f] = pwelch(xh, hanning(ceil(1*fs)), [], [], fs);
  [psd_xv, ~] = pwelch(xv, hanning(ceil(1*fs)), [], [], fs);

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

We compute the Power Spectral Density of the voltage across the inductance used for horizontal and vertical positioning of the Cercalo.

  [psd_cuh, f] = pwelch(cuh, hanning(ceil(1*fs)), [], [], fs);
  [psd_cuv, ~] = pwelch(cuv, hanning(ceil(1*fs)), [], [], fs);

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

Input / Output data

The identification data is loaded

  uh = load('mat/data_uh.mat', 't', 'xh', 'xv', 'uh', 'cuh', 'cuv');
  uv = load('mat/data_uv.mat', 't', 'xh', 'xv', 'uv', 'cuh', 'cuv');

We remove the first seconds where the Cercalo is turned on.

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

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

  <<plt-matlab>>

Identification signals when exciting the $x$ axis (png, pdf)

  <<plt-matlab>>

Identification signals when exciting the $y$ axis (png, pdf)

Estimation of the Frequency Response Function Matrix

We compute an estimate of the transfer functions.

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

  <<plt-matlab>>

Frequency Response Matrix (png, pdf)

Coherence

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

  <<plt-matlab>>

Coherence (png, pdf)

Extraction of a transfer function matrix

First we define the initial guess for the resonance frequencies and the weights associated.

  freqs_res = [410, 250]; % [Hz]
  freqs_res_weights = [10, 10]; % [Hz]

From the number of resonance frequency we want to fit, we define the order N of the system we want to obtain.

  N = 2*length(freqs_res);

We then make an initial guess on the complex values of the poles.

  xi = 0.001; % Approximate modal damping
  poles = [2*pi*freqs_res*(xi + 1i), 2*pi*freqs_res*(xi - 1i)];

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.

  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

Ignore data above some frequency.

  weight(f>1000) = 0;

  <<plt-matlab>>

Weights amplitude (png, pdf)

When we set some options for vfit3.

  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)

We define the number of iteration.

  Niter = 5;

An we run the vectfit3 algorithm.

  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

  <<plt-matlab>>

Transfer Function Extraction of the FRF matrix (png, pdf)

And finally, we create the identified state space model:

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

  save('mat/plant.mat', 'G');

Sensor Noise

Plant Analysis

Rotation Matrix

  G0 = freqresp(G, 0);

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