From 151ae070727ec852a46ff69ca6cdfce6fafd35fd Mon Sep 17 00:00:00 2001 From: Thomas Dehaeze Date: Wed, 30 Oct 2024 19:21:34 +0100 Subject: [PATCH] Add compliance measurement --- simscape-micro-station.org | 586 ++++++++++++++++++++++++++++++++++++- 1 file changed, 582 insertions(+), 4 deletions(-) diff --git a/simscape-micro-station.org b/simscape-micro-station.org index ccc0db2..1601ed1 100644 --- a/simscape-micro-station.org +++ b/simscape-micro-station.org @@ -164,9 +164,8 @@ initializeController( 'type', 'open-loop'); #+end_src -** TODO [#B] Make good "init" for the Simscape model - -** TODO [#B] Just keep smallest number of variants for each stage +** DONE [#B] Just keep smallest number of variants for each stage +CLOSED: [2024-10-30 Wed 16:15] - [ ] none - [ ] rigid @@ -178,10 +177,41 @@ initializeController( 'type', 'open-loop'); SCHEDULED: <2024-10-30 Wed> - [ ] Find the compliance measurements + - The one from modal analysis + - [[file:~/Cloud/work-projects/ID31-NASS/matlab/nass-measurements/2018-01-12 - Marc]] : Compliance measurement in X,Y,Z with geophone, marble not glued + - [[file:~/Cloud/work-projects/ID31-NASS/matlab/nass-measurements/2018-10-12 - Marc]]: same but in the hutch with glued marble + - *08/2020*: [[file:~/Cloud/work-projects/ID31-NASS/matlab/nass-measurements/micro-station-compliance/index.org::+TITLE: Compliance Measurement of the Micro Station]] - [ ] See if it matches somehow the current model - [ ] If not, see if model parameters can be tuned to have better match For instance from values here: file:/home/thomas/Cloud/meetings/esrf-meetings/2018-04-24-Simscape-Model/2018-04-24-Simscape-Model.pdf +** TODO [#B] Make good "init" for the Simscape model + +** TODO [#C] Could see the effect of each stage on the compliance + +- [ ] Put =rigid= mode one by one to see the effect + +** TODO [#C] Make a comparison of the measured vibrations of the micro-station with the vibrations of the simscape model of the micro-station +Do we have a correlation? At least in the frequency domain? + +Procedure: +- [ ] Take the time domain measurement of the vibrations due to spindle or translation stage +- [ ] Take the estimated PSD of the estimated disturbance force +- [ ] Simulate the Simscape model without the nano-hexapod (same conditions as when measuring the disturbances) and add only the considered disturbance +- [ ] Save the position errors due to the disturbance +- [ ] Compare with the measured vibrations + + +** TODO [#C] Add two perturbations: =Frz_x= and =Frz_y= +Maybe I can estimate them from the measurements that was made on the spindle? + +** TODO [#C] Determine which Degree-Of-Freedom to keep and which to constrain + +For instance, for now the granite can not rotate, but in reality, the modes may be linked to the granite's rotation. +By constraining more DoF, the simulation will be faster and the obtain state space will have a lower order. + +** TODO [#B] Compute eigenvalues of the model to see if we have similar frequencies than the modal analysis? + * Introduction :ignore: @@ -894,7 +924,385 @@ end end #+end_src -** Obtained Compliance of the Micro-Station +** Measured micro-station compliance +*** Introduction :ignore: + +The most important dynamical characteristic of the micro-station that should be well modeled is its compliance. +This is what can impact the nano-hexapod dynamics. + +- [ ] Add schematic of the experiment with Micro-Hexapod top platform, location of accelerometers, of impacts, etc... + +- 4 3-axis accelerometers +- 10 hammer impacts on the micro-hexapod top plaftorm +- *Was the rotation compensation axis present?* (I don't think so) + +*** Position of inertial sensors on top of the micro-hexapod +Orientation is relative to the frame determined by the X-ray + +| *Num* | *Position* | *Orientation* | *Sensibility* | *Channels* | +|-------+------------+---------------+---------------+------------| +| 1 | [0, +A, 0] | [x, y, z] | 1V/g | 1-3 | +| 2 | [-B, 0, 0] | [x, y, z] | 1V/g | 4-6 | +| 3 | [0, -A, 0] | [x, y, z] | 0.1V/g | 7-9 | +| 4 | [+B, 0, 0] | [x, y, z] | 1V/g | 10-12 | + +Instrumented Hammer: +- Channel 13 +- Sensibility: 230 uV/N + +| Acc Number | Dir | Channel Number | +|------------+-----+----------------| +| 1 | x | 1 | +| 1 | y | 2 | +| 1 | z | 3 | +| 2 | x | 4 | +| 2 | y | 5 | +| 2 | z | 6 | +| 3 | x | 7 | +| 3 | y | 8 | +| 3 | z | 9 | +| 4 | x | 10 | +| 4 | y | 11 | +| 4 | z | 12 | +| Hammer | | 13 | + +From the acceleration measurement of the 4 accelerometers, we can compute the translations and rotations: +| | *Formula* | *Formula (numbers)* | +|-------+--------------------------+-----------------------| +| $D_x$ | (1x + 2x + 3x + 4x)/4 | (1 + 4 + 7 + 10)/4 | +| $D_y$ | (1y + 2y + 3y + 4y)/4 | (2 + 5 + 8 + 11)/4 | +| $D_z$ | (1z + 2z + 3z + 4z)/4 | (3 + 6 + 9 + 12)/4 | +| $R_x$ | (1z - 3z)/A | (1 - 9)/A | +| $R_y$ | (2z - 4z)/B | (6 - 12)/B | +| $R_z$ | (3x - 1x)/A, (4y - 2y)/B | (7 - 1)/A, (11 - 5)/B | + +dL = J X +#+begin_src matlab +d = 0.14; +J = [1 0 0 0 0 -d; + 0 1 0 0 0 0; + 0 0 1 d 0 0; + 1 0 0 0 0 0; + 0 1 0 0 0 -d; + 0 0 1 0 d 0; + 1 0 0 0 0 d; + 0 1 0 0 0 0; + 0 0 1 -d 0 0; + 1 0 0 0 0 0; + 0 1 0 0 0 d; + 0 0 1 0 -d 0]; + +J_inv = pinv(J); +#+end_src + +| | Dx | Dy | Dz | Rx | Ry | Rz | +|-----+----+----+----+----+----+----| +| a1x | 1 | 0 | 0 | 0 | 0 | -d | +| a1y | 0 | 1 | 0 | 0 | 0 | 0 | +| a1z | 0 | 0 | 1 | d | 0 | 0 | +| a2x | 1 | 0 | 0 | 0 | 0 | 0 | +| a2y | 0 | 1 | 0 | 0 | 0 | -d | +| a2z | 0 | 0 | 1 | 0 | d | 0 | +| a3x | 1 | 0 | 0 | 0 | 0 | d | +| a3y | 0 | 1 | 0 | 0 | 0 | 0 | +| a3z | 0 | 0 | 1 | -d | 0 | 0 | +| a4x | 1 | 0 | 0 | 0 | 0 | 0 | +| a4y | 0 | 1 | 0 | 0 | 0 | d | +| a4z | 0 | 0 | 1 | 0 | -d | 0 | + +#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*) +data2orgtable(J_inv, {'Dx', 'Dy', 'Dz', 'Rx', 'Ry', 'Rz'}, {'d1x', 'd1y', 'd1z', 'd2x', 'd2y', 'd2z', 'd3x', 'd3y', 'd3z', 'd4x', 'd4y', 'd4z'}, ' %.5f '); +#+end_src + +#+RESULTS: +| | d1x | d1y | d1z | d2x | d2y | d2z | d3x | d3y | d3z | d4x | d4y | d4z | +|----+----------+------+---------+------+----------+---------+---------+------+----------+------+---------+----------| +| Dx | 0.25 | 0.0 | 0.0 | 0.25 | 0.0 | 0.0 | 0.25 | 0.0 | 0.0 | 0.25 | 0.0 | 0.0 | +| Dy | 0.0 | 0.25 | 0.0 | 0.0 | 0.25 | 0.0 | 0.0 | 0.25 | 0.0 | 0.0 | 0.25 | 0.0 | +| Dz | 0.0 | 0.0 | 0.25 | 0.0 | 0.0 | 0.25 | 0.0 | 0.0 | 0.25 | 0.0 | 0.0 | 0.25 | +| Rx | 0.0 | 0.0 | 3.57143 | 0.0 | 0.0 | -0.0 | 0.0 | 0.0 | -3.57143 | 0.0 | 0.0 | -0.0 | +| Ry | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 3.57143 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | -3.57143 | +| Rz | -1.78571 | 0.0 | -0.0 | 0.0 | -1.78571 | 0.0 | 1.78571 | 0.0 | -0.0 | 0.0 | 1.78571 | -0.0 | + +*** Hammer blow position/orientation + +| *Num* | *Direction* | *Position* | Accelerometer position | Jacobian number | +|-------+-------------+------------+------------------------+-----------------| +| 1 | -Y | [0, +A, 0] | 1 | -2 | +| 2 | -Z | [0, +A, 0] | 1 | -3 | +| 3 | X | [-B, 0, 0] | 2 | 4 | +| 4 | -Z | [-B, 0, 0] | 2 | -6 | +| 5 | Y | [0, -A, 0] | 3 | 8 | +| 6 | -Z | [0, -A, 0] | 3 | -9 | +| 7 | -X | [+B, 0, 0] | 4 | -10 | +| 8 | -Z | [+B, 0, 0] | 4 | -12 | +| 9 | -X | [0, -A, 0] | 3 | -7 | +| 10 | -X | [0, +A, 0] | 1 | -1 | + +From hammer blows to pure forces / torques: +| | *Formula* | Alternative | +|-------+--------------+-------------| +| $F_x$ | +3 | -7 | +| $F_y$ | -1 | +5 | +| $F_z$ | -(2 + 6)/2 | -(4 + 8)/2 | +| $M_x$ | A/2*(2 - 6) | | +| $M_y$ | B/2*(8 - 4) | | +| $M_z$ | A/2*(10 - 9) | | + +F = J' tau +#+begin_src matlab +Jf = [0 -1 0 0 0 0; + 0 0 -1 -d 0 0; + 1 0 0 0 0 0; + 0 0 -1 0 -d 0; + 0 1 0 0 0 0; + 0 0 -1 d 0 0; + -1 0 0 0 0 0; + 0 0 -1 0 d 0; + -1 0 0 0 0 -d; + -1 0 0 0 0 d]'; +Jf_inv = pinv(Jf); +#+end_src + +#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*) +data2orgtable(Jf, {'Fx', 'Fy', 'Fz', 'Mx', 'My', 'Mz'}, {'-F1y', '-F1z', '+F2x', '-F2z', '+F3y', '-F3z', '-F4x', '-F4z', '-F3x', '-F1x'}, ' %.2f '); +#+end_src + +#+RESULTS: +| | -F1y | -F1z | +F2x | -F2z | +F3y | -F3z | -F4x | -F4z | -F3x | -F1x | +|----+------+-------+------+-------+------+------+------+------+-------+------| +| Fx | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | -1.0 | 0.0 | -1.0 | -1.0 | +| Fy | -1.0 | 0.0 | 0.0 | 0.0 | 1.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | +| Fz | 0.0 | -1.0 | 0.0 | -1.0 | 0.0 | -1.0 | 0.0 | -1.0 | 0.0 | 0.0 | +| Mx | 0.0 | -0.14 | 0.0 | 0.0 | 0.0 | 0.14 | 0.0 | 0.0 | 0.0 | 0.0 | +| My | 0.0 | 0.0 | 0.0 | -0.14 | 0.0 | 0.0 | 0.0 | 0.14 | 0.0 | 0.0 | +| Mz | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | -0.14 | 0.14 | + +*** Compute FRF +Raw measurements are in file:/home/thomas/Cloud/work-projects/ID31-NASS/matlab/nass-measurements/micro-station-compliance/data/record + +For each measurement (10): +- Find the location of the 10 impacts based on "track13" +- Then for the 12 accelerometer data, compute the FRF, and average them for the 10 impacts +- Maybe have to take into account the sensitivity, etc... + +#+begin_src matlab +raw_data_path = '/home/thomas/Cloud/work-projects/ID31-NASS/matlab/nass-measurements/micro-station-compliance/data/record/'; +data = [ + load(sprintf('%s/Measurement1.mat', raw_data_path)), ... + load(sprintf('%s/Measurement2.mat', raw_data_path)), ... + load(sprintf('%s/Measurement3.mat', raw_data_path)), ... + load(sprintf('%s/Measurement4.mat', raw_data_path)), ... + load(sprintf('%s/Measurement5.mat', raw_data_path)), ... + load(sprintf('%s/Measurement6.mat', raw_data_path)), ... + load(sprintf('%s/Measurement7.mat', raw_data_path)), ... + load(sprintf('%s/Measurement8.mat', raw_data_path)), ... + load(sprintf('%s/Measurement9.mat', raw_data_path)), ... + load(sprintf('%s/Measurement10.mat', raw_data_path))]; +#+end_src + +#+begin_src matlab +Ts = 1e-3; % Sampling Time [s] +Nfft = floor(1/Ts); % Number of points for the FFT computation +% win = hanning(Nfft); % Hanning window +win = ones(Nfft, 1); % Rectangular window +% Get the frequency vector +[~, f] = tfestimate(data(1).Track13, data(1).Track1, win, [], Nfft, 1/Ts); +#+end_src + +#+begin_src matlab +pre_n = floor(0.1/Ts); +post_n = Nfft - pre_n - 1; +#+end_src + +#+begin_src matlab +G_raw = zeros(12,10,length(f)); +for i = 1:10 + % Find the impacts + [~, impacts_i] = find(diff(data(i).Track13 > 50)==1); + % Only keep the first 10 impacts if there are more than 10 impacts + if length(impacts_i)>10 + impacts_i(11:end) = []; + end + % Initialize the FRF for the current experiment + G_impact = zeros(12,length(f)); + for impact_i = impacts_i + i_start = impacts_i - pre_n; + i_end = impacts_i + post_n; + data_in = data(i).Track13(i_start:i_end); % [N] + % Remove hammer DC offset + data_in = data_in - mean(data_in(end-pre_n:end)); + % Combine all outputs [m/s^2] + data_out = [data(i).Track1( i_start:i_end); ... + data(i).Track2( i_start:i_end); ... + data(i).Track3( i_start:i_end); ... + data(i).Track4( i_start:i_end); ... + data(i).Track5( i_start:i_end); ... + data(i).Track6( i_start:i_end); ... + data(i).Track7( i_start:i_end); ... + data(i).Track8( i_start:i_end); ... + data(i).Track9( i_start:i_end); ... + data(i).Track10(i_start:i_end); ... + data(i).Track11(i_start:i_end); ... + data(i).Track12(i_start:i_end)]; + + [frf, ~] = tfestimate(data_in, data_out', win, [], Nfft, 1/Ts); + G_raw(:,i,:) = frf'./(-(2*pi*f').^2); + end +end +#+end_src + +#+begin_src matlab +%% Compute transfer function in cartesian frame +G_compl_b = pagemtimes(J_inv, pagemtimes(G_raw, Jf_inv)); +#+end_src + +#+begin_src matlab +colors = colororder; +figure; +hold on; +% plot(freqs, abs(squeeze(G_compl(1,1,:))), 'color', colors(1,:), 'DisplayName', '$C_x/F_x$') +% plot(freqs, abs(squeeze(G_compl(2,2,:))), 'color', colors(2,:), 'DisplayName', '$C_y/F_y$') +% plot(freqs, abs(squeeze(G_compl(3,3,:))), 'color', colors(3,:), 'DisplayName', '$C_z/F_z$') +plot(f, abs(squeeze(G_compl_b(1,1,:))), '-', 'color', colors(1,:), 'DisplayName', '$C_x/F_x$') +plot(f, abs(squeeze(G_compl_b(2,2,:))), '-', 'color', colors(2,:), 'DisplayName', '$C_y/F_y$') +plot(f, abs(squeeze(G_compl_b(3,3,:))), '-', 'color', colors(3,:), 'DisplayName', '$C_z/F_z$') +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude [m/N]'); +xlim([30, 300]); ylim([1e-9, 2e-6]); +leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; +#+end_src + +*** Load Data +#+begin_src matlab +% Load data +data = [load('data/Measurement1.mat'), ... + load('data/Measurement2.mat'), ... + load('data/Measurement3.mat'), ... + load('data/Measurement4.mat'), ... + load('data/Measurement5.mat'), ... + load('data/Measurement6.mat'), ... + load('data/Measurement7.mat'), ... + load('data/Measurement8.mat'), ... + load('data/Measurement9.mat'), ... + load('data/Measurement10.mat')]; +#+end_src + +#+begin_src matlab +% Frequency Vector +freqs = m3.FFT1_H1_1_13_X_Val; +w = 2*pi*freqs; + +% 12 outputs, 10 inputs +G_raw = zeros(12, 10, length(freqs)); + +for j = 1:10 + for i = 1:12 + G_raw(i,j,:) = data(j).(sprintf("FFT1_H1_%i_13_Y_ReIm", i))./(-w.^2); + end +end +#+end_src + +#+begin_src matlab +%% Compute transfer function in cartesian frame +G_compl = pagemtimes(J_inv, pagemtimes(G_raw, Jf_inv)); +#+end_src + +#+begin_src matlab +colors = colororder; +figure; +hold on; +plot(freqs, abs(squeeze(G_compl(1,1,:))), 'color', colors(1,:), 'DisplayName', '$C_x/F_x$') +plot(freqs, abs(squeeze(G_compl(2,2,:))), 'color', colors(2,:), 'DisplayName', '$C_y/F_y$') +plot(freqs, abs(squeeze(G_compl(3,3,:))), 'color', colors(3,:), 'DisplayName', '$C_z/F_z$') +plot(f, abs(squeeze(G_compl_b(1,1,:))), '.', 'color', colors(1,:), 'DisplayName', '$C_x/F_x$') +plot(f, abs(squeeze(G_compl_b(2,2,:))), '.', 'color', colors(2,:), 'DisplayName', '$C_y/F_y$') +plot(f, abs(squeeze(G_compl_b(3,3,:))), '.', 'color', colors(3,:), 'DisplayName', '$C_z/F_z$') +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude [m/N]'); +xlim([30, 300]); ylim([1e-9, 2e-6]); +leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; +#+end_src + +#+begin_src matlab +colors = colororder; +figure; +hold on; +plot(freqs, abs(squeeze(G_compl(4,4,:))), 'color', colors(1,:), 'DisplayName', '$R_x/M_x$') +plot(freqs, abs(squeeze(G_compl(5,5,:))), 'color', colors(2,:), 'DisplayName', '$R_y/M_y$') +plot(freqs, abs(squeeze(G_compl(6,6,:))), 'color', colors(3,:), 'DisplayName', '$R_z/M_z$') +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude [m/N]'); +xlim([30, 300]); ylim([5e-7, 5e-5]); +leg = legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1); +leg.ItemTokenSize(1) = 15; +#+end_src + +*** Diagonal Dynamics +#+begin_src matlab +figure; +ax1 = subplot(2,1,1); +hold on; +plot(freqs, abs(squeeze(G(1,1,:)), '-', 'DisplayName', '$D_x/F_x$') +plot(freqs, abs(squeeze(G(2,2,:)), '-', 'DisplayName', '$D_y/F_y$') +plot(freqs, abs(squeeze(G(3,3,:)), '-', 'DisplayName', '$D_z/F_z$') +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude [m/N]'); +ylim([1e-9, 2e-6]); +legend('location', 'southwest'); +xlim([30, 300]); +#+end_src + +#+begin_src matlab +figure; +ax1 = subplot(2,1,1); +hold on; +plot(freqs, abs(squeeze(G(4,4,:)), '-', 'DisplayName', '$R_x/M_x$') +plot(freqs, abs(squeeze(G(5,5,:)), '-', 'DisplayName', '$R_y/M_y$') +plot(freqs, abs(squeeze(G(6,6,:)), '-', 'DisplayName', '$R_z/M_z$') +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude [rad/Nm]'); +ylim([1e-7, 2e-4]); +legend('location', 'southwest'); +xlim([30, 300]); +#+end_src + +*** Equivalent Stiffness and Mass Estimation + +#+begin_src matlab +K = [1e7, 1e7, 2e8, 5e7, 3e7, 2e7]; +f_res = [125, 135, 390, 335, 335, 160]; +#+end_src + +#+begin_src matlab +M = [20, 20, 20, 11, 7, 20]; +f_res_est = sqrt(K./M)./(2*pi); +#+end_src + +Here is the inertia / stiffness to the granite that can represent the micro-station compliance dynamics: +#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*) +data2orgtable([K'], {'x', 'y', 'z', 'Rx', 'Ry', 'Rz'}, {'Stiffness', 'Inertia'}, ' %.1g '); +#+end_src + + +#+RESULTS: +| Stiffness | Inertia | +|-----------+-------------| +| x | 10000000.0 | +| y | 10000000.0 | +| z | 200000000.0 | +| Rx | 50000000.0 | +| Ry | 30000000.0 | +| Rz | 20000000.0 | + +** Compare with the Model #+begin_src matlab %% Initialize simulation with default parameters (flexible elements) initializeGround(); @@ -946,6 +1354,176 @@ legend('location', 'northwest'); - [ ] Comparison with the estimated (or measured) compliance +#+begin_src matlab + figure; + ax1 = subplot(2,1,1); + hold on; + plot(freqs, abs(squeeze(G(1,1,:))./(-w.^2)), '.') + plot(freqs, abs(squeeze(G(2,2,:))./(-w.^2)), '.') + plot(freqs, abs(squeeze(G(3,3,:))./(-w.^2)), '.') + set(gca,'ColorOrderIndex',1); + plot(freqs, abs(squeeze(freqresp(Gm(1,1,:), freqs, 'Hz'))), '-') + plot(freqs, abs(squeeze(freqresp(Gm(2,2,:), freqs, 'Hz'))), '-') + plot(freqs, abs(squeeze(freqresp(Gm(3,3,:), freqs, 'Hz'))), '-') + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + ylabel('Magnitude [m/N]'); set(gca, 'XTickLabel',[]); + ylim([1e-9, 2e-6]); + + ax2 = subplot(2,1,2); + hold on; + plot(freqs, 180/pi*angle(squeeze(G(1,1,:))./(-w.^2)), '.', 'DisplayName', '$D_x/F_x$') + plot(freqs, 180/pi*angle(squeeze(G(2,2,:))./(-w.^2)), '.', 'DisplayName', '$D_y/F_y$') + plot(freqs, 180/pi*angle(squeeze(G(3,3,:))./(-w.^2)), '.', 'DisplayName', '$D_z/F_z$') + set(gca,'ColorOrderIndex',1); + plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(1,1,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off') + plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(2,2,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off') + plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(3,3,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off') + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); + xlabel('Freqency [Hz]'); ylabel('Phase [deg]'); + ylim([-180, 180]); + yticks([-180, -90, 0, 90, 180]); + legend('location', 'southwest'); + + linkaxes([ax1,ax2],'x'); + xlim([30, 300]); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/compliance_diagonal_translations_comp_model.pdf', 'width', 'full', 'height', 'full'); +#+end_src + +#+name: fig:compliance_diagonal_translations_comp_model +#+caption: Dynamics from Forces to Translations +#+RESULTS: +[[file:figs/compliance_diagonal_translations_comp_model.png]] + +#+begin_src matlab + figure; + ax1 = subplot(2,1,1); + hold on; + plot(freqs, abs(squeeze(G(4,4,:))./(-w.^2)), '.') + plot(freqs, abs(squeeze(G(5,5,:))./(-w.^2)), '.') + plot(freqs, abs(squeeze(G(6,6,:))./(-w.^2)), '.') + set(gca,'ColorOrderIndex',1); + plot(freqs, abs(squeeze(freqresp(Gm(4,4,:), freqs, 'Hz'))), '-') + plot(freqs, abs(squeeze(freqresp(Gm(5,5,:), freqs, 'Hz'))), '-') + plot(freqs, abs(squeeze(freqresp(Gm(6,6,:), freqs, 'Hz'))), '-') + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + ylabel('Magnitude [rad/Nm]'); set(gca, 'XTickLabel',[]); + % ylim([1e-9, 2e-6]); + + ax2 = subplot(2,1,2); + hold on; + plot(freqs, 180/pi*angle(squeeze(G(4,4,:))./(-w.^2)), '.', 'DisplayName', '$R_x/M_x$') + plot(freqs, 180/pi*angle(squeeze(G(5,5,:))./(-w.^2)), '.', 'DisplayName', '$R_y/M_y$') + plot(freqs, 180/pi*angle(squeeze(G(6,6,:))./(-w.^2)), '.', 'DisplayName', '$R_z/M_z$') + set(gca,'ColorOrderIndex',1); + plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(4,4,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off') + plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(5,5,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off') + plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(6,6,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off') + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); + xlabel('Freqency [Hz]'); ylabel('Phase [deg]'); + ylim([-180, 180]); + yticks([-180, -90, 0, 90, 180]); + legend('location', 'southwest'); + + linkaxes([ax1,ax2],'x'); + xlim([30, 300]); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/compliance_diagonal_rotations_comp_model.pdf', 'width', 'full', 'height', 'full'); +#+end_src + +#+name: fig:compliance_diagonal_rotations_comp_model +#+caption: Dynamics from Torques to Rotations +#+RESULTS: +[[file:figs/compliance_diagonal_rotations_comp_model.png]] + +| | Stiffness | Unit | +|-----------+-----------+----------| +| $K_x$ | 1e7 | [N/m] | +| $K_y$ | 1e7 | [N/m] | +| $K_z$ | 2e8 | [N/m] | +| $K_{R_x}$ | 5e7 | [Nm/rad] | +| $K_{R_y}$ | 3e7 | [Nm/rad] | +| $K_{R_z}$ | 2e7 | [Nm/rad] | + +*** Coupling Dynamics +#+begin_src matlab + figure; + ax1 = subplot(2,1,1); + hold on; + plot(freqs, abs(squeeze(G(1,1,:))./(-w.^2)), '.') + plot(freqs, abs(squeeze(G(2,1,:))./(-w.^2)), '.') + plot(freqs, abs(squeeze(G(3,1,:))./(-w.^2)), '.') + set(gca,'ColorOrderIndex',1); + plot(freqs, abs(squeeze(freqresp(Gm(1,1,:), freqs, 'Hz'))), '-') + plot(freqs, abs(squeeze(freqresp(Gm(2,1,:), freqs, 'Hz'))), '-') + plot(freqs, abs(squeeze(freqresp(Gm(3,1,:), freqs, 'Hz'))), '-') + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + ylabel('Magnitude [m/N]'); set(gca, 'XTickLabel',[]); + ylim([1e-9, 2e-6]); + + ax2 = subplot(2,1,2); + hold on; + plot(freqs, 180/pi*angle(squeeze(G(1,1,:))./(-w.^2)), '.', 'DisplayName', '$D_x/F_x$') + plot(freqs, 180/pi*angle(squeeze(G(2,1,:))./(-w.^2)), '.', 'DisplayName', '$D_y/F_x$') + plot(freqs, 180/pi*angle(squeeze(G(3,1,:))./(-w.^2)), '.', 'DisplayName', '$D_z/F_x$') + set(gca,'ColorOrderIndex',1); + plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(1,1,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off') + plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(2,1,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off') + plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(3,1,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off') + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); + xlabel('Freqency [Hz]'); ylabel('Phase [deg]'); + ylim([-180, 180]); + yticks([-180, -90, 0, 90, 180]); + legend('location', 'southwest'); + + linkaxes([ax1,ax2],'x'); + xlim([30, 300]); +#+end_src + +#+begin_src matlab + figure; + ax1 = subplot(2,1,1); + hold on; + plot(freqs, abs(squeeze(G(5,1,:))./(-w.^2)), '.') + plot(freqs, abs(squeeze(G(4,2,:))./(-w.^2)), '.') + set(gca,'ColorOrderIndex',1); + plot(freqs, abs(squeeze(freqresp(Gm(5,1,:), freqs, 'Hz'))), '-') + plot(freqs, abs(squeeze(freqresp(Gm(4,2,:), freqs, 'Hz'))), '-') + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); + ylabel('Magnitude [m/N]'); set(gca, 'XTickLabel',[]); + ylim([1e-9, 2e-6]); + + ax2 = subplot(2,1,2); + hold on; + plot(freqs, 180/pi*angle(squeeze(G(5,1,:))./(-w.^2)), '.', 'DisplayName', '$R_y/F_x$') + plot(freqs, 180/pi*angle(squeeze(G(4,2,:))./(-w.^2)), '.', 'DisplayName', '$R_x/F_y$') + set(gca,'ColorOrderIndex',1); + plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(5,1,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off') + plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(4,2,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off') + hold off; + set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); + xlabel('Freqency [Hz]'); ylabel('Phase [deg]'); + ylim([-180, 180]); + yticks([-180, -90, 0, 90, 180]); + legend('location', 'southwest'); + + linkaxes([ax1,ax2],'x'); + xlim([30, 300]); +#+end_src + + + ** Conclusion For such a complex system, we believe that the Simscape Model represents the dynamics of the system with enough fidelity.