diff --git a/matlab/test_apa_2_dynamics.m b/matlab/test_apa_2_dynamics.m index 90974f3..9631b5e 100644 --- a/matlab/test_apa_2_dynamics.m +++ b/matlab/test_apa_2_dynamics.m @@ -149,9 +149,6 @@ add_mass_cc.de = add_mass_cc.de - mean(add_mass_cc.de(add_mass_cc.t<11)); apa_k_oc = 9.8 * added_mass / (mean(add_mass_oc.de(add_mass_oc.t > 12 & add_mass_oc.t < 12.5)) - mean(add_mass_oc.de(add_mass_oc.t > 20 & add_mass_oc.t < 20.5))); apa_k_sc = 9.8 * added_mass / (mean(add_mass_cc.de(add_mass_cc.t > 12 & add_mass_cc.t < 12.5)) - mean(add_mass_cc.de(add_mass_cc.t > 20 & add_mass_cc.t < 20.5))); -%% Estimated coupling factor -sqrt(1 - apa_k_sc/apa_k_oc) - % Dynamics % <> @@ -493,7 +490,6 @@ 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) - Niter = 100; % Number of iteration. N = 2; % Order of approximation poles = [-25 - 1i*60, -25 + 1i*60]; % First get for the pole location @@ -501,7 +497,7 @@ poles = [-25 - 1i*60, -25 + 1i*60]; % First get for the pole location G_dL_id = {zeros(1,length(i_kept))}; % Identification just between two frequencies -f_keep = (f>20 & f<200); +f_keep = (f>5 & f<200); for i = 1:length(i_kept) %% Estimate resonance frequency and damping @@ -543,12 +539,11 @@ legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1); %% Root Locus of the APA300ML with Integral Force Feedback % Comparison between the computed root locus from the plant model and the root locus estimated from the damped plant pole identification -figure; gains = logspace(-1, 3, 1000); figure; hold on; -G_iff_poles = pole(G_iff_model); +G_iff_poles = pole(pade(G_iff_model)); i = imag(G_iff_poles) > 100; % Only keep relevant poles plot(real(G_iff_poles(i)), imag(G_iff_poles(i)), 'kx', ... 'DisplayName', '$g = 0$'); @@ -558,7 +553,7 @@ plot(real(G_iff_zeros(i)), imag(G_iff_zeros(i)), 'ko', ... 'HandleVisibility', 'off'); for g = gains - clpoles = pole(feedback(G_iff_model, g*K_iff, 1)); + clpoles = pole(feedback(pade(G_iff_model), g*K_iff, 1)); i = imag(clpoles) > 100; % Only keep relevant poles plot(real(clpoles(i)), imag(clpoles(i)), 'k.', ... 'HandleVisibility', 'off'); diff --git a/matlab/test_apa_3_model_2dof.m b/matlab/test_apa_3_model_2dof.m index 94efc21..e95bdb0 100644 --- a/matlab/test_apa_3_model_2dof.m +++ b/matlab/test_apa_3_model_2dof.m @@ -84,7 +84,7 @@ ga = -abs(enc_frf(i_f,1))./abs(evalfr(G_norm('de', 'u'), 1i*2*pi*fa)); % Estimation ot the Sensor Gain fs = 600; % Frequency where the two FRF should match [Hz] -[~, i_f] = min(abs(f - fs)) +[~, i_f] = min(abs(f - fs)); gs = -abs(iff_frf(i_f,1))./abs(evalfr(G_norm('Vs', 'u'), 1i*2*pi*fs))/ga; % Obtained Dynamics diff --git a/matlab/test_apa_3_simscape.m b/matlab/test_apa_3_simscape.m deleted file mode 100644 index 15d8db9..0000000 --- a/matlab/test_apa_3_simscape.m +++ /dev/null @@ -1,566 +0,0 @@ -%% Clear Workspace and Close figures -clear; close all; clc; - -%% Intialize Laplace variable -s = zpk('s'); - -%% Path for functions, data and scripts -addpath('./src/'); % Path for scripts -addpath('./mat/'); % Path for data - -addpath('./STEPS/'); % Path for Simscape Model - -%% Open Simscape Model -mdl = 'test_apa_simscape'; % Name of the Simulink File -open(mdl); % Open Simscape Model - -%% Colors for the figures -colors = colororder; - -% First Identification -% <> - -% The APA is first initialized with default parameters: - -%% Initialize the structure with default values -n_hexapod = struct(); -n_hexapod.actuator = initializeAPA(... - 'type', '2dof', ... - 'Ga', 1, ... % Actuator constant [N/V] - 'Gs', 1); % Sensor constant [V/m] - - - -% The transfer function from excitation voltage $V_a$ (before the amplification of $20$ due to the PD200 amplifier) to: -% 1. the sensor stack voltage $V_s$ -% 2. the measured displacement by the encoder $d_e$ - - -%% Input/Output definition -clear io; io_i = 1; -io(io_i) = linio([mdl, '/Va'], 1, 'openinput'); io_i = io_i + 1; % DAC Voltage -io(io_i) = linio([mdl, '/Vs'], 1, 'openoutput'); io_i = io_i + 1; % Sensor Voltage -io(io_i) = linio([mdl, '/de'], 1, 'openoutput'); io_i = io_i + 1; % Encoder - -%% Linearization options -opts = linearizeOptions; -opts.SampleTime = 0; - -%% Run the linearization -Ga = linearize(mdl, io, 0.0, opts); -Ga.InputName = {'Va'}; -Ga.OutputName = {'Vs', 'de', 'da'}; - - - -% The obtain dynamics are shown in Figure ref:fig:apa_model_bench_bode_vs and ref:fig:apa_model_bench_bode_dl_z. -% It can be seen that: -% - the shape of these bode plots are very similar to the one measured in Section ref:sec:dynamical_meas_apa expect from a change in gain and exact location of poles and zeros -% - there is a sign error for the transfer function from $V_a$ to $V_s$. -% This will be corrected by taking a negative "sensor gain". -% - the low frequency zero of the transfer function from $V_a$ to $V_s$ is minimum phase as expected. -% The measured FRF are showing non-minimum phase zero, but it is most likely due to measurements artifacts. - - -%% Bode plot of the transfer function from u to taum -freqs = logspace(1, 3, 1000); - -figure; -tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); - -ax1 = nexttile([2,1]); -hold on; -plot(freqs, abs(squeeze(freqresp(Ga('Vs', 'Va'), freqs, 'Hz'))), 'k-') -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]); -hold off; - -ax2 = nexttile; -hold on; -plot(freqs, 180/pi*angle(squeeze(freqresp(Ga('Vs', 'Va'), freqs, 'Hz'))), 'k-') -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); -hold off; -yticks(-360:45:360); -ylim([-180, 0]) - -linkaxes([ax1,ax2],'x'); -xlim([freqs(1), freqs(end)]); - - - -% #+name: fig:apa_model_bench_bode_vs -% #+caption: Bode plot of the transfer function from $V_a$ to $V_s$ -% #+RESULTS: -% [[file:figs/apa_model_bench_bode_vs.png]] - - -%% Bode plot of the transfer function from Va to de and da -freqs = logspace(1, 3, 1000); - -figure; -tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); - -ax1 = nexttile([2,1]); -hold on; -plot(freqs, abs(squeeze(freqresp(Ga('de', 'Va'), freqs, 'Hz'))), 'DisplayName', 'Encoder') -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -ylabel('Amplitude $d/V_a$ [m/V]'); set(gca, 'XTickLabel',[]); -hold off; -legend('location', 'southwest'); - -ax2 = nexttile; -hold on; -plot(freqs, 180/pi*angle(squeeze(freqresp(Ga('de', '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, 0]) - -linkaxes([ax1,ax2],'x'); - -% Identification Data -% Let's load the measured FRF from the DAC voltage to the measured encoder and to the sensor stack voltage. - -%% Load Data -load('meas_apa_frf.mat', 'f', 'Ts', 'enc_frf', 'iff_frf', 'apa_nums'); - -% 2DoF APA -% Let's initialize the APA as a 2DoF model with unity sensor and actuator gains. - -%% Initialize a 2DoF APA with Ga=Gs=1 -n_hexapod.actuator = initializeAPA(... - 'type', '2dof', ... - 'ga', 1, ... - 'gs', 1); - -% Identification without actuator or sensor constants -% The transfer function from $V_a$ to $V_s$, $d_e$ and $d_a$ is identified. - -%% 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, '/de'], 1, 'openoutput'); io_i = io_i + 1; % Encoder -io(io_i) = linio([mdl, '/da'], 1, 'openoutput'); io_i = io_i + 1; % Attocube - -%% Identification -Gs = linearize(mdl, io, 0.0, options); -Gs.InputName = {'Va'}; -Gs.OutputName = {'Vs', 'de', 'da'}; - -% Actuator Constant -% Then, the actuator constant can be computed as shown in Eq. eqref:eq:actuator_constant_formula by dividing the measured DC gain of the transfer function from $V_a$ to $d_e$ by the estimated DC gain of the transfer function from $V_a$ (in truth the actuator force called $F_a$) to $d_e$ using the Simscape model. - -%% Estimated Actuator Constant -ga = -mean(abs(enc_frf(f>10 & f<20)))./dcgain(Gs('de', 'Va')); % [N/V] - -% Sensor Constant -% Similarly, the sensor constant can be estimated using Eq. eqref:eq:sensor_constant_formula. - - -%% Estimated Sensor Constant -gs = -mean(abs(iff_frf(f>400 & f<500)))./(ga*abs(squeeze(freqresp(Gs('Vs', 'Va'), 1e3, 'Hz')))); % [V/m] - -% Comparison -% Let's now initialize the APA with identified sensor and actuator constant: - -%% Set the identified constants -n_hexapod.actuator = initializeAPA(... - 'type', '2dof', ... - 'ga', ga, ... % Actuator gain [N/V] - 'gs', gs); % Sensor gain [V/m] - - - -% And identify the dynamics with included constants. - -%% Identify again the dynamics with correct Ga,Gs -Gs = linearize(mdl, io, 0.0, options); -Gs = Gs*exp(-Ts*s); -Gs.InputName = {'Va'}; -Gs.OutputName = {'Vs', 'de', 'da'}; - - - -% The transfer functions from $V_a$ to $d_e$ are compared in Figure ref:fig:apa_act_constant_comp and the one from $V_a$ to $V_s$ are compared in Figure ref:fig:apa_sens_constant_comp. - - -%% Bode plot of the transfer function from u to de -freqs = logspace(1,4,1000); - -figure; -tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); - -ax1 = nexttile([2,1]); -hold on; -for i = 1:length(apa_nums) - plot(f, abs(enc_frf(:, i)), 'color', [0,0,0,0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, abs(squeeze(freqresp(Gs('de', 'Va'), freqs, 'Hz')))) -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -ylabel('Amplitude $d\mathcal{L}_m/u$ [m/V]'); set(gca, 'XTickLabel',[]); -hold off; -ylim([1e-8, 1e-3]); - -ax2 = nexttile; -hold on; -for i = 1:length(apa_nums) - plot(f, 180/pi*angle(enc_frf(:,1)), 'color', [0,0,0,0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, 180/pi*angle(squeeze(freqresp(Gs('de', 'Va'), freqs, 'Hz')))) -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); -hold off; -yticks(-360:90:360); ylim([-180, 180]); - -linkaxes([ax1,ax2],'x'); -xlim([10, 2e3]); - - - -% #+name: fig:apa_act_constant_comp -% #+caption: Comparison of the experimental data and Simscape model ($V_a$ to $d_e$) -% #+RESULTS: -% [[file:figs/apa_act_constant_comp.png]] - - -%% Bode plot of the transfer function from Va to Vs (both Simscape and measured FRF) -freqs = logspace(1,4,1000); - -figure; -tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); - -ax1 = nexttile([2,1]); -hold on; -for i = 1:length(apa_nums) - plot(f, abs(iff_frf(:, i)), 'color', [0,0,0,0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, abs(squeeze(freqresp(Gs('Vs', 'Va'), freqs, 'Hz')))) -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -ylabel('Amplitude $\tau_m/u$ [V/V]'); set(gca, 'XTickLabel',[]); -hold off; -ylim([1e-2, 1e2]); - -ax2 = nexttile; -hold on; -for i = 1:length(apa_nums) - plot(f, 180/pi*angle(iff_frf(:,1)), 'color', [0,0,0,0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, 180/pi*angle(squeeze(freqresp(Gs('Vs', 'Va'), freqs, 'Hz')))) -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); -hold off; -yticks(-360:90:360); ylim([-180, 180]); - -linkaxes([ax1,ax2],'x'); -xlim([10, 2e3]); - - - -% #+name: fig:apa_sens_constant_comp -% #+caption: Comparison of the experimental data and Simscape model ($V_a$ to $V_s$) -% #+RESULTS: -% [[file:figs/apa_sens_constant_comp.png]] - - -%% Compare the FRF and identified dynamics from Va to Vs and da -colors = colororder; - -figure; -tiledlayout(3, 2, 'TileSpacing', 'Compact', 'Padding', 'None'); - -ax1 = nexttile([2,1]); -hold on; -plot(f, abs(enc_frf(:, 1)), 'color', [0,0,0,0.2], ... - 'DisplayName', 'FRF'); -for i = 2:length(apa_nums) - plot(f, abs(enc_frf(:, i)), 'color', [0,0,0, 0.2], ... - 'HandleVisibility', 'off'); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, abs(squeeze(freqresp(Gs('da', 'Va'), freqs, 'Hz'))), '--', ... - 'DisplayName', 'Model') -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -ylabel('Amplitude $d_a/V_a$ [m/V]'); set(gca, 'XTickLabel',[]); -hold off; -ylim([1e-8, 1e-3]); -legend('location', 'southwest'); - -ax1b = nexttile([2,1]); -hold on; -plot(f, abs(iff_frf(:, i)), 'color', [0,0,0,0.2], ... - 'DisplayName', 'FRF'); -for i = 1:length(apa_nums) - plot(f, abs(iff_frf(:, i)), 'color', [0,0,0,0.2], ... - 'HandleVisibility', 'off'); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, abs(squeeze(freqresp(Gs('Vs', 'Va'), freqs, 'Hz'))), '--', ... - 'DisplayName', 'Model') -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]); -hold off; -ylim([1e-2, 1e2]); -legend('location', 'southeast'); - -ax2 = nexttile; -hold on; -for i = 1:length(apa_nums) - plot(f, 180/pi*angle(enc_frf(:, i)), 'color', [0,0,0, 0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, 180/pi*angle(squeeze(freqresp(Gs('da', 'Va'), freqs, 'Hz'))), '--') -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); -hold off; -yticks(-360:90:360); ylim([-180, 180]); - -ax2b = nexttile; -hold on; -for i = 1:length(apa_nums) - plot(f, 180/pi*angle(iff_frf(:, i)), 'color', [0,0,0, 0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, 180/pi*angle(squeeze(freqresp(Gs('Vs', 'Va'), freqs, 'Hz'))), '--') -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); -hold off; -yticks(-360:90:360); ylim([-180, 180]); - -linkaxes([ax1,ax2,ax1b,ax2b],'x'); -xlim([10, 2e3]); - -% Flexible APA -% The Simscape APA model is initialized as a flexible one with unity "constants". - -%% Initialize the APA as a flexible body -n_hexapod.actuator = initializeAPA(... - 'type', 'flexible', ... - 'ga', 1, ... - 'gs', 1); - -% Identification without actuator or sensor constants -% The dynamics from $V_a$ to $V_s$, $d_e$ and $d_a$ is identified. - -%% Identify the dynamics -Gs = linearize(mdl, io, 0.0, options); -Gs.InputName = {'Va'}; -Gs.OutputName = {'Vs', 'de', 'da'}; - -% Actuator Constant -% Then, the actuator constant can be computed as shown in Eq. eqref:eq:actuator_constant_formula: - -%% Actuator Constant -ga = -mean(abs(enc_frf(f>10 & f<20)))./dcgain(Gs('de', 'Va')); % [N/V] - -% Sensor Constant - -%% Sensor Constant -gs = -mean(abs(iff_frf(f>400 & f<500)))./(ga*abs(squeeze(freqresp(Gs('Vs', 'Va'), 1e3, 'Hz')))); % [V/m] - -% Comparison -% Let's now initialize the flexible APA with identified sensor and actuator constant: - -%% Set the identified constants -n_hexapod.actuator = initializeAPA(... - 'type', 'flexible', ... - 'ga', ga, ... % Actuator gain [N/V] - 'gs', gs); % Sensor gain [V/m] - - - -% And identify the dynamics with included constants. - -%% Identify with updated constants -Gs = linearize(mdl, io, 0.0, options); -Gs = Gs*exp(-Ts*s); -Gs.InputName = {'Va'}; -Gs.OutputName = {'Vs', 'de', 'da'}; - - - -% The obtained dynamics is compared with the measured one in Figures ref:fig:apa_act_constant_comp_flex and ref:fig:apa_sens_constant_comp_flex. - - -%% Bode plot of the transfer function from V_a to d_e (both Simscape and measured FRF) -figure; -tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); - -ax1 = nexttile([2,1]); -hold on; -for i = 1:length(apa_nums) - plot(f, abs(enc_frf(:, i)), 'color', [0,0,0,0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, abs(squeeze(freqresp(Gs('de', 'Va'), freqs, 'Hz')))) -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -ylabel('Amplitude $d\mathcal{L}_m/u$ [m/V]'); set(gca, 'XTickLabel',[]); -hold off; -ylim([1e-9, 1e-3]); - -ax2 = nexttile; -hold on; -for i = 1:length(apa_nums) - plot(f, 180/pi*angle(enc_frf(:,1)), 'color', [0,0,0,0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, 180/pi*angle(squeeze(freqresp(Gs('de', 'Va'), freqs, 'Hz')))) -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); -hold off; -yticks(-360:90:360); ylim([-180, 180]); - -linkaxes([ax1,ax2],'x'); -xlim([10, 2e3]); - - - -% #+name: fig:apa_act_constant_comp_flex -% #+caption: Comparison of the experimental data and Simscape model ($u$ to $d\mathcal{L}_m$) -% #+RESULTS: -% [[file:figs/apa_act_constant_comp_flex.png]] - - -%% Bode plot of the transfer function from Va to Vs (both Simscape and measured FRF) -figure; -tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None'); - -ax1 = nexttile([2,1]); -hold on; -for i = 1:length(apa_nums) - plot(f, abs(iff_frf(:, i)), 'color', [0,0,0,0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, abs(squeeze(freqresp(Gs('Vs', 'Va'), freqs, 'Hz')))) -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -ylabel('Amplitude $\tau_m/u$ [V/V]'); set(gca, 'XTickLabel',[]); -hold off; -ylim([1e-2, 1e2]); - -ax2 = nexttile; -hold on; -for i = 1:length(apa_nums) - plot(f, 180/pi*angle(iff_frf(:,1)), 'color', [0,0,0,0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, 180/pi*angle(squeeze(freqresp(Gs('Vs', 'Va'), freqs, 'Hz')))) -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); -hold off; -yticks(-360:90:360); ylim([-180, 180]); - -linkaxes([ax1,ax2],'x'); -xlim([10, 2e3]); - -% Optimize 2-DoF model to fit the experimental Data -% <> -% The parameters of the 2DoF model presented in Section ref:sec:apa_2dof_model are now optimize such that the model best matches the measured FRF. - -% After optimization, the following parameters are used: - -%% Optimized parameters -n_hexapod.actuator = initializeAPA('type', '2dof', ... - 'Ga', -32.2, ... - 'Gs', 0.088, ... - 'k', ones(6,1)*0.38e6, ... - 'ke', ones(6,1)*1.75e6, ... - 'ka', ones(6,1)*3e7, ... - 'c', ones(6,1)*1.3e2, ... - 'ce', ones(6,1)*1e1, ... - 'ca', ones(6,1)*1e1 ... - ); - -%% 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, '/de'], 1, 'openoutput'); io_i = io_i + 1; % Encoder - -%% Identification with optimized parameters -Gs = exp(-s*Ts)*linearize(mdl, io, 0.0, options); -Gs.InputName = {'Va'}; -Gs.OutputName = {'Vs', 'de'}; - - - -% The dynamics is identified using the Simscape model and compared with the measured FRF in Figure ref:fig:comp_apa_plant_after_opt. - -%% Comparison of the experimental data and Simscape Model -freqs = 5*logspace(0, 3, 1000); -figure; -tiledlayout(3, 2, 'TileSpacing', 'Compact', 'Padding', 'None'); - -ax1 = nexttile([2,1]); -hold on; -for i = 1:length(apa_nums) - plot(f, abs(enc_frf(:, i)), 'color', [0,0,0,0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, abs(squeeze(freqresp(Gs('de', 'Va'), freqs, 'Hz')))) -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -ylabel('Amplitude $d_e/V_a$ [m/V]'); set(gca, 'XTickLabel',[]); -hold off; -ylim([1e-8, 1e-3]); - -ax1b = nexttile([2,1]); -hold on; -for i = 1:length(apa_nums) - plot(f, abs(iff_frf(:, i)), 'color', [0,0,0,0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, abs(squeeze(freqresp(Gs('Vs', 'Va'), freqs, 'Hz')))) -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); -ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]); -hold off; -ylim([1e-2, 1e2]); - -ax2 = nexttile; -hold on; -for i = 1:length(apa_nums) - plot(f, 180/pi*angle(enc_frf(:, i)), 'color', [0,0,0,0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, 180/pi*angle(squeeze(freqresp(Gs('de', 'Va'), freqs, 'Hz')))) -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); -hold off; -yticks(-360:90:360); ylim([-180, 180]); - -ax2b = nexttile; -hold on; -for i = 1:length(apa_nums) - plot(f, 180/pi*angle(iff_frf(:, i)), 'color', [0,0,0,0.2]); -end -set(gca,'ColorOrderIndex',1); -plot(freqs, 180/pi*angle(squeeze(freqresp(Gs('Vs', 'Va'), freqs, 'Hz')))) -hold off; -set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin'); -xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); -hold off; -yticks(-360:90:360); ylim([-180, 180]); - -linkaxes([ax1,ax2,ax1b,ax2b],'x'); -xlim([10, 2e3]); diff --git a/matlab/test_apa_4_model_flexible.m b/matlab/test_apa_4_model_flexible.m index b07cf48..1d538f5 100644 --- a/matlab/test_apa_4_model_flexible.m +++ b/matlab/test_apa_4_model_flexible.m @@ -116,6 +116,7 @@ ka = cE*A/L; % Stiffness of the two stacks [N/m] ga = d33*n*ka; % Actuator Constant [N/V] % Comparison of the obtained dynamics +% <> % The obtained dynamics using the /super element/ with the tuned "sensor gain" and "actuator gain" are compared with the experimentally identified frequency response functions in Figure ref:fig:test_apa_super_element_comp_frf. diff --git a/test-bench-apa.org b/test-bench-apa.org index 8311dc0..43ac64a 100644 --- a/test-bench-apa.org +++ b/test-bench-apa.org @@ -699,9 +699,6 @@ add_mass_cc.de = add_mass_cc.de - mean(add_mass_cc.de(add_mass_cc.t<11)); %% Estimation of the stiffness in Open Circuit and Closed-Circuit apa_k_oc = 9.8 * added_mass / (mean(add_mass_oc.de(add_mass_oc.t > 12 & add_mass_oc.t < 12.5)) - mean(add_mass_oc.de(add_mass_oc.t > 20 & add_mass_oc.t < 20.5))); apa_k_sc = 9.8 * added_mass / (mean(add_mass_cc.de(add_mass_cc.t > 12 & add_mass_cc.t < 12.5)) - mean(add_mass_cc.de(add_mass_cc.t > 20 & add_mass_cc.t < 20.5))); - -%% Estimated coupling factor -sqrt(1 - apa_k_sc/apa_k_oc) #+end_src ** Dynamics @@ -1051,7 +1048,7 @@ Noverlap = floor(Nfft/2); [~, f] = tfestimate(data.data(1).id_plant(1:end), data.data(1).dL(1:end), win, Noverlap, Nfft, 1/Ts); %% Gains used for analysis are between 1 and 50 -i_kept = [5:10] +i_kept = [5:10]; %% Identify the damped plant from u' to de for different IFF gains G_dL_frf = {zeros(1,length(i_kept))}; @@ -1080,7 +1077,6 @@ 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) - Niter = 100; % Number of iteration. N = 2; % Order of approximation poles = [-25 - 1i*60, -25 + 1i*60]; % First get for the pole location @@ -1088,7 +1084,7 @@ poles = [-25 - 1i*60, -25 + 1i*60]; % First get for the pole location G_dL_id = {zeros(1,length(i_kept))}; % Identification just between two frequencies -f_keep = (f>20 & f<200); +f_keep = (f>5 & f<200); for i = 1:length(i_kept) %% Estimate resonance frequency and damping @@ -1135,12 +1131,11 @@ The two obtained root loci are compared in Figure ref:fig:test_apa_iff_root_locu #+begin_src matlab :exports none :results none %% Root Locus of the APA300ML with Integral Force Feedback % Comparison between the computed root locus from the plant model and the root locus estimated from the damped plant pole identification -figure; gains = logspace(-1, 3, 1000); figure; hold on; -G_iff_poles = pole(G_iff_model); +G_iff_poles = pole(pade(G_iff_model)); i = imag(G_iff_poles) > 100; % Only keep relevant poles plot(real(G_iff_poles(i)), imag(G_iff_poles(i)), 'kx', ... 'DisplayName', '$g = 0$'); @@ -1150,7 +1145,7 @@ plot(real(G_iff_zeros(i)), imag(G_iff_zeros(i)), 'ko', ... 'HandleVisibility', 'off'); for g = gains - clpoles = pole(feedback(G_iff_model, g*K_iff, 1)); + clpoles = pole(feedback(pade(G_iff_model), g*K_iff, 1)); i = imag(clpoles) > 100; % Only keep relevant poles plot(real(clpoles(i)), imag(clpoles(i)), 'k.', ... 'HandleVisibility', 'off'); @@ -1318,7 +1313,7 @@ ga = -abs(enc_frf(i_f,1))./abs(evalfr(G_norm('de', 'u'), 1i*2*pi*fa)); % Estimation ot the Sensor Gain fs = 600; % Frequency where the two FRF should match [Hz] -[~, i_f] = min(abs(f - fs)) +[~, i_f] = min(abs(f - fs)); gs = -abs(iff_frf(i_f,1))./abs(evalfr(G_norm('Vs', 'u'), 1i*2*pi*fs))/ga; #+end_src