nass-metrology-test-bench/index.org

Cercalo Test Bench

• Introduction
  • Block Diagram
  • Amplifier for the Cercalo
  • Cercalo
  • Optical Setup
  • 4 quadrant Diode
  • ADC/DAC
• Identification
  • Introduction
  • ZIP file containing the data and matlab files
  • Excitation Data
  • Signals
  • Huddle Test
  • Input / Output data
  • Estimation of the Frequency Response Function Matrix
  • Coherence
  • Extraction of a transfer function matrix
• Plant Scaling
• Plant Analysis
  • Load Plant
  • RGA-Number
  • Rotation Matrix
• Control Objective
• Control Design
• Measurement of the non-repeatability

Introduction

Block Diagram

The block diagram of the setup to be controlled is shown in Fig. fig:block_diagram_simplify.

Block Diagram of the Experimental Setup

The transfer functions in the system are:

• Current Amplifier: from the voltage set by the DAC to the voltage across the Cercalo inductors \[ G_i = \begin{bmatrix} G_{i,h} & 0 \\ 0 & G_{i,v} \end{bmatrix} \]
• Voltage Amplifier: from the voltage across the Cercalo inductors to the measured voltage \[ G_a = \begin{bmatrix} G_{a,h} & 0 \\ 0 & G_{a,v} \end{bmatrix} \]
• Cercalo: Transfer function from the Voltage across the cercalo inductors to the 4 quadrant measurement \[ G_c = \begin{bmatrix} G_{\frac{V_{p,h}}{\tilde{U}_{c,h}}} & G_{\frac{V_{p,h}}{\tilde{U}_{c,v}}} \\ G_{\frac{V_{p,v}}{\tilde{U}_{c,h}}} & G_{\frac{V_{p,v}}{\tilde{U}_{c,v}}} \end{bmatrix} \]
• Newport Transfer function from the command signal of the Newport to the 4 quadrant measurement \[ G_n = \begin{bmatrix} G_{\frac{V_{p,h}}{U_{n,h}}} & G_{\frac{V_{p,h}}{U_{n,v}}} \\ G_{\frac{V_{p,v}}{U_{n,h}}} & G_{\frac{V_{n,v}}{U_{n,v}}} \end{bmatrix} \]

The block diagram with each transfer function is shown in Fig. fig:block_diagram.

Block Diagram of the Experimental Setup with detailed dynamics

Amplifier for the Cercalo

The value of the resistor in series with the buffer have been measured for both axis.

• $R_h = 41 \Omega$
• $L_{c,h} = 0.1 mH$
• $R_{c,h} = 9.3 \Omega$
• $R_v = 41 \Omega$
• $L_{c,v} = 0.1 mH$
• $R_{c,v} = 8.3 \Omega$

We want to find the transfer function from $U_c$ to $V_L$ and from $U_c$ to $i_c$.

We have that:

\begin{align*} V_C &= R_c i + L_c s i \\ U_c &= (R + R_c) i + L_c s i \end{align*}

Thus:

\begin{align} \frac{i_c}{U_c} &= \frac{1}{(R + R_c) + L_c s} \\ &= \frac{G_0}{1 + s/\omega_0} \end{align}

with

• $G_{0,i} = \frac{1}{R + R_c}$
• $\omega_0 = \frac{R + R_c}{L_c}$

And

\begin{align} \frac{V_c}{U_c} &= \frac{R_c + L_c s}{(R + R_c) + L_c s} \\ &= \frac{\frac{R_c}{R + R_c} + \frac{L_c}{R + R_c} s}{1 + \frac{L_c}{R + R_c} s} \\ &= \frac{G_0 + s/\omega_0}{1 + s/\omega_0} \\ \end{align}

with

• $G_0 = \frac{R_c}{R + R_c}$
• $\omega_0 = \frac{R + R_c}{L_c}$

