-
-
- 1. Estimation of Sensitivity Function +
- 1. Estimation of Sensitivity Function -
- 2. System Identification +
- 2. System Identification
-
-
- 2.1. Identification -
- 2.2. Identification - New -
- 2.3. Identification - White noise -
- 2.4. Test -
- 2.5. New controller - Higher bandwidth -
- 2.6. Added gain +
- 2.1. Identification +
- 2.2. Identification - New +
- 2.3. Identification - White noise +
- 2.4. Test +
- 2.5. New controller - Higher bandwidth +
- 2.6. Added gain
- - 3. Noise Budgeting +
- 3. Noise Budgeting -
- 4. Test Mode C +
- 4. Test Mode C -
- 5. Export numerator and denominator +
- 5. Export numerator and denominator
This report is also available as a pdf.
-
-
1. Estimation of Sensitivity Function
+
+
-
1. Estimation of Sensitivity Function
-
1.1. Load Data
+
+
-
1.1. Load Data
Two scans are performed: @@ -99,7 +99,7 @@ t = linspace(Ts, Ts*length(ol_drx), length(ol_d
-By comparison the frequency content of the crystal orientation errors between mode B and mode C, it is possible to estimate the Sensitivity transfer function (Figure 1). +By comparison the frequency content of the crystal orientation errors between mode B and mode C, it is possible to estimate the Sensitivity transfer function (Figure 1).
win = hanning(ceil(1/Ts));
@@ -113,7 +113,7 @@ By comparison the frequency content of the crystal orientation errors between mo
+
-
Figure 1: Estimation of the sensitivity transfer function magnitude
@@ -122,8 +122,8 @@ By comparison the frequency content of the crystal orientation errors between mo
-
1.2. Controller
+
+
-
1.2. Controller
load('X_tal_cage_PID.mat', 'K'); @@ -131,7 +131,7 @@ By comparison the frequency content of the crystal orientation errors between mo
+
-
Figure 2: Bode Plot of the Controller
@@ -139,8 +139,8 @@ By comparison the frequency content of the crystal orientation errors between mo
-
1.3. Test
+
+
-
1.3. Test
Ts = 5e-3; @@ -170,12 +170,12 @@ t = linspace(Ts, Ts*length(ol_drx), length(ol_d
-
2. System Identification
+
+
-
2. System Identification
-
2.1. Identification
+
+
-
2.1. Identification
ur = load('FJPUR_step.mat'); @@ -321,7 +321,7 @@ win = hanning(ceil(5*Fs));
+
Figure 3: Coherence
@@ -329,7 +329,7 @@ win = hanning(ceil(5*Fs)); -
+
Figure 4: Bode Plot of the DCM dynamics in the frame of the fast jack.
@@ -359,7 +359,7 @@ Kb = eye(3)*(2* -
+
Figure 5: Loop gain
@@ -367,7 +367,7 @@ Kb = eye(3)*(2* +
Figure 6: Loop gain
@@ -377,7 +377,7 @@ Kb = eye(3)*(2* -
+
-
Figure 7: Comparison of sensitivity functions
@@ -396,7 +396,7 @@ Lb = zeros(3, 3, length(f));
+
-
Figure 8: Root Locus
@@ -405,8 +405,8 @@ Lb = zeros(3, 3, length(f));
-
2.2. Identification - New
+
+
2.2. Identification - New
ur = load('FJPUR_step_new.mat'); @@ -577,7 +577,7 @@ Kb = eye(3)*(2* -+-
Figure 9: Comparison of sensitivity functions
@@ -596,7 +596,7 @@ Lb = zeros(3, 3, length(f));+-
Figure 10: Root Locus
@@ -605,8 +605,8 @@ Lb = zeros(3, 3, length(f));-2.3. Identification - White noise
++2.3. Identification - White noise
-ur = load('fjpur_white_noise.mat'); @@ -725,14 +725,14 @@ win = hanning(ceil(0.5*Fs));+-
Figure 11: description
+
Figure 12: Bode Plot of the DCM dynamics in the frame of the fast jack.
@@ -762,14 +762,14 @@ Kb = eye(3)*(2* -+-
Figure 13: Loop gain
+
Figure 14: Loop gain
@@ -779,7 +779,7 @@ Kb = eye(3)*(2* -+-
Figure 15: Comparison of sensitivity functions
@@ -798,7 +798,7 @@ Lb = zeros(3, 3, length(f));+-
Figure 16: Root Locus
@@ -807,8 +807,8 @@ Lb = zeros(3, 3, length(f));-2.4. Test
++2.4. Test
-%% Notch @@ -851,8 +851,8 @@ Kb = 0.8*eye(3)*-2.5. New controller - Higher bandwidth
++2.5. New controller - Higher bandwidth
%% Previously used controller @@ -916,7 +916,7 @@ Kb = 0.9*eye(3)* -+-
Figure 17: description
@@ -937,14 +937,14 @@ Lb_new = zeros(3, 3, length(f));+-
Figure 18: n
+
Figure 19: description
@@ -954,8 +954,8 @@ Lb_new = zeros(3, 3, length(f)); --2.6. Added gain
++2.6. Added gain
-%% Notch @@ -1017,21 +1017,21 @@ Kb_gain = 0.9*eye(3)+-
Figure 20: description
+-
Figure 21: description
+-
Figure 22: nyquist plot
@@ -1041,12 +1041,12 @@ Kb_gain = 0.9*eye(3)-3. Noise Budgeting
++3. Noise Budgeting
--3.1. No Displacement
++3.1. No Displacement
@@ -1275,8 +1275,8 @@ data_70_deg.drx = data_70_deg.allValues(:,6)
-3.2. Scans
++3.2. Scans
@@ -1437,8 +1437,8 @@ CPS_70_10_d = flip(-cumtrapz(flip(f), flip(pxx -
-3.3. Noise budgeting - No rotation
++3.3. Noise budgeting - No rotation
-First, we look at the position errors when the bragg axis is not moving @@ -1488,7 +1488,7 @@ win = hanning(ceil(1/Ts));
+-
Figure 23: Amplitude Spectral Density
@@ -1511,7 +1511,7 @@ CPS_dz = cumtrapz(f, pxx_ol_dz);+-
Figure 24: Cumulative Amplitude Spectrum
@@ -1520,17 +1520,17 @@ CPS_dz = cumtrapz(f, pxx_ol_dz);--3.4. Noise budgeting - Bragg rotation
++3.4. Noise budgeting - Bragg rotation
-4. Test Mode C
++4. Test Mode C
--4.1. Mode B and Mode C
++4.1. Mode B and Mode C
-data_B = extractDatData(sprintf("/home/thomas/mnt/data_id21/22Jan/blc13491/id21/test_regul_220119/%s","lut_const_fj_vel_19012022_1450.dat"), ... @@ -1620,12 +1620,12 @@ win = hanning(ceil(1/Ts));-5. Export numerator and denominator
++5. Export numerator and denominator
--5.1. Export
++5.1. Export
-K_order = 10; @@ -1673,8 +1673,8 @@ fclose(fileID);-5.2. Verify
++5.2. Verify
diff --git a/dcm-feedback-control.org b/dcm-feedback-control.org index 629e5f6..a5f3a09 100644 --- a/dcm-feedback-control.org +++ b/dcm-feedback-control.org @@ -235,6 +235,125 @@ legend('location', 'northwest'); <K_data = importdata('X_tal_cage_PID_20Hz.dat'); @@ -1687,7 +1687,7 @@ K = tf(K_data(1,:), K_data(2, -Created: 2022-02-15 mar. 14:18
+Created: 2022-02-20 dim. 23:36
> #+end_src +** Identification ID24 +#+begin_src matlab +load('test_id_id24_3.mat') +#+end_src + +#+begin_src matlab +t = 1e-4*ones(size(fjpur, 1), 1); + +ur.dz = fjpur(:,1) - mean(fjpur(:,1)); +ur.dry = fjpur(:,2) - mean(fjpur(:,2)); +ur.drx = fjpur(:,3) - mean(fjpur(:,3)); +ur.u = fjpur(:,7) - mean(fjpur(:,7)); + +uh.dz = fjpuh(:,1) - mean(fjpuh(:,1)); +uh.dry = fjpuh(:,2) - mean(fjpuh(:,2)); +uh.drx = fjpuh(:,3) - mean(fjpuh(:,3)); +uh.u = fjpuh(:,8) - mean(fjpuh(:,8)); + +d.dz = fjpd(:,1) - mean(fjpd(:,1)); +d.dry = fjpd(:,2) - mean(fjpd(:,2)); +d.drx = fjpd(:,3) - mean(fjpd(:,3)); +d.u = fjpd(:,9) - mean(fjpd(:,9)); +#+end_src + +#+begin_src matlab +J_a_311 = [1, 0.14, -0.0675 + 1, 0.14, 0.1525 + 1, -0.14, 0.0425]; + +J_a_111 = [1, 0.14, -0.1525 + 1, 0.14, 0.0675 + 1, -0.14, -0.0425]; + +ur.y = [J_a_311 * [-ur.dz, ur.dry,-ur.drx]']'; +uh.y = [J_a_311 * [-uh.dz, uh.dry,-uh.drx]']'; +d.y = [J_a_311 * [-d.dz, d.dry, -d.drx]']'; +#+end_src + +#+begin_src matlab +%% Sampling Time and Frequency +Ts = 1e-4; % [s] +Fs = 1/Ts; % [Hz] + +% Hannning Windows +win = hanning(ceil(1*Fs)); +#+end_src + +#+begin_src matlab +%% And we get the frequency vector +[G_ur, f] = tfestimate(ur.u, ur.y, win, [], [], 1/Ts); +[G_uh, ~] = tfestimate(uh.u, uh.y, win, [], [], 1/Ts); +[G_d, ~] = tfestimate(d.u, d.y, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab +[coh_ur, ~] = mscohere(ur.u, ur.y, win, [], [], 1/Ts); +[coh_uh, ~] = mscohere(uh.u, uh.y, win, [], [], 1/Ts); +[coh_d, ~] = mscohere(d.u, d.y, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab :exports none +%% Coherence +figure; +tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); +hold on; +plot(f, coh_ur(:,1), 'DisplayName', '$u_r$'); +plot(f, coh_uh(:,2), 'DisplayName', '$u_h$'); +plot(f, coh_d( :,3), 'DisplayName', '$u_d$'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +xlabel('Frequency [Hz]'); ylabel('Coherence'); +legend('location', 'southeast'); +xlim([1, 1e3]); +#+end_src + +#+begin_src matlab :exports none +%% Bode plot of the DCM dynamics in the frame of the Fast Jacks +figure; +tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); + +ax1 = nexttile([2,1]); +hold on; +plot(f, abs(G_ur(:,1)), 'DisplayName', '$u_r$'); +plot(f, abs(G_uh(:,2)), 'DisplayName', '$u_h$'); +plot(f, abs(G_d( :,3)), 'DisplayName', '$u_d$'); +plot(f, abs(G_ur(:,2)), 'color', [0, 0, 0, 0.2], ... + 'DisplayName', 'Off Diagonal'); +plot(f, abs(G_ur(:,3)), 'color', [0, 0, 0, 0.2], ... + 'HandleVisibility', 'off'); +plot(f, abs(G_uh(:,1)), 'color', [0, 0, 0, 0.2], ... + 'HandleVisibility', 'off'); +plot(f, abs(G_uh(:,3)), 'color', [0, 0, 0, 0.2], ... + 'HandleVisibility', 'off'); +plot(f, abs(G_d(:,1)), 'color', [0, 0, 0, 0.2], ... + 'HandleVisibility', 'off'); +plot(f, abs(G_d(:,2)), 'color', [0, 0, 0, 0.2], ... + 'HandleVisibility', 'off'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Amplitude'); set(gca, 'XTickLabel',[]); +ylim([1e-2, 1e2]); +legend('location', 'northwest'); + +ax2 = nexttile; +hold on; +plot(f, 180/pi*angle(G_ur(:,1))); +plot(f, 180/pi*angle(G_uh(:,2))); +plot(f, 180/pi*angle(G_d(:,3))); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +hold off; +yticks(-360:90:360); +% ylim([-45, 45]); + +linkaxes([ax1,ax2],'x'); +xlim([1, 1e3]); +#+end_src + ** Identification #+begin_src matlab ur = load('FJPUR_step.mat'); @@ -1364,7 +1483,7 @@ exportFig('figs/nyquist_compare.pdf', 'width', 'wide', 'height', 'tall'); #+end_src #+name: fig:nyquist_compare -#+caption: n +#+caption: Nyquist Plot #+RESULTS: [[file:figs/nyquist_compare.png]] @@ -1375,12 +1494,12 @@ tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); hold on; plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(K(1,1), f, 'Hz'))), 'DisplayName', 'PID'); -plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(Kb_old(1,1), f, 'Hz'))), 'DisplayName', 'Double Integrator + Lead'); -plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'DisplayName', 'Double Integrator + Lead + Notch'); +plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(Kb_old(1,1), f, 'Hz'))), 'DisplayName', '$K_{10Hz}$'); +plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'DisplayName', '$K_{20Hz}$'); hold off; set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -xlabel('Frequency [Hz]'); ylabel('Amplitude'); -ylim([1e-2, 1e1]); xlim([1, 1e3]); +xlabel('Frequency [Hz]'); ylabel('Sensitivity Magnitude'); +ylim([1e-3, 5]); xlim([1, 5e2]); legend('location', 'southeast'); #+end_src @@ -2257,7 +2376,7 @@ Fs = 1e4; % Sampling Frequency [Hz] fir_order = args.fir_order; % Filter's order delay = fir_order/2; % Delay induced by the filter B_fir = firls(args.fir_order, ... % Filter's order - [0 500/(Fs/2) 1e3/(Fs/2) 1], ... % Frequencies [Hz] + [0 25/(Fs/2) 34/(Fs/2) 1], ... % Frequencies [Hz] [1 1 0 0]); % Wanted Magnitudes #+end_src diff --git a/dcm-feedback-control.pdf b/dcm-feedback-control.pdf new file mode 100644 index 0000000..a1197e5 Binary files /dev/null and b/dcm-feedback-control.pdf differ diff --git a/figs/sensitivity_function_compare.pdf b/figs/sensitivity_function_compare.pdf index 44a6823..0a1f1e2 100644 Binary files a/figs/sensitivity_function_compare.pdf and b/figs/sensitivity_function_compare.pdf differ diff --git a/figs/sensitivity_function_compare.png b/figs/sensitivity_function_compare.png index 47cb5ab..58c1690 100644 Binary files a/figs/sensitivity_function_compare.png and b/figs/sensitivity_function_compare.png differ diff --git a/matlab/DCM_RT_control_v09.slx b/matlab/DCM_RT_control_v09.slx new file mode 100755 index 0000000..e4196d1 Binary files /dev/null and b/matlab/DCM_RT_control_v09.slx differ diff --git a/matlab/DCM_RT_control_v09.slx.r2019b b/matlab/DCM_RT_control_v09.slx.r2019b new file mode 100755 index 0000000..24da707 Binary files /dev/null and b/matlab/DCM_RT_control_v09.slx.r2019b differ diff --git a/matlab/calc_kinematics.m b/matlab/calc_kinematics.m new file mode 100755 index 0000000..a62574a --- /dev/null +++ b/matlab/calc_kinematics.m @@ -0,0 +1,32 @@ +%% this script computes the transforamtion matrices for piezo to pseudo coordinates + +fjpur_111 = [-140; -67.5; 0]*0.001; +fjpuh_111 = [-140; (220-67.5); 0]*0.001; +fjpd_111 = [140; (110-67.5); 0]*0.001; + +fjpur_311 = [-140; -(220-67.5); 0]*0.001; +fjpuh_311 = [-140; 67.5; 0]*0.001; +fjpd_311 = [140; -(110-67.5); 0]*0.001; + +A_p2r_111 = [1 -1*fjpur_111(1) fjpur_111(2); 1 -1*fjpuh_111(1) fjpuh_111(2); 1 -1*fjpd_111(1) fjpd_111(2)]; +A_p2r_311 = [1 -1*fjpur_311(1) fjpur_311(2); 1 -1*fjpuh_311(1) fjpuh_311(2); 1 -1*fjpd_311(1) fjpd_311(2)]; + +A_r2p_111 = inv(A_p2r_111); +A_r2p_311 = inv(A_p2r_311); + +A_x2p = [-1 0 0; 0 1 0; 0 0 1]; +B_x2p = [0 0 0]; + +motor_pos2mv = 0.6667; %modified by Ludovic 12/07/2021 +% % no calibr +calibr_fjpur = (1); +calibr_fjpuh = (1); +calibr_fjpd = (1); + +%calibr +% calibr_fjpur = (1-0.016); +% calibr_fjpuh = (1+0.052); +% calibr_fjpd = (1+0.0065); + +calibr_mat = [calibr_fjpur 0 0;0 calibr_fjpuh 0;0 0 calibr_fjpd]; + diff --git a/matlab/figs b/matlab/figs new file mode 120000 index 0000000..8ea5186 --- /dev/null +++ b/matlab/figs @@ -0,0 +1 @@ +../figs \ No newline at end of file diff --git a/matlab/src/extractDatData.m b/matlab/src/extractDatData.m new file mode 100644 index 0000000..8f83886 --- /dev/null +++ b/matlab/src/extractDatData.m @@ -0,0 +1,22 @@ +function [data_struct] = extractDatData(dat_file, names, scale) +% extractDatData - +% +% Syntax: extractDatData(data_file, lut_file, args) +% +% Inputs: +% - data_file - Where to load the .mat file +% - lut_file - Where to save the .dat file + +%% Load Data +data_array = importdata(dat_file); + +%% Initialize Struct +data_struct = struct(); + +%% Populate Struct +for i = 1:length(names) + data_struct.(names{i}) = scale(i)*data_array(:,i); +end + +%% Add Time +data_struct.time = 1e-4*[1:1:length(data_struct.(names{1}))]; diff --git a/matlab/src/processMeasData.m b/matlab/src/processMeasData.m new file mode 100644 index 0000000..3b4ac7d --- /dev/null +++ b/matlab/src/processMeasData.m @@ -0,0 +1,49 @@ +function [data] = processMeasData(data, args) +% processMeasData - +% +% Syntax: processMeasData(data_file, lut_file, args) +% +% Inputs: +% - data + +arguments + data + args.fir_order (1,1) double {mustBeNumericOrLogical} = 5000 +end + +%% Actuator Jacobian +J_a_111 = [1, 0.14, -0.0675 + 1, 0.14, 0.1525 + 1, -0.14, 0.0425]; + +%% FIR Filter +Fs = 1e4; % Sampling Frequency [Hz] +fir_order = args.fir_order; % Filter's order +delay = fir_order/2; % Delay induced by the filter +B_fir = firls(args.fir_order, ... % Filter's order + [0 50/(Fs/2) 0.1e3/(Fs/2) 1], ... % Frequencies [Hz] + [1 1 0 0]); % Wanted Magnitudes + +data.ddz = 10.5e-3./(2*cos(data.bragg)) - data.dz; + +%% Computation of the position of the FJ as measured by the interferometers +error = J_a_111 * [data.ddz, data.dry, data.drx]'; + +data.fjur_e = error(1,:)'; % [m] +data.fjuh_e = error(2,:)'; % [m] +data.fjd_e = error(3,:)'; % [m] + +%% Filtering all measured Fast Jack Position using the FIR filter +data.fjur_e_filt = filter(B_fir, 1, data.fjur_e); +data.fjuh_e_filt = filter(B_fir, 1, data.fjuh_e); +data.fjd_e_filt = filter(B_fir, 1, data.fjd_e); + +%% Compensation of the delay introduced by the FIR filter +data.fjur_e_filt(1:end-delay) = data.fjur_e_filt(delay+1:end); +data.fjuh_e_filt(1:end-delay) = data.fjuh_e_filt(delay+1:end); +data.fjd_e_filt( 1:end-delay) = data.fjd_e_filt( delay+1:end); + +%% Re-sample data to have same data points in FJUR +[data.fjur_e_resampl, data.fjur_resampl] = resample(data.fjur_e_filt, data.fjur, 1/100e-9); +[data.fjuh_e_resampl, data.fjuh_resampl] = resample(data.fjuh_e_filt, data.fjuh, 1/100e-9); +[data.fjd_e_resampl, data.fjd_resampl] = resample(data.fjd_e_filt, data.fjd, 1/100e-9);