#+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="./css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="./css/readtheorg.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="./css/zenburn.css"/>
#+HTML_HEAD: <script type="text/javascript" src="./js/jquery.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="./js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="./js/jquery.stickytableheaders.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="./js/readtheorg.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

#+PROPERTY: header-args:latex  :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/Cloud/tikz/org/}{config.tex}")
#+PROPERTY: header-args:latex+ :imagemagick t :fit yes
#+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150
#+PROPERTY: header-args:latex+ :imoutoptions -quality 100
#+PROPERTY: header-args:latex+ :results raw replace :buffer no
#+PROPERTY: header-args:latex+ :eval no-export
#+PROPERTY: header-args:latex+ :exports both
#+PROPERTY: header-args:latex+ :mkdirp yes
#+PROPERTY: header-args:latex+ :output-dir figs
: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]]

* Setup
:PROPERTIES:
:header-args:matlab+: :tangle matlab/setup_experiment.m
:header-args:matlab+: :comments org :mkdirp yes
:END:

** Parameters
#+begin_src matlab
  Ts = 1e-4;
#+end_src

** Filter White Noise
#+begin_src matlab
  Glpf = 1/(1 + s/2/pi/500);

  Gz = c2d(Glpf, Ts, 'tustin');
#+end_src

* Run Experiment and Save Data
:PROPERTIES:
:header-args:matlab+: :tangle matlab/run_experiment.m
:header-args:matlab+: :comments org :mkdirp yes
:END:

** Load Data
#+begin_src matlab
  data = SimulinkRealTime.utils.getFileScopeData('data/apa95ml.dat').data;
#+end_src

** Save Data
#+begin_src matlab
  u = data(:, 1); % Input Voltage [V]
  y = data(:, 2); % Output Displacement [m]
  t = data(:, 3); % Time [s]
#+end_src

#+begin_src matlab
  save('./mat/huddle_test.mat', 't', 'u', 'y', 'Glpf');
#+end_src

* Huddle Test
:PROPERTIES:
:header-args:matlab+: :tangle matlab/huddle_test.m
:header-args:matlab+: :comments org :mkdirp yes
:END:

** 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

** Load Data                                                       :noexport:
#+begin_src matlab
  load('./mat/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]]

* Transfer Function Estimation using the DAC as the driver
:PROPERTIES:
:header-args:matlab+: :tangle matlab/tf_estimation.m
:header-args:matlab+: :comments org :mkdirp yes
:END:

** Introduction                                                      :ignore:
#+begin_important
Results presented in this sections are wrong as the ADC cannot deliver enought current to the piezoelectric actuator.
#+end_important

** 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

** Load Data                                                       :noexport:
#+begin_src matlab
  ht = load('./mat/huddle_test.mat', 't', 'u', 'y');
  load('./mat/apa95ml_5kg_10V.mat', 't', 'u', 'y');
#+end_src

** Time Domain Data
#+begin_src matlab :exports none
  figure;

  subplot(1,2,1);
  plot(t, u)
  ylabel('Input Voltage [V]'); xlabel('Time [s]');


  subplot(1,2,2);
  plot(t, y)
  ylabel('Output Displacement [m]'); xlabel('Time [s]');
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/apa95ml_5kg_10V_time_domain.pdf', 'width', 'full', 'height', 'tall');
#+end_src

#+name: fig:apa95ml_5kg_10V_time_domain
#+caption: Time domain signals during the test
#+RESULTS:
[[file:figs/apa95ml_5kg_10V_time_domain.png]]

** 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', 'northeast');
  xlim([1, Fs/2]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/apa95ml_5kg_10V_pdf_comp_huddle.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:apa95ml_5kg_10V_pdf_comp_huddle
#+caption: Comparison of the ASD for the identification test and the huddle test
#+RESULTS:
[[file:figs/apa95ml_5kg_10V_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);
  [co_est, ~] = mscohere(  u, -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_10V_coh.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:apa95ml_5kg_10V_coh
#+caption: Coherence
#+RESULTS:
[[file:figs/apa95ml_5kg_10V_coh.png]]

#+begin_src matlab :exports none
  figure;
  ax1 = subplot(2, 1, 1);
  hold on;
  plot(f, abs(tf_est), 'k-')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
  ylabel('Amplitude'); xlabel('Frequency [Hz]');
  hold off;

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(f, 180/pi*angle(tf_est), 'k-')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
  ylabel('Phase'); xlabel('Frequency [Hz]');
  hold off;

  linkaxes([ax1,ax2], 'x');
  xlim([10, 5e3]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/apa95ml_5kg_10V_tf.pdf', 'width', 'full', 'height', 'full');
#+end_src

#+name: fig:apa95ml_5kg_10V_tf
#+caption: Estimation of the transfer function from input voltage to displacement
#+RESULTS:
[[file:figs/apa95ml_5kg_10V_tf.png]]

** Comparison with the FEM model
#+begin_src matlab
  load('mat/fem_model_5kg.mat', 'Ghm');
#+end_src

#+begin_src matlab :exports none
  figure;
  ax1 = subplot(2, 1, 1);
  hold on;
  plot(f, abs(tf_est), 'DisplayName', 'Identification')
  plot(f, abs(squeeze(freqresp(Ghm, f, 'Hz'))), 'DisplayName', 'FEM')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
  ylabel('Amplitude'); xlabel('Frequency [Hz]');
  legend('location', 'northeast')
  hold off;

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(f, 180/pi*angle(tf_est))
  plot(f, 180/pi*angle(squeeze(freqresp(Ghm, f, 'Hz'))))
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
  ylabel('Phase'); xlabel('Frequency [Hz]');
  hold off;

  linkaxes([ax1,ax2], 'x');
  xlim([10, 5e3]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/apa95ml_5kg_comp_fem.pdf', 'width', 'full', 'height', 'full');
#+end_src

#+name: fig:apa95ml_5kg_comp_fem
#+caption: Comparison of the identified transfer function and the one estimated from the FE model
#+RESULTS:
[[file:figs/apa95ml_5kg_comp_fem.png]]

** Conclusion                                                        :ignore:
#+begin_important
The problem comes from the fact that the piezo is driven directly by the DAC that cannot deliver enought current.
In the next section, a current amplifier is used.
#+end_important

* Transfer Function Estimation using the PI Amplifier
** 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

** Load Data
#+begin_src matlab
  ht = load('./mat/huddle_test.mat', 't', 'u', 'y');
  load('./mat/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('mat/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]]

* Transfer function from force actuator to force sensor
** 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

** Load Data                                                         :ignore:
The data are loaded:
#+begin_src matlab
  a_ss = load('mat/apa95ml_5kg_1a_2s.mat', 't', 'u', 'y', 'v');
  aa_s = load('mat/apa95ml_5kg_2a_1s.mat', 't', 'u', 'y', 'v');
  load('mat/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]]

* IFF Tests
** 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

** First tests with few gains
#+begin_src matlab
  iff_g10 = load('./mat/apa95ml_iff_g10_res.mat', 'u', 't', 'y', 'v');
  iff_g100 = load('./mat/apa95ml_iff_g100_res.mat', 'u', 't', 'y', 'v');
  iff_of = load('./mat/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('./mat/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]]