Analysis of the resonance frequencies of the APA

test-bench-apa300ml.org

@@ -1612,5 +1612,872 @@ It is the expected behavior as shown in the Figure [[fig:souleille18_results]] (
 #+caption: Results obtained in cite:souleille18_concep_activ_mount_space_applic
 [[file:figs/souleille18_results.png]]
 
-* Bibliography                                                       :ignore:
+* Test Bench APA300ML - Simscape Model
+** Introduction
+** Matlab Init                                             :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+<<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+<<matlab-init>>
+#+end_src
+
+#+begin_src matlab :tangle no
+addpath('matlab/test_bench_apa300ml/');
+#+end_src
+
+#+begin_src matlab :eval no
+addpath('test_bench_apa300ml/');
+#+end_src
+
+#+begin_src matlab
+open('test_bench_apa300ml.slx')
+#+end_src
+
+** Nano Hexapod object
+#+begin_src matlab
+n_hexapod = struct();
+#+end_src
+
+*** APA - 2 DoF
+#+begin_src matlab
+n_hexapod.actuator = struct();
+
+n_hexapod.actuator.type = 1;
+
+n_hexapod.actuator.k  = ones(6,1)*0.35e6; % [N/m]
+n_hexapod.actuator.ke = ones(6,1)*1.5e6; % [N/m]
+n_hexapod.actuator.ka = ones(6,1)*43e6; % [N/m]
+
+n_hexapod.actuator.c  = ones(6,1)*3e1; % [N/(m/s)]
+n_hexapod.actuator.ce = ones(6,1)*1e1; % [N/(m/s)]
+n_hexapod.actuator.ca = ones(6,1)*1e1; % [N/(m/s)]
+
+n_hexapod.actuator.Leq = ones(6,1)*0.056; % [m]
+
+n_hexapod.actuator.Ga = ones(6,1)*1; % Actuator gain [N/V]
+n_hexapod.actuator.Gs = ones(6,1)*1; % Sensor gain [V/m]
+#+end_src
+
+*** APA - Flexible Frame
+#+begin_src matlab
+n_hexapod.actuator.type = 2;
+
+n_hexapod.actuator.K = readmatrix('APA300ML_b_mat_K.CSV'); % Stiffness Matrix
+n_hexapod.actuator.M = readmatrix('APA300ML_b_mat_M.CSV'); % Mass Matrix
+n_hexapod.actuator.xi = 0.01; % Damping ratio
+n_hexapod.actuator.P = extractNodes('APA300ML_b_out_nodes_3D.txt'); % Node coordinates [m]
+
+n_hexapod.actuator.ks = 235e6; % Stiffness of one stack [N/m]
+n_hexapod.actuator.cs = 1e1; % Stiffness of one stack [N/m]
+
+n_hexapod.actuator.Ga = ones(6,1)*1; % Actuator gain [N/V]
+n_hexapod.actuator.Gs = ones(6,1)*1; % Sensor gain [V/m]
+#+end_src
+
+*** APA - Fully Flexible
+#+begin_src matlab
+n_hexapod.actuator.type = 3;
+
+n_hexapod.actuator.K = readmatrix('APA300ML_full_mat_K.CSV'); % Stiffness Matrix
+n_hexapod.actuator.M = readmatrix('APA300ML_full_mat_M.CSV'); % Mass Matrix
+n_hexapod.actuator.xi = 0.01; % Damping ratio
+n_hexapod.actuator.P = extractNodes('APA300ML_full_out_nodes_3D.txt'); % Node coordiantes [m]
+
+n_hexapod.actuator.Ga = ones(6,1)*1; % Actuator gain [N/V]
+n_hexapod.actuator.Gs = ones(6,1)*1; % Sensor gain [V/m]
+#+end_src
+
+** Identification
+#+begin_src matlab
+%% Options for Linearized
+options = linearizeOptions;
+options.SampleTime = 0;
+
+%% Name of the Simulink File
+mdl = 'test_bench_apa300ml';
+
+%% Input/Output definition
+clear io; io_i = 1;
+io(io_i) = linio([mdl, '/Va'],  1, 'openinput');  io_i = io_i + 1; % Actuator Voltage
+io(io_i) = linio([mdl, '/Vs'],  1, 'openoutput'); io_i = io_i + 1; % Sensor Voltage
+io(io_i) = linio([mdl, '/dL'],  1, 'openoutput'); io_i = io_i + 1; % Relative Motion Outputs
+io(io_i) = linio([mdl, '/z'],   1, 'openoutput'); io_i = io_i + 1; % Vertical Motion
+
+%% Run the linearization
+Ga = linearize(mdl, io, 0.0, options);
+Ga.InputName  = {'Va'};
+Ga.OutputName = {'Vs', 'dL', 'z'};
+#+end_src
+
+#+begin_src matlab :exports none
+freqs = logspace(1, 3, 1000);
+
+figure;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+plot(freqs, abs(squeeze(freqresp(Ga('Vs', 'Va'), freqs, 'Hz'))), 'DisplayName', '')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]);
+hold off;
+legend('location', 'southwest');
+
+ax2 = nexttile;
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(Ga('Vs', 'Va'), freqs, 'Hz'))))
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+yticks(-360:45:360);
+ylim([-180, 180])
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+begin_src matlab :exports none
+freqs = logspace(1, 3, 1000);
+
+figure;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+plot(freqs, abs(squeeze(freqresp(Ga('dL', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Encoder')
+plot(freqs, abs(squeeze(freqresp(Ga('z', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Interferometer')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]);
+hold off;
+legend('location', 'southwest');
+
+ax2 = nexttile;
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(Ga('dL', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(Ga('z', 'Va'), freqs, 'Hz'))))
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+yticks(-360:45:360);
+ylim([-180, 180])
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+** Compare 2-DoF with flexible
+*** APA - 2 DoF
+#+begin_src matlab
+n_hexapod = struct();
+
+n_hexapod.actuator = struct();
+
+n_hexapod.actuator.type = 1;
+
+n_hexapod.actuator.k  = ones(6,1)*0.35e6; % [N/m]
+n_hexapod.actuator.ke = ones(6,1)*1.5e6; % [N/m]
+n_hexapod.actuator.ka = ones(6,1)*43e6; % [N/m]
+
+n_hexapod.actuator.c  = ones(6,1)*3e1; % [N/(m/s)]
+n_hexapod.actuator.ce = ones(6,1)*1e1; % [N/(m/s)]
+n_hexapod.actuator.ca = ones(6,1)*1e1; % [N/(m/s)]
+
+n_hexapod.actuator.Leq = ones(6,1)*0.056; % [m]
+
+n_hexapod.actuator.Ga = ones(6,1)*-2.15; % Actuator gain [N/V]
+n_hexapod.actuator.Gs = ones(6,1)*2.305e-08; % Sensor gain [V/m]
+#+end_src
+
+#+begin_src matlab
+G_2dof = linearize(mdl, io, 0.0, options);
+G_2dof.InputName  = {'Va'};
+G_2dof.OutputName = {'Vs', 'dL', 'z'};
+#+end_src
+
+*** APA - Fully Flexible
+#+begin_src matlab
+n_hexapod = struct();
+
+n_hexapod.actuator.type = 3;
+
+n_hexapod.actuator.K = readmatrix('APA300ML_full_mat_K.CSV'); % Stiffness Matrix
+n_hexapod.actuator.M = readmatrix('APA300ML_full_mat_M.CSV'); % Mass Matrix
+n_hexapod.actuator.xi = 0.01; % Damping ratio
+n_hexapod.actuator.P = extractNodes('APA300ML_full_out_nodes_3D.txt'); % Node coordiantes [m]
+
+n_hexapod.actuator.Ga = ones(6,1)*1; % Actuator gain [N/V]
+n_hexapod.actuator.Gs = ones(6,1)*1; % Sensor gain [V/m]
+#+end_src
+
+#+begin_src matlab
+G_flex = linearize(mdl, io, 0.0, options);
+G_flex.InputName  = {'Va'};
+G_flex.OutputName = {'Vs', 'dL', 'z'};
+#+end_src
+
+*** Comparison
+
+#+begin_src matlab :exports none
+freqs = logspace(1, 4, 1000);
+
+figure;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+plot(freqs, abs(squeeze(freqresp(G_2dof('Vs', 'Va'), freqs, 'Hz'))), 'DisplayName', '$G_a$')
+plot(freqs, abs(squeeze(freqresp(G_flex('Vs', 'Va'), freqs, 'Hz'))), 'DisplayName', '$G_s$')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]);
+hold off;
+legend('location', 'southwest');
+
+ax2 = nexttile;
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(G_2dof('Vs', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(G_flex('Vs', 'Va'), freqs, 'Hz'))))
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+yticks(-360:45:360);
+ylim([-180, 180])
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+begin_src matlab :exports none
+freqs = logspace(1, 4, 1000);
+
+figure;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+plot(freqs, abs(squeeze(freqresp(G_2dof('dL', 'Va'), freqs, 'Hz'))), 'DisplayName', '$G_a$')
+plot(freqs, abs(squeeze(freqresp(G_flex('dL', 'Va'), freqs, 'Hz'))), 'DisplayName', '$G_s$')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ylabel('Amplitude $d_L/V_a$ [m/V]'); set(gca, 'XTickLabel',[]);
+hold off;
+legend('location', 'southwest');
+
+ax2 = nexttile;
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(G_2dof('dL', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(G_flex('dL', 'Va'), freqs, 'Hz'))))
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+yticks(-360:45:360);
+ylim([-180, 180])
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+begin_src matlab :exports none
+freqs = logspace(1, 4, 1000);
+
+figure;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+plot(freqs, abs(squeeze(freqresp(G_2dof('z', 'Va'), freqs, 'Hz'))), 'DisplayName', '$G_a$')
+plot(freqs, abs(squeeze(freqresp(G_flex('z', 'Va'), freqs, 'Hz'))), 'DisplayName', '$G_s$')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ylabel('Amplitude $z/V_a$ [m/V]'); set(gca, 'XTickLabel',[]);
+hold off;
+legend('location', 'southwest');
+
+ax2 = nexttile;
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(G_2dof('z', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(G_flex('z', 'Va'), freqs, 'Hz'))))
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+yticks(-360:45:360);
+ylim([-180, 180])
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+* Test Bench Struts - Simscape Model
+** Introduction
+** Matlab Init                                             :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+<<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+<<matlab-init>>
+#+end_src
+
+#+begin_src matlab :tangle no
+addpath('matlab/');
+addpath('matlab/test_bench_struts/');
+addpath('matlab/png/');
+#+end_src
+
+#+begin_src matlab :eval no
+addpath('test_bench_struts/');
+addpath('png/');
+#+end_src
+
+#+begin_src matlab
+open('test_bench_struts.slx')
+#+end_src
+
+** Nano Hexapod object
+#+begin_src matlab
+n_hexapod = struct();
+#+end_src
+
+*** Flexible Joint - Bot
+#+begin_src matlab
+n_hexapod.flex_bot = struct();
+
+n_hexapod.flex_bot.type = 1; % 1: 2dof / 2: 3dof / 3: 4dof
+
+n_hexapod.flex_bot.kRx = ones(6,1)*5; % X bending stiffness [Nm/rad]
+n_hexapod.flex_bot.kRy = ones(6,1)*5; % Y bending stiffness [Nm/rad]
+n_hexapod.flex_bot.kRz = ones(6,1)*260; % Torsionnal stiffness [Nm/rad]
+n_hexapod.flex_bot.kz  = ones(6,1)*1e8; % Axial stiffness [N/m]
+
+n_hexapod.flex_bot.cRx = ones(6,1)*0.1; % [Nm/(rad/s)]
+n_hexapod.flex_bot.cRy = ones(6,1)*0.1; % [Nm/(rad/s)]
+n_hexapod.flex_bot.cRz = ones(6,1)*0.1; % [Nm/(rad/s)]
+n_hexapod.flex_bot.cz  = ones(6,1)*1e2; %[N/(m/s)]
+#+end_src
+
+*** Flexible Joint - Top
+#+begin_src matlab
+n_hexapod.flex_top = struct();
+
+n_hexapod.flex_top.type = 2; % 1: 2dof / 2: 3dof / 3: 4dof
+
+n_hexapod.flex_top.kRx = ones(6,1)*5; % X bending stiffness [Nm/rad]
+n_hexapod.flex_top.kRy = ones(6,1)*5; % Y bending stiffness [Nm/rad]
+n_hexapod.flex_top.kRz = ones(6,1)*260; % Torsionnal stiffness [Nm/rad]
+n_hexapod.flex_top.kz  = ones(6,1)*1e8; % Axial stiffness [N/m]
+
+n_hexapod.flex_top.cRx = ones(6,1)*0.1; % [Nm/(rad/s)]
+n_hexapod.flex_top.cRy = ones(6,1)*0.1; % [Nm/(rad/s)]
+n_hexapod.flex_top.cRz = ones(6,1)*0.1; % [Nm/(rad/s)]
+n_hexapod.flex_top.cz  = ones(6,1)*1e2; %[N/(m/s)]
+#+end_src
+
+*** APA - 2 DoF
+#+begin_src matlab
+n_hexapod.actuator = struct();
+
+n_hexapod.actuator.type = 1;
+
+n_hexapod.actuator.k  = ones(6,1)*0.35e6; % [N/m]
+n_hexapod.actuator.ke = ones(6,1)*1.5e6; % [N/m]
+n_hexapod.actuator.ka = ones(6,1)*43e6; % [N/m]
+
+n_hexapod.actuator.c  = ones(6,1)*3e1; % [N/(m/s)]
+n_hexapod.actuator.ce = ones(6,1)*1e1; % [N/(m/s)]
+n_hexapod.actuator.ca = ones(6,1)*1e1; % [N/(m/s)]
+
+n_hexapod.actuator.Leq = ones(6,1)*0.056; % [m]
+
+n_hexapod.actuator.Ga = ones(6,1)*1; % Actuator gain [N/V]
+n_hexapod.actuator.Gs = ones(6,1)*1; % Sensor gain [V/m]
+#+end_src
+
+*** APA - Flexible Frame
+#+begin_src matlab
+n_hexapod.actuator.type = 2;
+
+n_hexapod.actuator.K = readmatrix('APA300ML_b_mat_K.CSV'); % Stiffness Matrix
+n_hexapod.actuator.M = readmatrix('APA300ML_b_mat_M.CSV'); % Mass Matrix
+n_hexapod.actuator.xi = 0.01; % Damping ratio
+n_hexapod.actuator.P = extractNodes('APA300ML_b_out_nodes_3D.txt'); % Node coordinates [m]
+
+n_hexapod.actuator.ks = 235e6; % Stiffness of one stack [N/m]
+n_hexapod.actuator.cs = 1e1; % Stiffness of one stack [N/m]
+
+n_hexapod.actuator.Ga = ones(6,1)*1; % Actuator gain [N/V]
+n_hexapod.actuator.Gs = ones(6,1)*1; % Sensor gain [V/m]
+#+end_src
+
+*** APA - Fully Flexible
+#+begin_src matlab
+n_hexapod.actuator.type = 3;
+
+n_hexapod.actuator.K = readmatrix('APA300ML_full_mat_K.CSV'); % Stiffness Matrix
+n_hexapod.actuator.M = readmatrix('APA300ML_full_mat_M.CSV'); % Mass Matrix
+n_hexapod.actuator.xi = 0.01; % Damping ratio
+n_hexapod.actuator.P = extractNodes('APA300ML_full_out_nodes_3D.txt'); % Node coordiantes [m]
+
+n_hexapod.actuator.Ga = ones(6,1)*1; % Actuator gain [N/V]
+n_hexapod.actuator.Gs = ones(6,1)*1; % Sensor gain [V/m]
+#+end_src
+
+
+** Identification
+#+begin_src matlab
+%% Options for Linearized
+options = linearizeOptions;
+options.SampleTime = 0;
+
+%% Name of the Simulink File
+mdl = 'test_bench_struts';
+
+%% Input/Output definition
+clear io; io_i = 1;
+io(io_i) = linio([mdl, '/Va'],  1, 'openinput');  io_i = io_i + 1; % Actuator Voltage
+io(io_i) = linio([mdl, '/Vs'],  1, 'openoutput'); io_i = io_i + 1; % Sensor Voltage
+io(io_i) = linio([mdl, '/dL'],  1, 'openoutput'); io_i = io_i + 1; % Relative Motion Outputs
+io(io_i) = linio([mdl, '/z'],   1, 'openoutput'); io_i = io_i + 1; % Vertical Motion
+
+%% Run the linearization
+Gs = linearize(mdl, io, 0.0, options);
+Gs.InputName  = {'Va'};
+Gs.OutputName = {'Vs', 'dL', 'z'};
+#+end_src
+
+#+begin_src matlab :exports none
+freqs = logspace(1, 3, 1000);
+
+figure;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+plot(freqs, abs(squeeze(freqresp(Gs('Vs', 'Va'), freqs, 'Hz'))), 'DisplayName', '')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]);
+hold off;
+legend('location', 'southwest');
+
+ax2 = nexttile;
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gs('Vs', 'Va'), freqs, 'Hz'))))
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+yticks(-360:45:360);
+ylim([-180, 180])
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+begin_src matlab :exports none
+freqs = logspace(1, 4, 1000);
+
+figure;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+plot(freqs, abs(squeeze(freqresp(Gs('dL', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Encoder')
+plot(freqs, abs(squeeze(freqresp(Gs('z', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Interferometer')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]);
+hold off;
+legend('location', 'southwest');
+
+ax2 = nexttile;
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gs('dL', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gs('z', 'Va'), freqs, 'Hz'))))
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+yticks(-360:45:360);
+ylim([-180, 180])
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+** Compare flexible joints
+*** Perfect
+#+begin_src matlab
+n_hexapod.flex_bot.type = 1; % 1: 2dof / 2: 3dof / 3: 4dof
+n_hexapod.flex_top.type = 2; % 1: 2dof / 2: 3dof / 3: 4dof
+#+end_src
+
+#+begin_src matlab
+Gp = linearize(mdl, io, 0.0, options);
+Gp.InputName  = {'Va'};
+Gp.OutputName = {'Vs', 'dL', 'z'};
+#+end_src
+
+*** Top Flexible
+#+begin_src matlab
+n_hexapod.flex_bot.type = 1; % 1: 2dof / 2: 3dof / 3: 4dof
+n_hexapod.flex_top.type = 3; % 1: 2dof / 2: 3dof / 3: 4dof
+#+end_src
+
+#+begin_src matlab
+Gt = linearize(mdl, io, 0.0, options);
+Gt.InputName  = {'Va'};
+Gt.OutputName = {'Vs', 'dL', 'z'};
+#+end_src
+
+*** Bottom Flexible
+#+begin_src matlab
+n_hexapod.flex_bot.type = 3; % 1: 2dof / 2: 3dof / 3: 4dof
+n_hexapod.flex_top.type = 2; % 1: 2dof / 2: 3dof / 3: 4dof
+#+end_src
+
+#+begin_src matlab
+Gb = linearize(mdl, io, 0.0, options);
+Gb.InputName  = {'Va'};
+Gb.OutputName = {'Vs', 'dL', 'z'};
+#+end_src
+
+*** Both Flexible
+#+begin_src matlab
+n_hexapod.flex_bot.type = 3; % 1: 2dof / 2: 3dof / 3: 4dof
+n_hexapod.flex_top.type = 3; % 1: 2dof / 2: 3dof / 3: 4dof
+#+end_src
+
+#+begin_src matlab
+Gf = linearize(mdl, io, 0.0, options);
+Gf.InputName  = {'Va'};
+Gf.OutputName = {'Vs', 'dL', 'z'};
+#+end_src
+
+*** Comparison
+#+begin_src matlab :exports none
+freqs = logspace(1, 4, 1000);
+
+figure;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+plot(freqs, abs(squeeze(freqresp(Gp('Vs', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Perfect')
+plot(freqs, abs(squeeze(freqresp(Gt('Vs', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Top')
+plot(freqs, abs(squeeze(freqresp(Gb('Vs', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Bot')
+plot(freqs, abs(squeeze(freqresp(Gf('Vs', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Flex')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]);
+hold off;
+legend('location', 'southwest');
+
+ax2 = nexttile;
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gp('Vs', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gt('Vs', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gb('Vs', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gf('Vs', 'Va'), freqs, 'Hz'))))
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+yticks(-360:45:360);
+ylim([-180, 180])
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+begin_src matlab :exports none
+freqs = logspace(1, 4, 1000);
+
+figure;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+plot(freqs, abs(squeeze(freqresp(Gp('dL', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Perfect')
+plot(freqs, abs(squeeze(freqresp(Gt('dL', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Top')
+plot(freqs, abs(squeeze(freqresp(Gb('dL', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Bot')
+plot(freqs, abs(squeeze(freqresp(Gf('dL', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Flex')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ylabel('Amplitude $d_L/V_a$ [m/V]'); set(gca, 'XTickLabel',[]);
+hold off;
+legend('location', 'southwest');
+
+ax2 = nexttile;
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gp('dL', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gt('dL', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gb('dL', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gf('dL', 'Va'), freqs, 'Hz'))))
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+yticks(-360:45:360);
+ylim([-180, 180])
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+begin_src matlab :exports none
+freqs = logspace(1, 4, 1000);
+
+figure;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+plot(freqs, abs(squeeze(freqresp(Gp('z', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Perfect')
+plot(freqs, abs(squeeze(freqresp(Gt('z', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Top')
+plot(freqs, abs(squeeze(freqresp(Gb('z', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Bot')
+plot(freqs, abs(squeeze(freqresp(Gf('z', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Flex')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ylabel('Amplitude $z/V_a$ [m/V]'); set(gca, 'XTickLabel',[]);
+hold off;
+legend('location', 'southwest');
+
+ax2 = nexttile;
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gp('z', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gt('z', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gb('z', 'Va'), freqs, 'Hz'))))
+plot(freqs, 180/pi*angle(squeeze(freqresp(Gf('z', 'Va'), freqs, 'Hz'))))
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+yticks(-360:45:360);
+ylim([-180, 180])
+
+linkaxes([ax1,ax2],'x');
+#+end_src
+
+* Resonance frequencies - APA300ML
+** Introduction
+
+Three main resonances are foreseen to be problematic for the control of the APA300ML:
+- Mode in X-bending at 189Hz (Figure [[fig:mode_bending_x]])
+- Mode in Y-bending at 285Hz (Figure [[fig:mode_bending_y]])
+- Mode in Z-torsion at 400Hz (Figure [[fig:mode_torsion_z]])
+
+#+name: fig:mode_bending_x
+#+caption: X-bending mode (189Hz)
+#+attr_latex: :width 0.9\linewidth
+[[file:figs/mode_bending_x.gif]]
+
+#+name: fig:mode_bending_y
+#+caption: Y-bending mode (285Hz)
+#+attr_latex: :width 0.9\linewidth
+[[file:figs/mode_bending_y.gif]]
+
+#+name: fig:mode_torsion_z
+#+caption: Z-torsion mode (400Hz)
+#+attr_latex: :width 0.9\linewidth
+[[file:figs/mode_torsion_z.gif]]
+
+These modes are present when flexible joints are fixed to the ends of the APA300ML.
+
+In this section, we try to find the resonance frequency of these modes when one end of the APA is fixed and the other is free.
+
+** Matlab Init                                              :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+<<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+<<matlab-init>>
+#+end_src
+
+#+begin_src matlab :tangle no
+addpath('matlab/');
+addpath('matlab/mat/');
+#+end_src
+
+#+begin_src matlab :eval no
+addpath('mat/');
+#+end_src
+
+** Setup
+
+The measurement setup is shown in Figure [[fig:measurement_setup_torsion]].
+A Laser vibrometer is measuring the difference of motion of two points.
+The APA is excited with an instrumented hammer and the transfer function from the hammer to the measured rotation is computed.
+
+#+begin_note
+- Laser Doppler Vibrometer Polytec OFV512
+- Instrumented hammer
+#+end_note
+
+#+name: fig:measurement_setup_torsion
+#+caption: Measurement setup with a Laser Doppler Vibrometer and one instrumental hammer
+#+attr_latex: :width 0.7\linewidth
+[[file:figs/measurement_setup_torsion.png]]
+
+** Bending - X
+
+The setup to measure the X-bending motion is shown in Figure [[fig:measurement_setup_X_bending]].
+The APA is excited with an instrumented hammer having a solid metallic tip.
+The impact point is on the back-side of the APA aligned with the top measurement point.
+
+#+name: fig:measurement_setup_X_bending
+#+caption: X-Bending measurement setup
+#+attr_latex: :width 0.7\linewidth
+[[file:figs/measurement_setup_X_bending.png]]
+
+The data is loaded.
+#+begin_src matlab
+bending_X = load('apa300ml_bending_X_top.mat')
+#+end_src
+
+#+begin_src matlab
+Ts = bending_X.Track1_X_Resolution; % Sampling frequency [Hz]
+#+end_src
+
+The transfer function from the input force to the output "rotation" (difference between the two measured distances).
+#+begin_src matlab
+win = hann(ceil(1/Ts));
+[G_bending_X, f]  = tfestimate(bending_X.Track1, bending_X.Track2, win, [], [], 1/Ts);
+#+end_src
+
+The result is shown in Figure [[fig:apa300ml_meas_freq_bending_x]].
+
+The can clearly observe a nice peak at 280Hz, and then peaks at the odd "harmonics" (third "harmonic" at 840Hz, and fifth "harmonic" at 1400Hz).
+#+begin_src matlab :exports none
+figure;
+hold on;
+plot(f, abs(G_bending_X), 'k-');
+hold off;
+set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Amplitude');
+xlim([50, 2e3]); ylim([1e-5, 2e-1]);
+text(280, 5.5e-2,{'280Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
+text(840, 2.0e-3,{'840Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
+text(1400, 7.0e-3,{'1400Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+exportFig('figs/apa300ml_meas_freq_bending_x.pdf', 'width', 'wide', 'height', 'normal');
+#+end_src
+
+#+name: fig:apa300ml_meas_freq_bending_x
+#+caption: Obtained FRF for the X-bending
+#+RESULTS:
+[[file:figs/apa300ml_meas_freq_bending_x.png]]
+
+** Bending - Y
+
+The setup to measure the Y-bending is shown in Figure [[fig:measurement_setup_Y_bending]].
+
+The impact point of the instrumented hammer is located on the back surface of the top interface (on the back of the 2 measurements points).
+
+#+name: fig:measurement_setup_Y_bending
+#+caption: Y-Bending measurement setup
+#+attr_latex: :width 0.7\linewidth
+[[file:figs/measurement_setup_Y_bending.png]]
+
+The data is loaded, and the transfer function from the force to the measured rotation is computed.
+#+begin_src matlab
+bending_Y = load('apa300ml_bending_Y_top.mat')
+[G_bending_Y, ~]  = tfestimate(bending_Y.Track1, bending_Y.Track2, win, [], [], 1/Ts);
+#+end_src
+
+The results are shown in Figure [[fig:apa300ml_meas_freq_bending_y]].
+The main resonance is at 412Hz, and we also see the third "harmonic" at 1220Hz.
+
+#+begin_src matlab :exports none
+figure;
+hold on;
+plot(f, abs(G_bending_Y), 'k-');
+hold off;
+set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Amplitude');
+xlim([50, 2e3]); ylim([1e-5, 3e-2])
+text(412, 1.5e-2,{'412Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
+text(1218, 1.5e-2,{'1220Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+exportFig('figs/apa300ml_meas_freq_bending_y.pdf', 'width', 'wide', 'height', 'normal');
+#+end_src
+
+#+name: fig:apa300ml_meas_freq_bending_y
+#+caption: Obtained FRF for the Y-bending
+#+RESULTS:
+[[file:figs/apa300ml_meas_freq_bending_y.png]]
+
+** Torsion - Z
+
+Finally, we measure the Z-torsion resonance as shown in Figure [[fig:measurement_setup_torsion_bis]].
+
+The excitation is shown on the other side of the APA, on the side to excite the torsion motion.
+
+#+name: fig:measurement_setup_torsion_bis
+#+caption: Z-Torsion measurement setup
+#+attr_latex: :width 0.7\linewidth
+[[file:figs/measurement_setup_torsion_bis.png]]
+
+The data is loaded, and the transfer function computed.
+#+begin_src matlab
+torsion = load('apa300ml_torsion_left.mat')
+[G_torsion, ~]  = tfestimate(torsion_left.Track1, torsion_left.Track2, win, [], [], 1/Ts);
+#+end_src
+
+The results are shown in Figure [[fig:apa300ml_meas_freq_torsion_z]].
+We observe a first peak at 267Hz, which corresponds to the X-bending mode that was measured at 280Hz.
+And then a second peak at 415Hz, which corresponds to the X-bending mode that was measured at 412Hz.
+The mode in pure torsion is probably at higher frequency (peak around 1kHz?).
+#+begin_src matlab :exports none
+figure;
+hold on;
+plot(f, abs(G_torsion), 'k-');
+hold off;
+set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Amplitude');
+xlim([50, 2e3]); ylim([1e-5, 2e-2])
+text(415, 4.3e-3,{'415Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
+text(267, 8e-4,{'267Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+exportFig('figs/apa300ml_meas_freq_torsion_z.pdf', 'width', 'wide', 'height', 'normal');
+#+end_src
+
+#+name: fig:apa300ml_meas_freq_torsion_z
+#+caption: Obtained FRF for the Z-torsion
+#+RESULTS:
+[[file:figs/apa300ml_meas_freq_torsion_z.png]]
+
+** Compare
+The three measurements are shown in Figure [[fig:apa300ml_meas_freq_compare]].
+#+begin_src matlab :exports none
+figure;
+hold on;
+plot(f, abs(G_torsion),    'DisplayName', 'Torsion');
+plot(f, abs(G_bending_X),  'DisplayName', 'Bending - X');
+plot(f, abs(G_bending_Y),  'DisplayName', 'Bending - Y');
+hold off;
+set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Amplitude');
+xlim([50, 2e3]); ylim([1e-5, 1e-1]);
+legend('location', 'southeast');
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+exportFig('figs/apa300ml_meas_freq_compare.pdf', 'width', 'full', 'height', 'tall');
+#+end_src
+
+#+name: fig:apa300ml_meas_freq_compare
+#+caption: Obtained FRF - Comparison
+#+RESULTS:
+[[file:figs/apa300ml_meas_freq_compare.png]]
+
+** Conclusion
+
+#+name: tab:apa300ml_measured_modes_freq
+#+caption: Measured frequency of the modes
+#+attr_latex: :environment tabularx :width 0.3\linewidth :align lX
+#+attr_latex: :center t :booktabs t :float t
+| Mode      | Measured Frequency |
+|-----------+--------------------|
+| X-Bending | 280Hz              |
+| Y-Bending | 410Hz              |
+| Z-Torsion | ?                  |
+
+* Bibliography                                                        :ignore:
 #+latex: \printbibliography