Let's verify that the electrical circuit behaves as a constant current amplifier in the frequency band of interest.

  Rh = 41; % [Ohm]
  Lch = 0.1e-3; % [H]
  Rch = 9.3; % [Ohm]

  Rv = 41; % [Ohm]
  Lcv = 0.1e-3; % [H]
  Rcv = 8.3; % [Ohm]

  Gih = 1/(Rh + Rch + Lch * s);
  Gvh = (Rch + Lch * s)/(Rh + Rch + Lch * s);

  Giv = 1/(Rv + Rcv + Lcv * s);
  Gvv = (Rcv + Lcv * s)/(Rv + Rcv + Lcv * s);

  Gih0 = freqresp(Gih, 0);
  Gvh0 = freqresp(Gvh, 0);
  Giv0 = freqresp(Giv, 0);
  Gvv0 = freqresp(Gvv, 0);

  <<plt-matlab>>

Transfer function for the current amplifier (png, pdf)

The current amplifier has a constant gain over all the frequency band of interest. \[ G_i(s) = \begin{bmatrix} 0.02 & 0 \\ 0 & 0.02 \end{bmatrix}\quad \left[\frac{A}{V}\right] \] \[ G_a(s) = \begin{bmatrix} 0.185 & 0 \\ 0 & 0.168 \end{bmatrix} \left[\frac{V}{V}\right] \]

Cercalo

From the Cercalo documentation, we have the parameters shown on table tab:cercalo_parameters.

	Maximum Stroke [deg]	Resonance Frequency [Hz]	DC Gain [mA/deg]	Gain at resonance [deg/V]	RC Resistance [Ohm]
AX1 (Horizontal)	5	411.13	28.4	382.9	9.41
AX2 (Vertical)	5	252.5	35.2	350.4

Cercalo Parameters

The Inductance and DC resistance of the two axis of the Cercalo have been measured:

• $L_{c,h} = 0.1\ \text{mH}$
• $L_{c,v} = 0.1\ \text{mH}$
• $R_{c,h} = 9.3\ \Omega$
• $R_{c,v} = 8.3\ \Omega$

Let's first consider the horizontal direction and we try to model the Cercalo by a spring/mass/damper system (Fig. fig:mech_cercalo).

1 degree-of-freedom model of the Cercalo

The equation of motion is:

\begin{align*} \frac{x}{F} &= \frac{1}{k + c s + m s^2} \\ &= \frac{G_0}{1 + 2 \xi \frac{s}{\omega_0} + \frac{s^2}{\omega_0^2}} \end{align*}

with:

• $G_0 = 1/k$ is the gain at DC in rad/N
• $\xi = \frac{c}{2 \sqrt{km}}$ is the damping ratio of the system
• $\omega_0 = \sqrt{\frac{k}{m}}$ is the resonance frequency in rad

The force $F$ applied to the mass is proportional to the current $I$ flowing through the voice coils: \[ \frac{F}{I} = \alpha \] with $\alpha$ is in $N/A$ and is to be determined.

The current $I$ is also proportional to the voltage at the output of the buffer:

\begin{align*} \frac{I_c}{U_c} &= \frac{1}{(R + R_c) + L_c s} \\ &\approx 0.02 \left[ \frac{A}{V} \right] \end{align*}

Let's try to determine the equivalent mass and spring values. From table tab:cercalo_parameters, for the horizontal direction: \[ \left| \frac{x}{I} \right|(0) = \left| \alpha \frac{x}{F} \right|(0) = 28.4\ \frac{mA}{deg} = 1.63\ \frac{A}{rad} \]

So: \[ \alpha \frac{1}{k} = 1.63 \Longleftrightarrow k = \frac{\alpha}{1.63} \left[\frac{N}{rad}\right] \]

We also know the resonance frequency: \[ \omega_0 = 411.1\ \text{Hz} = 2583\ \frac{rad}{s} \]

And the gain at resonance:

