Sercalo Test Bench
Table of ContentsClose
This report is also available as a pdf.
1 Introduction
1.1 Block Diagram
The block diagram of the setup to be controlled is shown in Fig. 1.
Figure 1: Block Diagram of the Experimental Setup
The transfer functions in the system are:
- Current Amplifier: from the voltage set by the DAC to the current going to the Sercalo’s inductors
- Impedance of the Sercalo that converts the current going to the sercalo to the voltage across the sercalo:
- Voltage Amplifier: from the voltage across the Sercalo inductors to the measured voltage
- Sercalo: Transfer function from the current going through the sercalo inductors to the 4 quadrant measurement
- Newport Transfer function from the command signal of the Newport to the 4 quadrant measurement
- 4 Quadrant Diode: the gain of the 4 quadrant diode in [V/rad] is inverse in order to obtain the physical angle of the beam
The block diagram with each transfer function is shown in Fig. 2.
Figure 2: Block Diagram of the Experimental Setup with detailed dynamics
1.2 Sercalo
From the Sercalo documentation, we have the parameters shown on table 1.
Max. Stroke | Res. Freq. | DC Gain | Gain at res. | RC Res. | |
---|---|---|---|---|---|
[deg] | [Hz] | [mA/deg] | [deg/V] | [Ohm] | |
AX1 (Horizontal) | 5 | 411.13 | 28.4 | 382.9 | 9.41 |
AX2 (Vertical) | 5 | 252.5 | 35.2 | 350.4 |
The Inductance and DC resistance of the two axis of the Sercalo have been measured:
Let’s first consider the horizontal direction and we try to model the Sercalo by a spring/mass/damper system (Fig. 3).
Figure 3: 1 degree-of-freedom model of the Sercalo
The equation of motion is:
with:
is the gain at DC in rad/N is the damping ratio of the system is the resonance frequency in rad
The force
The current
Let’s try to determine the equivalent mass and spring values.
From table 1, for the horizontal direction:
So:
We also know the resonance frequency:
And the gain at resonance:
Thus:
Important
and in terms of the physical properties:
Thus, we have to determine
1.3 Optical Setup
1.4 Newport
1.5 4 quadrant Diode
1.6 ADC/DAC
Let’s compute the theoretical noise of the ADC/DAC.
with
2 Identification of the system dynamics
In this section, we seek to identify all the blocks as shown in Fig. 1.
Signal | Name | Unit |
---|---|---|
Voltage Sent to Sercalo - Horizontal | Uch |
[V] |
Voltage Sent to Sercalo - 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] |
2.1 Calibration of the 4 Quadrant Diode
Prior to any dynamic identification, we would like to be able to determine the meaning of the 4 quadrant diode measurement. For instance, instead of obtaining transfer function in [V/V] from the input of the sercalo to the measurement voltage of the 4QD, we would like to obtain the transfer function in [rad/V]. This will give insight to physical interpretation.
To calibrate the 4 quadrant photo-diode, we can use the metrology included in the Newport. We can choose precisely the angle of the Newport mirror and see what is the value measured by the 4 Quadrant Diode. We then should be able to obtain the “gain” of the 4QD in [V/rad].
2.1.1 Input / Output data
The identification data is loaded
Copyuh = load('mat/data_cal_pd_h.mat', 't', 'Vph', 'Vpv', 'Vnh'); uv = load('mat/data_cal_pd_v.mat', 't', 'Vph', 'Vpv', 'Vnv');
We remove the first seconds where the Sercalo is turned on.
Copyt0 = 1; uh.Vph(uh.t<t0) = []; uh.Vpv(uh.t<t0) = []; uh.Vnh(uh.t<t0) = []; uh.t(uh.t<t0) = []; uh.t = uh.t - uh.t(1); % We start at t=0 t0 = 1; uv.Vph(uv.t<t0) = []; uv.Vpv(uv.t<t0) = []; uv.Vnv(uv.t<t0) = []; uv.t(uv.t<t0) = []; uv.t = uv.t - uv.t(1); % We start at t=0
2.1.2 Linear Regression to obtain the gain of the 4QD
We plot the angle of mirror
Gain of the Newport metrology in [rad/V].
Copygn0 = 2.62e-3;
The angular displacement of the beam is twice the angular displacement of the Newport mirror.
We do a linear regression
is the measured voltage of the 4QD in [V] is the beam angle (twice the mirror angle) in [rad] is the identified gain of the 4QD in [rad/V]
The linear regression is shown in Fig. 10.
Copybh = [ones(size(uh.Vnh)) 2*gn0*uh.Vnh]\uh.Vph; bv = [ones(size(uv.Vnv)) 2*gn0*uv.Vnv]\uv.Vpv;
Thus, we obtain the “gain of the 4 quadrant photo-diode as shown on table 2.
Horizontal [V/rad] | Vertical [V/rad] |
---|---|
-31.0 | 36.3 |
CopyGd = tf([bh(2) 0 ; 0 bv(2)]);
We obtain:
2.2 Identification of the Sercalo Impedance, Current Amplifier and Voltage Amplifier dynamics
We wish here to determine
We ignore the electro-mechanical coupling.
2.2.1 Electrical Schematic
The schematic of the electrical circuit used to drive the Sercalo is shown in Fig. 11.
Figure 11: Current Amplifier Schematic
The elements are:
: the voltage generated by the DAC- BUF: is a unity-gain open-loop buffer that allows to increase the output current
: a chosen resistor that will determine the gain of the current amplifier : inductor present in the Sercalo : resistance of the inductor : voltage measured across the Sercalo’s inductor : amplified voltage measured across the Sercalo’s inductor is the current going through the Sercalo’s inductor
The values of the components have been measured for the horizontal and vertical directions:
Let’s first determine the transfer function from
We have that:
Thus:
with
Now, determine the transfer function from
Finally, the transfer function of the voltage amplifier
with
is the gain 1000 (60dB) is the cut-off frequency of the voltage amplifier set to 1000Hz
2.2.2 Theoretical Transfer Functions
The values of the components in the current amplifier have been measured.
CopyRh = 41; % [Ohm] Lch = 0.1e-3; % [H] Rch = 9.3; % [Ohm] Rv = 41; % [Ohm] Lcv = 0.1e-3; % [H] Rcv = 8.3; % [Ohm]
CopyGi = blkdiag(1/(Rh + Rch + Lch * s), 1/(Rv + Rcv + Lcv * s)); Zc = blkdiag(Rch+Lch*s, Rcv+Lcv*s); Ga = blkdiag(1000/(1 + s/2/pi/1000), 1000/(1 + s/2/pi/1000));
Important
Over the frequency band of interest, the current amplifier transfer function
CopyGi = tf(blkdiag(1/(Rh + Rch), 1/(Rv + Rcv))); Zc = tf(blkdiag(Rch, Rcv));
2.2.3 Identified Transfer Functions
Noise is generated using the DAC (
The identification data is loaded.
Copyuh = load('mat/data_uch.mat', 't', 'Uch', 'Vch'); uv = load('mat/data_ucv.mat', 't', 'Ucv', 'Vcv');
We remove the first seconds where the Sercalo is turned on.
Copywin = hanning(ceil(1*fs)); [GaZcGi_h, f] = tfestimate(uh.Uch, uh.Vch, win, [], [], fs); [GaZcGi_v, ~] = tfestimate(uv.Ucv, uv.Vcv, win, [], [], fs);
There is a gain mismatch, that is probably due to bad identification of the inductance and resistance measurement of the sercalo inductors.
Thus, we suppose
CopyGi_resp_h = abs(GaZcGi_h)./squeeze(abs(freqresp(Ga(1,1)*Zc(1,1), f, 'Hz'))); Gi_resp_v = abs(GaZcGi_v)./squeeze(abs(freqresp(Ga(2,2)*Zc(2,2), f, 'Hz'))); Gi = tf(blkdiag(mean(Gi_resp_h(f>20 & f<200)), mean(Gi_resp_v(f>20 & f<200))));
Finally, we have the following transfer functions:
ans = filepath; if ischar(ans), fid = fopen('/tmp/babel-ZKMGJu/matlab-FA7h5L', 'w'); fprintf(fid, '%s\n', ans); fclose(fid); else, dlmwrite('/tmp/babel-ZKMGJu/matlab-FA7h5L', ans, '\t') end 'org_babel_eoe' Gi,Zc,Ga 'org_babel_eoe' ans = filepath; if ischar(ans), fid = fopen('/tmp/babel-ZKMGJu/matlab-FA7h5L', 'w'); fprintf(fid, '%s\n', ans); fclose(fid); else, dlmwrite('/tmp/babel-ZKMGJu/matlab-FA7h5L', ans, '\t') end 'org_babel_eoe' ans = 'org_babel_eoe' Gi,Zc,Ga Gi = From input 1 to output... 1: 0.01275 2: 0 From input 2 to output... 1: 0 2: 0.01382 Static gain. Zc = From input 1 to output... 1: 9.3 2: 0 From input 2 to output... 1: 0 2: 8.3 Static gain. Ga = From input 1 to output... 6.2832e+06 1: ---------- (s+6283) 2: 0 From input 2 to output... 1: 0 6.2832e+06 2: ---------- (s+6283) Continuous-time zero/pole/gain model.
2.3 Identification of the Sercalo Dynamics
We now wish to identify the dynamics of the Sercalo identified by
To do so, we inject some noise at the input of the current amplifier
The transfer function obtained will be
2.3.1 Input / Output data
The identification data is loaded
Copyuh = load('mat/data_uch.mat', 't', 'Uch', 'Vph', 'Vpv'); uv = load('mat/data_ucv.mat', 't', 'Ucv', 'Vph', 'Vpv');
We remove the first seconds where the Sercalo is turned on.
Copyt0 = 1; uh.Uch(uh.t<t0) = []; uh.Vph(uh.t<t0) = []; uh.Vpv(uh.t<t0) = []; uh.t(uh.t<t0) = []; uh.t = uh.t - uh.t(1); % We start at t=0 t0 = 1; uv.Ucv(uv.t<t0) = []; uv.Vph(uv.t<t0) = []; uv.Vpv(uv.t<t0) = []; uv.t(uv.t<t0) = []; uv.t = uv.t - uv.t(1); % We start at t=0
2.3.2 Coherence
The window used for the spectral analysis is an hanning
windows with temporal size equal to 1 second.
Copywin = hanning(ceil(1*fs));
Copy[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);
2.3.3 Estimation of the Frequency Response Function Matrix
We compute an estimate of the transfer functions.
Copy[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);
2.3.4 Time Delay
Now, we would like to remove the time delay included in the FRF prior to the model extraction.
Estimation of the time delay:
CopyTs_delay = Ts; % [s] G_delay = tf(1, 1, 'InputDelay', Ts_delay); G_delay_resp = squeeze(freqresp(G_delay, f, 'Hz'));
We then remove the time delay from the frequency response function.
Copytf_Uch_Vph = tf_Uch_Vph./G_delay_resp; tf_Uch_Vpv = tf_Uch_Vpv./G_delay_resp; tf_Ucv_Vph = tf_Ucv_Vph./G_delay_resp; tf_Ucv_Vpv = tf_Ucv_Vpv./G_delay_resp;
2.3.5 Extraction of a transfer function matrix
First we define the initial guess for the resonance frequencies and the weights associated.
Copyfreqs_res_uh = [410]; % [Hz] freqs_res_uv = [250]; % [Hz]
We then make an initial guess on the complex values of the poles.
Copyxi = 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. Finally, we ignore data above some frequency.
Copyweight_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_uh) 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 weight_Uch_Vph(f>1000) = 0; weight_Uch_Vpv(f>1000) = 0; weight_Ucv_Vph(f>1000) = 0; weight_Ucv_Vpv(f>1000) = 0;
The weights are shown in Fig. 20.
When we set some options for vfit3
.
Copyopts = 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.
CopyNiter = 5;
An we run the vectfit3
algorithm.
Copyfor 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
And finally, we create the identified
CopyG_Uch_Vph = tf(minreal(ss(full(SER_Uch_Vph.A),SER_Uch_Vph.B,SER_Uch_Vph.C,SER_Uch_Vph.D))); G_Ucv_Vph = tf(minreal(ss(full(SER_Ucv_Vph.A),SER_Ucv_Vph.B,SER_Ucv_Vph.C,SER_Ucv_Vph.D))); G_Uch_Vpv = tf(minreal(ss(full(SER_Uch_Vpv.A),SER_Uch_Vpv.B,SER_Uch_Vpv.C,SER_Uch_Vpv.D))); G_Ucv_Vpv = tf(minreal(ss(full(SER_Ucv_Vpv.A),SER_Ucv_Vpv.B,SER_Ucv_Vpv.C,SER_Ucv_Vpv.D))); Gc = [G_Uch_Vph, G_Ucv_Vph; G_Uch_Vpv, G_Ucv_Vpv]*inv(Gi);
2.4 Identification of the Newport Dynamics
We here identify the transfer function from a reference sent to the Newport
To do so, we inject noise to the Newport
2.4.1 Input / Output data
The identification data is loaded
Copyuh = load('mat/data_unh.mat', 't', 'Unh', 'Vph', 'Vpv'); uv = load('mat/data_unv.mat', 't', 'Unv', 'Vph', 'Vpv');
We remove the first seconds where the Sercalo is turned on.
Copyt0 = 3; uh.Unh(uh.t<t0) = []; uh.Vph(uh.t<t0) = []; uh.Vpv(uh.t<t0) = []; uh.t(uh.t<t0) = []; uh.t = uh.t - uh.t(1); % We start at t=0 t0 = 1.5; uv.Unv(uv.t<t0) = []; uv.Vph(uv.t<t0) = []; uv.Vpv(uv.t<t0) = []; uv.t(uv.t<t0) = []; uv.t = uv.t - uv.t(1); % We start at t=0
2.4.2 Coherence
The window used for the spectral analysis is an hanning
windows with temporal size equal to 1 second.
Copywin = hanning(ceil(1*fs));
Copy[coh_Unh_Vph, f] = mscohere(uh.Unh, uh.Vph, win, [], [], fs); [coh_Unh_Vpv, ~] = mscohere(uh.Unh, uh.Vpv, win, [], [], fs); [coh_Unv_Vph, ~] = mscohere(uv.Unv, uv.Vph, win, [], [], fs); [coh_Unv_Vpv, ~] = mscohere(uv.Unv, uv.Vpv, win, [], [], fs);
2.4.3 Estimation of the Frequency Response Function Matrix
We compute an estimate of the transfer functions.
Copy[tf_Unh_Vph, f] = tfestimate(uh.Unh, uh.Vph, win, [], [], fs); [tf_Unh_Vpv, ~] = tfestimate(uh.Unh, uh.Vpv, win, [], [], fs); [tf_Unv_Vph, ~] = tfestimate(uv.Unv, uv.Vph, win, [], [], fs); [tf_Unv_Vpv, ~] = tfestimate(uv.Unv, uv.Vpv, win, [], [], fs);
2.4.4 Time Delay
Now, we would like to remove the time delay included in the FRF prior to the model extraction.
Estimation of the time delay:
CopyTs_delay = 0.0005; % [s] G_delay = tf(1, 1, 'InputDelay', Ts_delay); G_delay_resp = squeeze(freqresp(G_delay, f, 'Hz'));
We then remove the time delay from the frequency response function.
2.4.5 Extraction of a transfer function matrix
From Fig. 26, it seems reasonable to model the Newport dynamics as diagonal and constant.
CopyGn = blkdiag(tf(mean(abs(tf_Unh_Vph(f>10 & f<100)))), tf(mean(abs(tf_Unv_Vpv(f>10 & f<100)))));
2.5 Full System
We now have identified:
We name the input and output of each transfer function:
CopyGi.InputName = {'Uch', 'Ucv'}; Gi.OutputName = {'Ich', 'Icv'}; Zc.InputName = {'Ich', 'Icv'}; Zc.OutputName = {'Vtch', 'Vtcv'}; Ga.InputName = {'Vtch', 'Vtcv'}; Ga.OutputName = {'Vch', 'Vcv'}; Gc.InputName = {'Ich', 'Icv'}; Gc.OutputName = {'Vpch', 'Vpcv'}; Gn.InputName = {'Unh', 'Unv'}; Gn.OutputName = {'Vpnh', 'Vpnv'}; Gd.InputName = {'Rh', 'Rv'}; Gd.OutputName = {'Vph', 'Vpv'};
CopySh = sumblk('Vph = Vpch + Vpnh'); Sv = sumblk('Vpv = Vpcv + Vpnv');
Copyinputs = {'Uch', 'Ucv', 'Unh', 'Unv'}; outputs = {'Vch', 'Vcv', 'Ich', 'Icv', 'Rh', 'Rv', 'Vph', 'Vpv'}; sys = connect(Gi, Zc, Ga, Gc, Gn, inv(Gd), Sh, Sv, inputs, outputs);
The file mat/plant.mat
is accessible here.
Copysave('mat/plant.mat', 'sys', 'Gi', 'Zc', 'Ga', 'Gc', 'Gn', 'Gd');
3 Huddle Test
The goal is to determine the noise of the photodiodes as well as the noise of the Attocube interferometer.
Multiple measurements are done with different experimental configuration as follow:
Number | OL/CL | Compensation Unit | Aluminum |
---|---|---|---|
1 | Open Loop | ||
2 | Open Loop | Compensation Unit | |
3 | Closed Loop | Compensation Unit | |
4 | Open Loop | Compensation Unit | Aluminum |
5 | Closed Loop | Compensation Unit | Aluminum |
3.1 Load Data
Copyht_1 = load('./mat/data_huddle_test_1.mat', 't', 'Vph', 'Vpv', 'Va'); ht_2 = load('./mat/data_huddle_test_2.mat', 't', 'Vph', 'Vpv', 'Va'); ht_3 = load('./mat/data_huddle_test_3.mat', 't', 'Uch', 'Ucv', 'Vph', 'Vpv', 'Va'); ht_4 = load('./mat/data_huddle_test_4.mat', 't', 'Vph', 'Vpv', 'Va'); % ht_5 = load('./mat/data_huddle_test_5.mat', 't', 'Uch', 'Ucv', 'Vph', 'Vpv', 'Va');
Copyfs = 1e4;
3.2 Pre-processing
Copyt0 = 1; % [s] tend = 100; % [s] ht_s = {ht_1 ht_2 ht_3 ht_4} for i = 1:length(ht_s) ht_s{i}.Vph(ht_s{i}.t<t0) = []; ht_s{i}.Vpv(ht_s{i}.t<t0) = []; ht_s{i}.Va(ht_s{i}.t<t0) = []; ht_s{i}.t(ht_s{i}.t<t0) = []; ht_s{i}.t = ht_s{i}.t - ht_s{i}.t(1); % We start at t=0 ht_s{i}.Vph(tend*fs+1:end) = []; ht_s{i}.Vpv(tend*fs+1:end) = []; ht_s{i}.Va(tend*fs+1:end) = []; ht_s{i}.t(tend*fs+1:end) = []; ht_s{i}.Va = ht_s{i}.Va - mean(ht_s{i}.Va); end ht_1 = ht_s{1}; ht_2 = ht_s{2}; ht_3 = ht_s{3}; ht_4 = ht_s{4};
3.3 Filter data with low pass filter
We filter the data with a first order low pass filter with a crossover frequency of
Copyw0 = 50; % [Hz] G_lpf = 1/(1 + s/2/pi/w0); ht_1.Vaf = lsim(G_lpf, ht_1.Va, ht_1.t); ht_2.Vaf = lsim(G_lpf, ht_2.Va, ht_2.t); ht_3.Vaf = lsim(G_lpf, ht_3.Va, ht_3.t); ht_4.Vaf = lsim(G_lpf, ht_4.Va, ht_4.t);
3.4 Time domain plots
3.5 Power Spectral Density
Copywin = hanning(ceil(1*fs));
Copy[psd_Va1, f] = pwelch(ht_1.Va, win, [], [], fs); [psd_Va2, ~] = pwelch(ht_2.Va, win, [], [], fs); [psd_Va3, ~] = pwelch(ht_3.Va, win, [], [], fs); [psd_Va4, ~] = pwelch(ht_4.Va, win, [], [], fs);
Copy[psd_Vph1, ~] = pwelch(ht_1.Vph, win, [], [], fs); [psd_Vph2, ~] = pwelch(ht_2.Vph, win, [], [], fs); [psd_Vph3, ~] = pwelch(ht_3.Vph, win, [], [], fs); [psd_Vph4, ~] = pwelch(ht_4.Vph, win, [], [], fs);
Copy[psd_Vpv1, ~] = pwelch(ht_1.Vpv, win, [], [], fs); [psd_Vpv2, ~] = pwelch(ht_2.Vpv, win, [], [], fs); [psd_Vpv3, ~] = pwelch(ht_3.Vpv, win, [], [], fs); [psd_Vpv4, ~] = pwelch(ht_4.Vpv, win, [], [], fs);
3.6 Conclusion
The Attocube’s “Environmental Compensation Unit” does not have a significant effect on the stability of the measurement.
4 Budget Error
Goals:
- List all sources of error and compute their effects on the Attocube measurement
- Think about how to determine the value of the individual sources of error
- Sum all the sources of error and determine the limiting ones
Sources of error for the Attocube measurement:
- Beam non-perpendicularity to the concave mirror is linked to the non-perfect feedback loop:
- We have only finite gain / limited bandwidth so the Sercalo mirror angle will not be perfect
- The non-perpendicularity is measured by the 4QD and is used as the feedback signal, however this signal is noisy and even with infinite gain, this noise will be transmitted to the angle of the beam
- Sercalo/Newport unwanted translation perpendicular to its surface.
This can be due to:
- Non idealities in the mechanics of the Sercalo
- Temperature variations
- The reproducible part of the perpendicular translation with respect to the angle of the Sercalo can be taken into account and subtracted from the Attocube measurement
- Temperature variations of the metrology frame
- Change in the refractive air index in the beam path. This can be due to change of Temperature, Pressure and Humidity of the air in the beam path
Procedure:
- in section 4.1: We estimate the effect of an angle error of the Sercalo mirror on the Attocube measurement
- in section 4.2: The effect of perpendicular motion of the Newport and Sercalo mirrors on the Attocube measurement is determined.
- in section 4.3: We estimate the expected change of refractive index of the air in the beam path and the resulting Attocube measurement error
- in section 4.5: The feedback system using the 4 quadrant diode and the Sercalo is studied. Sensor noise, actuator noise and their effects on the control error is discussed.
4.1 Effect of the Sercalo angle error on the measured distance by the Attocube
To simplify, we suppose that the Newport mirror is a flat mirror (instead of a concave one).
The geometry of the setup is shown in Fig. 33 where:
is the reference surface of the Attocube is the point where the beam first hits the Sercalo mirror is the point where the beam first hits the Newport mirror is the angle error from its ideal 45 degrees
We define the angle error range
Copythetas_c = logspace(-7, -4, 100); % [rad]
The geometrical parameters of the setup are defined below.
CopyH = 0.05; % [m] L = 0.05; % [m]
Figure 33: Schematic of the geometry used to evaluate the effect of
The nominal points
CopyO = [-L, 0]; S = [0, 0]; X = [0, H];
Thus, the initial path length
Copypath_nominal = norm(S-O) + norm(X-S) + norm(S-X) + norm(O-S);
We now compute the new path length when there is an error angle
Copypath_length = zeros(size(thetas_c)); for i = 1:length(thetas_c) theta_c = thetas_c(i); Y = [H*tan(2*theta_c), H]; M = 2*H/(tan(pi/4-theta_c)+1/tan(2*theta_c))*[1, tan(pi/4-theta_c)]; T = [-L, M(2)+(L+M(1))*tan(4*theta_c)]; path_length(i) = norm(S-O) + norm(Y-S) + norm(M-Y) + norm(T-M); end
We then compute the distance error and we plot it as a function of the Sercalo angle error (Fig. 34).
Copypath_error = path_length - path_nominal;
Figure 34: Effect of an angle error of the Sercalo on the distance error measured by the Attocube (png, pdf)
And we plot the beam path using Matlab for an high angle to verify that the code is working (Fig. 35).
Copytheta = 2*2*pi/360; % [rad] H = 0.05; % [m] L = 0.05; % [m] O = [-L, 0]; S = [0, 0]; X = [0, H]; Y = [H*tan(2*theta), H]; M = 2*H/(tan(pi/4-theta)+1/tan(2*theta))*[1, tan(pi/4-theta)]; T = [-L, M(2)+(L+M(1))*tan(4*theta)];
Important
Based on Fig. 34, we see that an angle error
with:
expressed in [m] in [rad]
Some example are shown in table 4.
The tracking error of the feedback system used to position the Sercalo mirror should thus be limited to few micro-meters.
Angle Error |
Distance measurement error |
---|---|
4.2 Unwanted motion of Sercalo/Newport mirrors perpendicular to its surface
From Figs 36 and 37, it is clear that perpendicular motions of the Sercalo mirror and of the Newport mirror have an impact on the measured distance by the Attocube interferometer.
More precisely, if the note:
the perpendicular motion of the Sercalo’s mirror the perpendicular motion of the Newport’s mirror
We have that:
Note here that
Figure 36: Effect of a Perpendicular motion of the Sercalo Mirror
Figure 37: Effect of a Perpendicular motion of the Newport Mirror
Important
The motion of the both Sercalo’s and Newport’s mirrors perpendicular to its surface is fully transmitted to the measured distance by the Attocube interferometer.
This motion can be measured and the repeatable part can be compensated. However, the non repeatability of this motion should be less than few nano-meters.
4.3 Change in refractive index of the air in the beam path
Three physical properties of the air makes change of the Attocube measurement:
- Temperature:
- Pressure:
- Humidity:
These physical properties should change relatively slowly, however, for a beam path of 100mm:
Air property Variations | Measurement error |
---|---|
100nm | |
27nm | |
10nm |
An Environmental Compensation Unit is used and can compensate for variations or air properties up to:
Air property Variations | Measurement error |
---|---|
Important
The total measurement error induced by air properties variations is then:
The beam path should be protected using aluminum to minimize the change in refractive index of the air in the beam path.
4.4 Thermal Expansion of the Metrology Frame
The material used for the metrology frame is Aluminum.
Its linear thermal expansion coefficient is
The distance between the Attocube head and the Attocube is approximatively equal to 5cm.
Thus, the temperature of the metrology frame should be kept constant to less than
4.5 Estimation of the Sercalo angle error due to Noise
In this section, we seek to estimate the angle error
Consider the block diagram in Fig. 38 with:
: represents the transfer function from a voltage applied by the Speedgoat DAC used for the Sercalo to the Beam angle : is the control law used
The signals are:
: is the “true” laser beam angle : is the measured beam angle ( ) : is the measurement noise of the laser beam angle using the 4 quadrant diode. It includes:- ADC noise
noise, Shot noise, Ambian noise, Intensity noise…
: is noise at the input of the Sercalo. It includes:- DAC noise of the speadgoat
: is disturbance on the angle of the beam. It includes:- Angle variations of the Newport mirror
Figure 38: Block Diagram of the Feedback system
4.5.1 Estimation of sources of noise and disturbances
Let’s estimate the values of
4.5.1.1 ADC Quantization Noise
The ADC quantization noise is:
with:
: is the range of the ADC in : is the number of ADC’s bits : is the sampling frequency in
For the ADC used:
Important
4.5.1.2 DAC Quantization Noise
The DAC quantization noise is:
with:
: is the range of the DAC in : is the number of DAC’s bits : is the sampling frequency in
For the DAC used:
Important
4.5.1.3 Noise of the Newport Mirror angle
Plus, we estimate the effect of DAC quantization noise on the angle error on the Newport mirror.
The gain of the Newport is:
If we integrate that to obtain an rms value:
Which is much less than the noise equivalent angle specified by Newport:
We expect the angle noise of the Newport mirror to be around
Important
If we suppose a white noise, the power spectral density of the beam angle due to the noise of the Newport mirror corresponds to:
4.5.1.4 Disturbances due the Newport Mirror Rotation
We will rotate the Newport mirror in order to simulate a displacement of the Sample:
- The angle range for the Newport mirror is
- The radius of the concave mirror is 200 mm
Figure 39: Rotation of the (concave) Newport mirror
If we suppose small angles, the corresponding beam deviation is:
4.5.2 Perfect Control
If the feedback is perfect, the Sercalo angle error will be equal to the 4 quadrant diode noise. Let’s estimate the 4QD noise in radians.
If we note
We should then have that the voltage of each quadrant is as large as possible.
Suppose here that
If we just consider the ADC noise:
- the ADC range is
with . - thus, the LSB corresponds to:
- this corresponds to an error
4.5.3 Error due to DAC noise used for the Sercalo
Copyload('./mat/plant.mat', 'Gi', 'Gc', 'Gd');
CopyG = inv(Gd)*Gc*Gi;
Dynamical estimation:
- ASD of DAC noise used for the Sercalo
- Multiply by transfer function from Sercalo voltage to angle estimation using the 4QD
Copyfreqs = logspace(1, 3, 1000); fs = 1e4; % ASD of the DAC voltage going to the Sercalo in [V/sqrt(Hz)] asd_uc = (20/2^16)/sqrt(12*fs)*ones(length(freqs), 1); % ASD of the measured angle by the QD in [rad/sqrt(Hz)] asd_theta = asd_uc.*abs(squeeze(freqresp(G(1,1), freqs, 'Hz'))); figure; loglog(freqs, asd_theta)
Then the corresponding ASD of the measured displacement by the interferometer is:
Copyasd_L = asd_theta*10^(-6); % [m/sqrt(Hz)]
And we integrate that to have the RMS value:
Copycps_L = 1/pi*cumtrapz(2*pi*freqs, (asd_L).^2);
The RMS value is:
Copysqrt(cps_L(end))
1.647e-11
Copyfigure; loglog(freqs, cps_L)
Let’s estimate the beam angle error corresponding to 1 LSB of the sercalo’s DAC. Gain of the Sercalo is approximatively 5 degrees for 10V. However the beam angle deviation is 4 times the angle deviation of the sercalo mirror, thus:
Copyd_alpha = 4*(20/2^16)*(5*pi/180)/10 % [rad]
1.0653e-05
This corresponds to a measurement error of the Attocube equals to (in [m])
Copy1e-6*d_alpha % [m]
1.0653e-11
Important
The DAC noise use for the Sercalo does not limit the performance of the system.
5 Plant Uncertainty
5.1 Coprime Factorization
Copyload('mat/plant.mat', 'sys', 'Gi', 'Zc', 'Ga', 'Gc', 'Gn', 'Gd');
Copy[fact, Ml, Nl] = lncf(Gc*Gi);
6 Plant Scaling
The goal is the scale the plant prior to control synthesis. This will simplify the choice of weighting functions and will yield useful insight on the controllability of the plant.
Value | Unit | Variable Name | |
---|---|---|---|
Expected perturbations | 1 | [V] | |
Maximum input usage | 10 | [V] | |
Maximum wanted error | 10 | [ |
|
Measured noise | 5 | [ |
6.1 Control Objective
The maximum expected stroke is
Thus, we require the sensitivity function at
In terms of loop gain, this is equivalent to:
6.2 General Configuration
The plant is put in a general configuration as shown in Fig. 40.
Figure 40: General Control Configuration
7 Plant Analysis
7.1 Load Plant
Copyload('mat/plant.mat', 'G');
7.2 RGA-Number
Copyfreqs = 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)';
Copyfigure; plot(freqs, RGAnum); set(gca, 'xscale', 'log');
CopyU = 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
7.3 Rotation Matrix
CopyG0 = freqresp(G, 0);
8 Active Damping
8.1 Load Plant
Copyload('mat/plant.mat', 'sys', 'Gi', 'Zc', 'Ga', 'Gc', 'Gn', 'Gd');
8.2 Integral Force Feedback
Copybode(sys({'Vch', 'Vcv'}, {'Uch', 'Ucv'}));
CopyKppf = blkdiag(-10000/s, tf(0)); Kppf.InputName = {'Vch', 'Vcv'}; Kppf.OutputName = {'Uch', 'Ucv'};
Copyinputs = {'Uch', 'Ucv', 'Unh', 'Unv'}; outputs = {'Ich', 'Icv', 'Rh', 'Rv', 'Vph', 'Vpv'}; sys_cl = connect(sys, Kppf, inputs, outputs); figure; bode(sys_cl({'Vph', 'Vpv'}, {'Uch', 'Ucv'}), sys({'Vph', 'Vpv'}, {'Uch', 'Ucv'}))
8.3 Conclusion
Active damping does not seems to be applicable here.
9 Decentralized Control of the Sercalo
9.1 Load Plant
Copyload('mat/plant.mat', 'sys', 'Gi', 'Zc', 'Ga', 'Gc', 'Gn', 'Gd');
9.2 Diagonal Controller
Using SISOTOOL
, a diagonal controller is designed.
The two SISO loop gains are shown in Fig. 41.
CopyKh = -0.25598*(s+112)*(s^2 + 15.93*s + 6.686e06)/((s^2*(s+352.5)*(1+s/2/pi/2000))); Kv = 10207*(s+55.15)*(s^2 + 17.45*s + 2.491e06)/(s^2*(s+491.2)*(s+7695)); K = blkdiag(Kh, Kv); K.InputName = {'Rh', 'Rv'}; K.OutputName = {'Uch', 'Ucv'};
We then close the loop and we look at the transfer function from the Newport rotation signal to the beam angle (Fig. 42).
Copyinputs = {'Uch', 'Ucv', 'Unh', 'Unv'}; outputs = {'Vch', 'Vcv', 'Ich', 'Icv', 'Rh', 'Rv', 'Vph', 'Vpv'}; sys_cl = connect(sys, -K, inputs, outputs);
9.3 Save the Controller
CopyKd = c2d(K, 1e-4, 'tustin');
The diagonal controller is accessible here.
Copysave('mat/K_diag.mat', 'K', 'Kd');
10 Newport Control
In this section, we try to implement a simple decentralized controller for the Newport. This can be used to align the 4QD:
- once there is a signal from the 4QD, the Newport feedback loop is closed
- thus, the Newport is positioned such that the beam hits the center of the 4QD
- then we can move the 4QD manually in X-Y plane in order to cancel the command signal of the Newport
- finally, we are sure to be aligned when the command signal of the Newport is 0
10.1 Load Plant
Copyload('mat/plant.mat', 'Gn', 'Gd');
10.2 Analysis
The plant is basically a constant until frequencies up to the required bandwidth.
We get that constant value.
CopyGn0 = freqresp(inv(Gd)*Gn, 0);
We design two controller containing 2 integrators and one lead near the crossover frequency set to 10Hz.
Copyh = 2; w0 = 2*pi*10; Knh = 1/Gn0(1,1) * (w0/s)^2 * (1 + s/w0*h)/(1 + s/w0/h)/h; Knv = 1/Gn0(2,2) * (w0/s)^2 * (1 + s/w0*h)/(1 + s/w0/h)/h;
10.3 Save
CopyKn = blkdiag(Knh, Knv); Knd = c2d(Kn, 1e-4, 'tustin');
The controllers can be downloaded here.
Copysave('mat/K_newport.mat', 'Kn', 'Knd');
11 Measurement of the non-repeatability
The goal here is the measure the non-repeatability of the setup.
All sources of error (detailed in the budget error in Section 4) will contribute to the non-repeatability of the system.
11.1 Data Load and pre-processing
Copyuh = load('mat/data_rep_h.mat', ... 't', 'Uch', 'Ucv', ... 'Unh', 'Unv', ... 'Vph', 'Vpv', ... 'Vnh', 'Vnv', ... 'Va'); uv = load('mat/data_rep_v.mat', ... 't', 'Uch', 'Ucv', ... 'Unh', 'Unv', ... 'Vph', 'Vpv', ... 'Vnh', 'Vnv', ... 'Va');
Copy% Let's start one second after the first command in the system i0 = find(uh.Unh ~= 0, 1) + fs; iend = i0+fs*floor((length(uh.t)-i0)/fs); uh.Uch([1:i0-1, iend:end]) = []; uh.Ucv([1:i0-1, iend:end]) = []; uh.Unh([1:i0-1, iend:end]) = []; uh.Unv([1:i0-1, iend:end]) = []; uh.Vph([1:i0-1, iend:end]) = []; uh.Vpv([1:i0-1, iend:end]) = []; uh.Vnh([1:i0-1, iend:end]) = []; uh.Vnv([1:i0-1, iend:end]) = []; uh.Va ([1:i0-1, iend:end]) = []; uh.t ([1:i0-1, iend:end]) = []; % We reset the time t uh.t = uh.t - uh.t(1);
Copy% Let's start one second after the first command in the system i0 = find(uv.Unv ~= 0, 1) + fs; iend = i0+fs*floor((length(uv.t)-i0)/fs); uv.Uch([1:i0-1, iend:end]) = []; uv.Ucv([1:i0-1, iend:end]) = []; uv.Unh([1:i0-1, iend:end]) = []; uv.Unv([1:i0-1, iend:end]) = []; uv.Vph([1:i0-1, iend:end]) = []; uv.Vpv([1:i0-1, iend:end]) = []; uv.Vnh([1:i0-1, iend:end]) = []; uv.Vnv([1:i0-1, iend:end]) = []; uv.Va ([1:i0-1, iend:end]) = []; uv.t ([1:i0-1, iend:end]) = []; % We reset the time t uv.t = uv.t - uv.t(1);
11.2 Some Time domain plots
11.3 Verify Tracking Angle Error
Let’s verify that the positioning error of the beam is small and what could be the effect on the distance measured by the intereferometer.
Copyload('./mat/plant.mat', 'Gd');
Let’s compute the PSD of the error to see the frequency content.
Copy[psd_UhRh, f] = pwelch(uh.Vph/freqresp(Gd(1,1), 0), hanning(ceil(1*fs)), [], [], fs); [psd_UhRv, ~] = pwelch(uh.Vpv/freqresp(Gd(2,2), 0), hanning(ceil(1*fs)), [], [], fs); [psd_UvRh, ~] = pwelch(uv.Vph/freqresp(Gd(1,1), 0), hanning(ceil(1*fs)), [], [], fs); [psd_UvRv, ~] = pwelch(uv.Vpv/freqresp(Gd(2,2), 0), hanning(ceil(1*fs)), [], [], fs);
Let’s convert that to errors in distance
is the nominal distance traveled by the beam is the distance traveled by the beam with an angle error is the angle error
CopyL = 0.1; % [m]
Copy[psd_UhLh, f] = pwelch(0.5*L*(uh.Vph/freqresp(Gd(1,1), 0)).^2, hanning(ceil(1*fs)), [], [], fs); [psd_UhLv, ~] = pwelch(0.5*L*(uh.Vpv/freqresp(Gd(2,2), 0)).^2, hanning(ceil(1*fs)), [], [], fs); [psd_UvLh, ~] = pwelch(0.5*L*(uv.Vph/freqresp(Gd(1,1), 0)).^2, hanning(ceil(1*fs)), [], [], fs); [psd_UvLv, ~] = pwelch(0.5*L*(uv.Vpv/freqresp(Gd(2,2), 0)).^2, hanning(ceil(1*fs)), [], [], fs);
Now, compare that with the PSD of the measured distance by the interferometer (Fig. 47).
Copy[psd_Lh, f] = pwelch(uh.Va, hanning(ceil(1*fs)), [], [], fs); [psd_Lv, ~] = pwelch(uv.Va, hanning(ceil(1*fs)), [], [], fs);
Figure 47: Comparison of the effect of tracking error on the measured distance and the measured distance by the Attocube (png, pdf)
Important
The tracking errors are a limiting factor.
11.4 Processing
First, we get the mean value as measured by the interferometer for each value of the Newport angle.
CopyVahm = mean(reshape(uh.Va, [fs floor(length(uh.t)/fs)]),2); Unhm = mean(reshape(uh.Unh, [fs floor(length(uh.t)/fs)]),2); Vavm = mean(reshape(uv.Va, [fs floor(length(uv.t)/fs)]),2); Unvm = mean(reshape(uv.Unv, [fs floor(length(uv.t)/fs)]),2);
And we can compute the RMS value of the non-repeatable part:
Va - Horizontal [nm rms] | Va - Vertical [nm rms] |
---|---|
19.6 | 13.9 |
11.5 Analysis of the non-repeatable contributions
Let’s know try to determine where does the non-repeatability comes from.
From the 4QD signal, we can compute the angle error of the beam and thus determine the corresponding displacement measured by the attocube.
We then take the non-repeatable part of this displacement and we compare that with the total non-repeatability.
We also plot the displacement measured during the huddle test.
All the signals are shown on Fig. 50.
11.6 Results with a low pass filter
We filter the data with a first order low pass filter with a crossover frequency of
Copyw0 = 10; % [Hz] G_lpf = 1/(1 + s/2/pi/w0); uh.Vaf = lsim(G_lpf, uh.Va, uh.t); uv.Vaf = lsim(G_lpf, uv.Va, uv.t);
11.7 Processing
First, we get the mean value as measured by the interferometer for each value of the Newport angle.
CopyVahm = mean(reshape(uh.Vaf, [fs floor(length(uh.t)/fs)]),2); Unhm = mean(reshape(uh.Unh, [fs floor(length(uh.t)/fs)]),2); Vavm = mean(reshape(uv.Vaf, [fs floor(length(uv.t)/fs)]),2); Unvm = mean(reshape(uv.Unv, [fs floor(length(uv.t)/fs)]),2);
And we can compute the RMS value of the non-repeatable part:
Va - Horizontal [nm rms] | Va - Vertical [nm rms] |
---|---|
22.9 | 13.9 |