tdehaeze/test-bench-apa
1
0
Fork 0
Code Issues Pull Requests Releases Activity
5734ae5d6b
test-bench-apa/index.org

765 lines
22 KiB
Org Mode
Raw Blame History

#+TITLE: Test Bench APA95ML
:DRAWER:
#+STARTUP: overview
#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas
#+HTML_LINK_HOME: ../index.html
#+HTML_LINK_UP: ../index.html
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="https://research.tdehaeze.xyz/css/style.css"/>
#+HTML_HEAD: <script type="text/javascript" src="https://research.tdehaeze.xyz/js/script.js"></script>
#+HTML_MATHJAX: align: center tagside: right font: TeX
#+PROPERTY: header-args:matlab :session *MATLAB*
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :results none
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :output-dir figs
#+PROPERTY: header-args:matlab+ :tangle no
#+PROPERTY: header-args:matlab+ :mkdirp yes
#+PROPERTY: header-args:shell :eval no-export
:END:
* Introduction :ignore:
#+name: fig:setup_picture
#+caption: Picture of the Setup
[[file:figs/setup_picture.png]]
#+name: fig:setup_zoom
#+caption: Zoom on the APA
[[file:figs/setup_zoom.png]]
* Huddle Test
:PROPERTIES:
:header-args:matlab+: :tangle matlab/huddle_test.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
<<sec:huddle_test>>
** Introduction :ignore:
** 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/mat/');
#+end_src
#+begin_src matlab :eval no
addpath('./mat/');
#+end_src
** Load Data :noexport:
#+begin_src matlab
load('huddle_test.mat', 't', 'y');
#+end_src
#+begin_src matlab
y = y - mean(y(1000:end));
#+end_src
** Time Domain Data
#+begin_src matlab :exports none
figure;
plot(t(1000:end), y(1000:end))
ylabel('Output Displacement [m]'); xlabel('Time [s]');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/huddle_test_time_domain.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:huddle_test_time_domain
#+caption: Measurement of the Mass displacement during Huddle Test
#+RESULTS:
[[file:figs/huddle_test_time_domain.png]]
** PSD of Measurement Noise
#+begin_src matlab
Ts = t(end)/(length(t)-1);
Fs = 1/Ts;
win = hanning(ceil(1*Fs));
#+end_src
#+begin_src matlab
[pxx, f] = pwelch(y(1000:end), win, [], [], Fs);
#+end_src
#+begin_src matlab :exports none
figure;
plot(f, sqrt(pxx));
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
xlim([1, Fs/2]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/huddle_test_pdf.pdf', 'width', 'wide', 'height', 'tall');
#+end_src
#+name: fig:huddle_test_pdf
#+caption: Amplitude Spectral Density of the Displacement during Huddle Test
#+RESULTS:
[[file:figs/huddle_test_pdf.png]]
* Identification of the dynamics from actuator to displacement
:PROPERTIES:
:header-args:matlab+: :tangle matlab/motion_identification.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
<<sec:motion_identification>>
** Introduction :ignore:
** 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/mat/');
#+end_src
#+begin_src matlab :eval no
addpath('./mat/');
#+end_src
** Load Data
#+begin_src matlab
ht = load('huddle_test.mat', 't', 'u', 'y');
load('apa95ml_5kg_Amp_E505.mat', 't', 'u', 'um', 'y');
#+end_src
#+begin_src matlab
u = 10*(u - mean(u)); % Input Voltage of Piezo [V]
um = 10*(um - mean(um)); % Monitor [V]
y = y - mean(y); % Mass displacement [m]
ht.u = 10*(ht.u - mean(ht.u));
ht.y = ht.y - mean(ht.y);
#+end_src
** Comparison of the PSD with Huddle Test
#+begin_src matlab
Ts = t(end)/(length(t)-1);
Fs = 1/Ts;
win = hanning(ceil(1*Fs));
#+end_src
#+begin_src matlab
[pxx, f] = pwelch(y, win, [], [], Fs);
[pht, ~] = pwelch(ht.y, win, [], [], Fs);
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
plot(f, sqrt(pxx), 'DisplayName', '5kg');
plot(f, sqrt(pht), 'DisplayName', 'Huddle Test');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
legend('location', 'southwest');
xlim([1, Fs/2]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/apa95ml_5kg_PI_pdf_comp_huddle.pdf', 'width', 'wide', 'height', 'tall');
#+end_src
#+name: fig:apa95ml_5kg_PI_pdf_comp_huddle
#+caption: Comparison of the ASD for the identification test and the huddle test
#+RESULTS:
[[file:figs/apa95ml_5kg_PI_pdf_comp_huddle.png]]
** Compute TF estimate and Coherence
#+begin_src matlab
Ts = t(end)/(length(t)-1);
Fs = 1/Ts;
#+end_src
#+begin_src matlab
win = hann(ceil(1/Ts));
[tf_est, f] = tfestimate(u, -y, win, [], [], 1/Ts);
[tf_um , ~] = tfestimate(um, -y, win, [], [], 1/Ts);
[co_est, ~] = mscohere( um, -y, win, [], [], 1/Ts);
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
plot(f, co_est, 'k-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Coherence'); xlabel('Frequency [Hz]');
hold off;
xlim([10, 5e3]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/apa95ml_5kg_PI_coh.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:apa95ml_5kg_PI_coh
#+caption: Coherence
#+RESULTS:
[[file:figs/apa95ml_5kg_PI_coh.png]]
#+begin_src matlab :exports none
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(tf_est), 'DisplayName', 'Input Voltage')
plot(f, abs(tf_um), 'DisplayName', 'Monitor Voltage')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude [m/V]'); xlabel('Frequency [Hz]');
hold off;
legend('location', 'southwest')
ax2 = subplot(2, 1, 2);
hold on;
plot(f, 180/pi*unwrap(angle(tf_est)))
plot(f, 180/pi*unwrap(angle(tf_um)))
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
ylim([-540, 0]);
yticks(-540:90:0);
linkaxes([ax1,ax2], 'x');
xlim([10, 5e3]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/apa95ml_5kg_PI_tf.pdf', 'width', 'full', 'height', 'full');
#+end_src
#+name: fig:apa95ml_5kg_PI_tf
#+caption: Estimation of the transfer function from input voltage to displacement
#+RESULTS:
[[file:figs/apa95ml_5kg_PI_tf.png]]
** Comparison with the FEM model
#+begin_src matlab
load('fem_model_5kg.mat', 'G');
#+end_src
#+begin_src matlab :exports none
freqs = logspace(0, 4, 1000);
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(tf_um), 'DisplayName', 'Identification')
plot(freqs, abs(squeeze(freqresp(G, freqs, 'Hz'))), 'DisplayName', 'FEM')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude [m/V]'); xlabel('Frequency [Hz]');
legend('location', 'northeast')
hold off;
ax2 = subplot(2, 1, 2);
hold on;
plot(f, 180/pi*unwrap(angle(tf_um)))
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G, freqs, 'Hz')))))
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
ylim([-540, 0]);
yticks(-540:90:0);
linkaxes([ax1,ax2], 'x');
xlim([10, 5e3]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/apa95ml_5kg_pi_comp_fem.pdf', 'width', 'full', 'height', 'full');
#+end_src
#+name: fig:apa95ml_5kg_pi_comp_fem
#+caption: Comparison of the identified transfer function and the one estimated from the FE model
#+RESULTS:
[[file:figs/apa95ml_5kg_pi_comp_fem.png]]
* Identification of the dynamics from actuator to force sensor
:PROPERTIES:
:header-args:matlab+: :tangle matlab/force_sensor_identification.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
<<sec:force_sensor_identification>>
** Introduction :ignore:
Two measurements are performed:
- Speedgoat DAC => Voltage Amplifier (x20) => 1 Piezo Stack => ... => 2 Stacks as Force Sensor (parallel) => Speedgoat ADC
- Speedgoat DAC => Voltage Amplifier (x20) => 2 Piezo Stacks (parallel) => ... => 1 Stack as Force Sensor => Speedgoat ADC
The obtained dynamics from force actuator to force sensor are compare with the FEM model.
** 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/mat/');
#+end_src
#+begin_src matlab :eval no
addpath('./mat/');
#+end_src
** Load Data :ignore:
The data are loaded:
#+begin_src matlab
a_ss = load('apa95ml_5kg_1a_2s.mat', 't', 'u', 'y', 'v');
aa_s = load('apa95ml_5kg_2a_1s.mat', 't', 'u', 'y', 'v');
load('G_force_sensor_5kg.mat', 'G');
#+end_src
** Adjust gain :ignore:
Let's use the amplifier gain to obtain the true voltage applied to the actuator stack(s)
The parameters of the piezoelectric stacks are defined below:
#+begin_src matlab
d33 = 3e-10; % Strain constant [m/V]
n = 80; % Number of layers per stack
eT = 1.6e-8; % Permittivity under constant stress [F/m]
sD = 2e-11; % Elastic compliance under constant electric displacement [m2/N]
ka = 235e6; % Stack stiffness [N/m]
#+end_src
From the FEM, we construct the transfer function from DAC voltage to ADC voltage.
#+begin_src matlab
Gfem_aa_s = exp(-s/1e4)*20*(2*d33*n*ka)*(G(3,1)+G(3,2))*d33/(eT*sD*n);
Gfem_a_ss = exp(-s/1e4)*20*( d33*n*ka)*(G(3,1)+G(2,1))*d33/(eT*sD*n);
#+end_src
** Compute TF estimate and Coherence :ignore:
The transfer function from input voltage to output voltage are computed and shown in Figure [[fig:bode_plot_force_sensor_voltage_comp_fem]].
#+begin_src matlab
Ts = a_ss.t(end)/(length(a_ss.t)-1);
Fs = 1/Ts;
win = hann(ceil(10/Ts));
[tf_a_ss, f] = tfestimate(a_ss.u, a_ss.v, win, [], [], 1/Ts);
[coh_a_ss, ~] = mscohere( a_ss.u, a_ss.v, win, [], [], 1/Ts);
[tf_aa_s, f] = tfestimate(aa_s.u, aa_s.v, win, [], [], 1/Ts);
[coh_aa_s, ~] = mscohere( aa_s.u, aa_s.v, win, [], [], 1/Ts);
#+end_src
#+begin_src matlab :exports none
freqs = logspace(1, 4, 1000);
figure;
ax1 = subplot(2, 1, 1);
hold on;
set(gca,'ColorOrderIndex',1)
plot(f, abs(tf_aa_s), '-')
set(gca,'ColorOrderIndex',1)
plot(freqs, abs(squeeze(freqresp(Gfem_aa_s, freqs, 'Hz'))), '--')
set(gca,'ColorOrderIndex',2)
plot(f, abs(tf_a_ss), '-')
set(gca,'ColorOrderIndex',2)
plot(freqs, abs(squeeze(freqresp(Gfem_a_ss, freqs, 'Hz'))), '--')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); xlabel('Frequency [Hz]');
hold off;
ylim([1e-2, 1e2]);
ax2 = subplot(2, 1, 2);
hold on;
set(gca,'ColorOrderIndex',1)
plot(f, 180/pi*angle(tf_aa_s), '-', 'DisplayName', '2 Act - 1 Sen')
set(gca,'ColorOrderIndex',1)
plot(freqs, 180/pi*angle(squeeze(freqresp(Gfem_aa_s, freqs, 'Hz'))), '--', 'DisplayName', '2 Act - 1 Sen, - FEM')
set(gca,'ColorOrderIndex',2)
plot(f, 180/pi*angle(tf_a_ss), '-', 'DisplayName', '1 Act - 2 Sen')
set(gca,'ColorOrderIndex',2)
plot(freqs, 180/pi*angle(squeeze(freqresp(Gfem_a_ss, freqs, 'Hz'))), '--', 'DisplayName', '1 Act - 2 Sen, - FEM')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
ylim([-180, 180]);
yticks(-180:90:180);
legend('location', 'northeast')
linkaxes([ax1,ax2], 'x');
xlim([10, 5e3]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/bode_plot_force_sensor_voltage_comp_fem.pdf', 'width', 'full', 'height', 'full');
#+end_src
#+name: fig:bode_plot_force_sensor_voltage_comp_fem
#+caption: Comparison of the identified dynamics from voltage output to voltage input and the FEM
#+RESULTS:
[[file:figs/bode_plot_force_sensor_voltage_comp_fem.png]]
** System Identification
#+begin_src matlab
w_z = 2*pi*111; % Zeros frequency [rad/s]
w_p = 2*pi*255; % Pole frequency [rad/s]
xi_z = 0.05;
xi_p = 0.015;
G_inf = 2;
Gi = G_inf*(s^2 - 2*xi_z*w_z*s + w_z^2)/(s^2 + 2*xi_p*w_p*s + w_p^2);
#+end_src
#+begin_src matlab :exports none
freqs = logspace(1, 4, 1000);
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(tf_aa_s), '-')
plot(freqs, abs(squeeze(freqresp(Gi, freqs, 'Hz'))), '--')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); xlabel('Frequency [Hz]');
hold off;
ylim([1e-2, 1e2]);
ax2 = subplot(2, 1, 2);
hold on;
plot(f, 180/pi*angle(tf_aa_s), '-', 'DisplayName', '2 Act - 1 Sen')
plot(freqs, 180/pi*angle(squeeze(freqresp(Gi, freqs, 'Hz'))), '--', 'DisplayName', '2 Act - 1 Sen, - FEM')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
ylim([-180, 180]);
yticks(-180:90:180);
legend('location', 'northeast')
linkaxes([ax1,ax2], 'x');
xlim([10, 5e3]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/iff_plant_identification_apa95ml.pdf', 'width', 'full', 'height', 'full');
#+end_src
#+name: fig:iff_plant_identification_apa95ml
#+caption: Identification of the IFF plant
#+RESULTS:
[[file:figs/iff_plant_identification_apa95ml.png]]
** Integral Force Feedback
#+begin_src matlab :exports none
gains = logspace(0, 5, 1000);
figure;
hold on;
plot(real(pole(Gi)), imag(pole(Gi)), 'kx');
plot(real(tzero(Gi)), imag(tzero(Gi)), 'ko');
for i = 1:length(gains)
cl_poles = pole(feedback(Gi, (gains(i)/(s + 2*2*pi)*s/(s + 0.5*2*pi))));
plot(real(cl_poles), imag(cl_poles), 'k.');
end
ylim([0, 1800]);
xlim([-1600,200]);
xlabel('Real Part')
ylabel('Imaginary Part')
axis square
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/root_locus_iff_apa95ml_identification.pdf', 'width', 'wide', 'height', 'tall');
#+end_src
#+name: fig:root_locus_iff_apa95ml_identification
#+caption: Root Locus for IFF
#+RESULTS:
[[file:figs/root_locus_iff_apa95ml_identification.png]]
* Integral Force Feedback
:PROPERTIES:
:header-args:matlab+: :tangle matlab/integral_force_feedback.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
<<sec:integral_force_feedback>>
** Introduction :ignore:
** 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/mat/');
#+end_src
#+begin_src matlab :eval no
addpath('./mat/');
#+end_src
** First tests with few gains
#+begin_src matlab
iff_g10 = load('apa95ml_iff_g10_res.mat', 'u', 't', 'y', 'v');
iff_g100 = load('apa95ml_iff_g100_res.mat', 'u', 't', 'y', 'v');
iff_of = load('apa95ml_iff_off_res.mat', 'u', 't', 'y', 'v');
#+end_src
#+begin_src matlab
Ts = 1e-4;
win = hann(ceil(10/Ts));
[tf_iff_g10, f] = tfestimate(iff_g10.u, iff_g10.y, win, [], [], 1/Ts);
[co_iff_g10, ~] = mscohere(iff_g10.u, iff_g10.y, win, [], [], 1/Ts);
[tf_iff_g100, f] = tfestimate(iff_g100.u, iff_g100.y, win, [], [], 1/Ts);
[co_iff_g100, ~] = mscohere(iff_g100.u, iff_g100.y, win, [], [], 1/Ts);
[tf_iff_of, ~] = tfestimate(iff_of.u, iff_of.y, win, [], [], 1/Ts);
[co_iff_of, ~] = mscohere(iff_of.u, iff_of.y, win, [], [], 1/Ts);
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
plot(f, co_iff_of, '-', 'DisplayName', 'g=0')
plot(f, co_iff_g10, '-', 'DisplayName', 'g=10')
plot(f, co_iff_g100, '-', 'DisplayName', 'g=100')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Coherence'); xlabel('Frequency [Hz]');
hold off;
legend();
xlim([60, 600])
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/iff_first_test_coherence.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:iff_first_test_coherence
#+caption: Coherence
#+RESULTS:
[[file:figs/iff_first_test_coherence.png]]
#+begin_src matlab :exports none
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(tf_iff_of), '-', 'DisplayName', 'g=0')
plot(f, abs(tf_iff_g10), '-', 'DisplayName', 'g=10')
plot(f, abs(tf_iff_g100), '-', 'DisplayName', 'g=100')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); xlabel('Frequency [Hz]');
hold off;
legend();
ax2 = subplot(2, 1, 2);
hold on;
plot(f, 180/pi*angle(-tf_iff_of), '-')
plot(f, 180/pi*angle(-tf_iff_g10), '-')
plot(f, 180/pi*angle(-tf_iff_g100), '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
xlim([60, 600]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/iff_first_test_bode_plot.pdf', 'width', 'full', 'height', 'full');
#+end_src
#+name: fig:iff_first_test_bode_plot
#+caption: Bode plot for different values of IFF gain
#+RESULTS:
[[file:figs/iff_first_test_bode_plot.png]]
** Second test with many Gains
#+begin_src matlab
load('apa95ml_iff_test.mat', 'results');
#+end_src
#+begin_src matlab
Ts = 1e-4;
win = hann(ceil(10/Ts));
#+end_src
#+begin_src matlab
tf_iff = {zeros(1, length(results))};
co_iff = {zeros(1, length(results))};
g_iff = [0, 1, 5, 10, 50, 100];
for i=1:length(results)
[tf_est, f] = tfestimate(results{i}.u, results{i}.y, win, [], [], 1/Ts);
[co_est, ~] = mscohere(results{i}.u, results{i}.y, win, [], [], 1/Ts);
tf_iff(i) = {tf_est};
co_iff(i) = {co_est};
end
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
for i = 1:length(results)
plot(f, co_iff{i}, '-', 'DisplayName', sprintf('g = %0.f', g_iff(i)))
end
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Coherence'); xlabel('Frequency [Hz]');
hold off;
legend();
xlim([60, 600])
#+end_src
#+begin_src matlab :exports none
figure;
ax1 = subplot(2, 1, 1);
hold on;
for i = 1:length(results)
plot(f, abs(tf_iff{i}), '-', 'DisplayName', sprintf('g = %0.f', g_iff(i)))
end
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); xlabel('Frequency [Hz]');
hold off;
legend();
ax2 = subplot(2, 1, 2);
hold on;
for i = 1:length(results)
plot(f, 180/pi*angle(-tf_iff{i}), '-')
end
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
xlim([60, 600]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/iff_results_bode_plots.pdf', 'width', 'full', 'height', 'full');
#+end_src
#+name: fig:iff_results_bode_plots
#+caption:
#+RESULTS:
[[file:figs/iff_results_bode_plots.png]]
#+begin_src matlab
G_id = {zeros(1,length(results))};
f_start = 70; % [Hz]
f_end = 500; % [Hz]
for i = 1:length(results)
tf_id = tf_iff{i}(sum(f<f_start):length(f)-sum(f>f_end));
f_id = f(sum(f<f_start):length(f)-sum(f>f_end));
gfr = idfrd(tf_id, 2*pi*f_id, Ts);
G_id(i) = {procest(gfr,'P2UDZ')};
end
#+end_src
#+begin_src matlab :exports none
figure;
ax1 = subplot(2, 1, 1);
hold on;
for i = 1:length(results)
set(gca,'ColorOrderIndex',i)
plot(f, abs(tf_iff{i}), '-', 'DisplayName', sprintf('g = %0.f', g_iff(i)))
set(gca,'ColorOrderIndex',i)
plot(f, abs(squeeze(freqresp(G_id{i}, f, 'Hz'))), '--', 'HandleVisibility', 'off')
end
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); xlabel('Frequency [Hz]');
hold off;
legend();
ax2 = subplot(2, 1, 2);
hold on;
for i = 1:length(results)
set(gca,'ColorOrderIndex',i)
plot(f, 180/pi*angle(tf_iff{i}), '-')
set(gca,'ColorOrderIndex',i)
plot(f, 180/pi*angle(squeeze(freqresp(G_id{i}, f, 'Hz'))), '--')
end
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
xlim([60, 600]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/iff_results_bode_plots_identification.pdf', 'width', 'full', 'height', 'full');
#+end_src
#+name: fig:iff_results_bode_plots_identification
#+caption:
#+RESULTS:
[[file:figs/iff_results_bode_plots_identification.png]]
#+begin_src matlab :exports none
w_z = 2*pi*111; % Zeros frequency [rad/s]
w_p = 2*pi*255; % Pole frequency [rad/s]
xi_z = 0.05;
xi_p = 0.015;
G_inf = 2;
Gi = G_inf*(s^2 - 2*xi_z*w_z*s + w_z^2)/(s^2 + 2*xi_p*w_p*s + w_p^2);
gains = logspace(0, 5, 1000);
figure;
hold on;
plot(real(pole(Gi)), imag(pole(Gi)), 'kx', 'HandleVisibility', 'off');
plot(real(tzero(Gi)), imag(tzero(Gi)), 'ko', 'HandleVisibility', 'off');
for i = 1:length(results)
set(gca,'ColorOrderIndex',i)
plot(real(pole(G_id{i})), imag(pole(G_id{i})), 'o', 'DisplayName', sprintf('g = %0.f', g_iff(i)));
end
for i = 1:length(gains)
cl_poles = pole(feedback(Gi, (gains(i)/(s + 2*pi*2))));
plot(real(cl_poles), imag(cl_poles), 'k.', 'HandleVisibility', 'off');
end
ylim([0, 1800]);
xlim([-1600,200]);
xlabel('Real Part')
ylabel('Imaginary Part')
axis square
legend('location', 'northwest');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/iff_results_root_locus.pdf', 'width', 'wide', 'height', 'tall');
#+end_src
#+name: fig:iff_results_root_locus
#+caption:
#+RESULTS:
[[file:figs/iff_results_root_locus.png]]
Reference in New Issue View Git Blame Copy Permalink