\begin{align*} \left| \frac{x}{U_c} \right|(j\omega_0) &= \left| 0.02 \frac{x}{I_c} \right| (j\omega_0) \\ &= \left| 0.02 \alpha \frac{x}{F} \right| (j\omega_0) \\ &= 0.02 \alpha \frac{1/k}{2\xi} \\ &= 282.9\ \left[\frac{deg}{V}\right] \\ &= 4.938\ \left[\frac{rad}{V}\right] \end{align*}

Thus:

\begin{align*} & \frac{\alpha}{2 \xi k} = 245 \\ \Leftrightarrow & \frac{1.63}{2 \xi} = 245 \\ \Leftrightarrow & \xi = 0.0033 \\ \Leftrightarrow & \xi = 0.33 \% \end{align*}

\begin{align*} G_0 &= \frac{1.63}{\alpha}\ \frac{rad}{N} \\ \xi &= 0.0033 \\ \omega_0 &= 2583\ \frac{rad}{s} \end{align*}

and in terms of the physical properties:

\begin{align*} k &= \frac{\alpha}{1.63}\ \frac{N}{rad} \\ \xi &= 0.0033 \\ m &= \frac{\alpha}{1.1 \cdot 10^7}\ \frac{kg}{m^2} \end{align*}

Thus, we have to determine $\alpha$. This can be done experimentally by determining the gain at DC or at resonance of the system. For that, we need to know the angle of the mirror, thus we need to calibrate the photo-diodes. This will be done using the Newport.

Optical Setup

4 quadrant Diode

The front view of the 4 quadrant photo-diode is shown in Fig. fig:4qd_naming.

Front view of the 4QD

Each of the photo-diode is amplified using a 4-channel amplifier as shown in Fig. fig:4qd_amplifier.

Wiring of the amplifier. The amplifier is located on the bottom right of the board

ADC/DAC

Let's compute the theoretical noise of the ADC/DAC.

\begin{align*} \Delta V &= 20 V \\ n &= 16bits \\ q &= \Delta V/2^n = 305 \mu V \\ f_N &= 10kHz \\ \Gamma_n &= \frac{q^2}{12 f_N} = 7.76 \cdot 10^{-13} \frac{V^2}{Hz} \end{align*}

with $\Delta V$ the total range of the ADC, $n$ its number of bits, $q$ the quantization, $f_N$ the sampling frequency and $\Gamma_n$ its theoretical Power Spectral Density.

Identification

<<sec:plant_identification>>

Introduction   ignore

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.

Signals

Signal	Name	Unit
Voltage Sent to Cercalo - Horizontal	Uch	[V]
Voltage Sent to Cercalo - Vertical	Ucv	[V]
Voltage Sent to Newport - Horizontal	Unh	[V]
Voltage Sent to Newport - Vertical	Unv	[V]
4Q Photodiode Measurement - Horizontal	Vph	[V]
4Q Photodiode Measurement - Vertical	Vpv	[V]
Measured Voltage across the Inductance - Horizontal	Vch	[V]
Measured Voltage across the Inductance - Vertical	Vcv	[V]
Newport Metrology - Horizontal	Vnh	[V]
Newport Metrology - Vertical	Vnv	[V]
Attocube Measurement	Va	[m]

Huddle Test

We load the data taken during the Huddle Test.

  load('mat/data_test.mat', ...
       't', 'Uch', 'Ucv', ...
       'Unh', 'Unv', ...
       'Vph', 'Vpv', ...
       'Vch', 'Vcv', ...
       'Vnh', 'Vnv', ...
       'Va');

