Add compliance measurement

This commit is contained in:
Thomas Dehaeze 2024-10-30 19:21:34 +01:00
parent f7c5716c04
commit 151ae07072

View File

@ -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.