We remove the first second of data where everything is settling down.

  t0 = 1;

  Uch(t<t0) = [];
  Ucv(t<t0) = [];
  Unh(t<t0) = [];
  Unv(t<t0) = [];
  Vph(t<t0) = [];
  Vpv(t<t0) = [];
  Vch(t<t0) = [];
  Vcv(t<t0) = [];
  Vnh(t<t0) = [];
  Vnv(t<t0) = [];
  Va(t<t0)  = [];
  t(t<t0)   = [];

  t = t - t(1); % We start at t=0

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

  [psd_Vph, f] = pwelch(Vph, hanning(ceil(1*fs)), [], [], fs);
  [psd_Vpv, ~] = pwelch(Vpv, hanning(ceil(1*fs)), [], [], fs);

  figure;
  hold on;
  plot(f, sqrt(psd_Vph), 'DisplayName', '$\Gamma_{Vp_h}$');
  plot(f, sqrt(psd_Vpv), 'DisplayName', '$\Gamma_{Vp_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_Vch, f] = pwelch(Vch, hanning(ceil(1*fs)), [], [], fs);
  [psd_Vcv, ~] = pwelch(Vcv, hanning(ceil(1*fs)), [], [], fs);

  figure;
  hold on;
  plot(f, sqrt(psd_Vch), 'DisplayName', '$\Gamma_{Vc_h}$');
  plot(f, sqrt(psd_Vcv), 'DisplayName', '$\Gamma_{Vc_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', 'Uch', 'Ucv', ...
            'Unh', 'Unv', ...
            'Vph', 'Vpv', ...
            'Vch', 'Vcv', ...
            'Vnh', 'Vnv', ...
            'Va');

  uv = load('mat/data_uh.mat', ...
            't', 'Uch', 'Ucv', ...
            'Unh', 'Unv', ...
            'Vph', 'Vpv', ...
            'Vch', 'Vcv', ...
            'Vnh', 'Vnv', ...
            'Va');

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

  t0 = 1;

  uh.Uch(t<t0) = [];
  uh.Ucv(t<t0) = [];
  uh.Unh(t<t0) = [];
  uh.Unv(t<t0) = [];
  uh.Vph(t<t0) = [];
  uh.Vpv(t<t0) = [];
  uh.Vch(t<t0) = [];
  uh.Vcv(t<t0) = [];
  uh.Vnh(t<t0) = [];
  uh.Vnv(t<t0) = [];
  uh.Va(t<t0)  = [];
  uh.t(t<t0)   = [];

  uh.t = uh.t - uh.t(1); % We start at t=0

  t0 = 1;

  uv.Uch(t<t0) = [];
  uv.Ucv(t<t0) = [];
  uv.Unh(t<t0) = [];
  uv.Unv(t<t0) = [];
  uv.Vph(t<t0) = [];
  uv.Vpv(t<t0) = [];
  uv.Vch(t<t0) = [];
  uv.Vcv(t<t0) = [];
  uv.Vnh(t<t0) = [];
  uv.Vnv(t<t0) = [];
  uv.Va(t<t0)  = [];
  uv.t(t<t0)   = [];

  uv.t = uv.t - uv.t(1); % We start at t=0

  <<plt-matlab>>

Identification signals when exciting the horizontal direction (png, pdf)

  <<plt-matlab>>

Identification signals when exciting in the vertical direction (png, pdf)

Estimation of the Frequency Response Function Matrix

  win = hanning(ceil(1*fs));

We compute an estimate of the transfer functions.

  [tf_Uch_Vph, f] = tfestimate(uh.Uch, uh.Vph, win, [], [], fs);
  [tf_Uch_Vpv, ~] = tfestimate(uh.Uch, uh.Vpv, win, [], [], fs);
  [tf_Ucv_Vph, ~] = tfestimate(uv.Ucv, uv.Vph, win, [], [], fs);
  [tf_Ucv_Vpv, ~] = tfestimate(uv.Ucv, uv.Vpv, win, [], [], fs);

  <<plt-matlab>>

Frequency Response Matrix (png, pdf)

Coherence

  [coh_Uch_Vph, f] = mscohere(uh.Uch, uh.Vph, win, [], [], fs);
  [coh_Uch_Vpv, ~] = mscohere(uh.Uch, uh.Vpv, win, [], [], fs);
  [coh_Ucv_Vph, ~] = mscohere(uv.Ucv, uv.Vph, win, [], [], fs);
  [coh_Ucv_Vpv, ~] = mscohere(uv.Ucv, uv.Vpv, win, [], [], 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_uh = [410]; % [Hz]
  freqs_res_uv = [250]; % [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;

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

  xi = 0.001; % Approximate modal damping
  poles_uh = [2*pi*freqs_res_uh*(xi + 1i), 2*pi*freqs_res_uh*(xi - 1i)];
  poles_uv = [2*pi*freqs_res_uv*(xi + 1i), 2*pi*freqs_res_uv*(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_Uch_Vph = coh_Uch_Vph';
  weight_Uch_Vpv = coh_Uch_Vpv';
  weight_Ucv_Vph = coh_Ucv_Vph';
  weight_Ucv_Vpv = coh_Ucv_Vpv';

  alpha = 0.1;

  for freq_i = 1:length(freqs_res)
    weight_Uch_Vph(f>(1-alpha)*freqs_res_uh(freq_i) & f<(1 + alpha)*freqs_res_uh(freq_i)) = 10;
    weight_Uch_Vpv(f>(1-alpha)*freqs_res_uh(freq_i) & f<(1 + alpha)*freqs_res_uh(freq_i)) = 10;
    weight_Ucv_Vph(f>(1-alpha)*freqs_res_uv(freq_i) & f<(1 + alpha)*freqs_res_uv(freq_i)) = 10;
    weight_Ucv_Vpv(f>(1-alpha)*freqs_res_uv(freq_i) & f<(1 + alpha)*freqs_res_uv(freq_i)) = 10;
  end

Ignore data above some frequency.

  weight_Uch_Vph(f>1000) = 0;
  weight_Uch_Vpv(f>1000) = 0;
  weight_Ucv_Vph(f>1000) = 0;
  weight_Ucv_Vpv(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_Uch_Vph, poles, ~, fit_Uch_Vph] = vectfit3(tf_Uch_Vph.', 1i*2*pi*f, poles_uh, weight_Uch_Vph, opts);
  end
  for iter = 1:Niter
    [SER_Uch_Vpv, poles, ~, fit_Uch_Vpv] = vectfit3(tf_Uch_Vpv.', 1i*2*pi*f, poles_uh, weight_Uch_Vpv, opts);
  end
  for iter = 1:Niter
    [SER_Ucv_Vph, poles, ~, fit_Ucv_Vph] = vectfit3(tf_Ucv_Vph.', 1i*2*pi*f, poles_uv, weight_Ucv_Vph, opts);
  end
  for iter = 1:Niter
    [SER_Ucv_Vpv, poles, ~, fit_Ucv_Vpv] = vectfit3(tf_Ucv_Vpv.', 1i*2*pi*f, poles_uv, weight_Ucv_Vpv, opts);
  end

  <<plt-matlab>>

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

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

Plant Scaling

• measured noise
• expected perturbations
• maximum input usage
• maximum wanted error

Plant Analysis

Load Plant

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

RGA-Number

  freqs = logspace(2, 4, 1000);
  G_resp = freqresp(G, freqs, 'Hz');
  A = zeros(size(G_resp));
  RGAnum = zeros(1, length(freqs));

  for i = 1:length(freqs)
    A(:, :, i) = G_resp(:, :, i).*inv(G_resp(:, :, i))';
    RGAnum(i) = sum(sum(abs(A(:, :, i)-eye(2))));
  end
  % RGA = G0.*inv(G0)';

  figure;
  plot(freqs, RGAnum);
  set(gca, 'xscale', 'log');

  U = zeros(2, 2, length(freqs));
  S = zeros(2, 2, length(freqs))
  V = zeros(2, 2, length(freqs));

  for i = 1:length(freqs)
    [Ui, Si, Vi] = svd(G_resp(:, :, i));
    U(:, :, i) = Ui;
    S(:, :, i) = Si;
    V(:, :, i) = Vi;
  end

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

Control Design

Measurement of the non-repeatability