#+TITLE: Amplifier Piezoelectric Actuator APA300ML - Test Bench
:DRAWER:
#+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>

#+BIND: org-latex-image-default-option "scale=1"
#+BIND: org-latex-image-default-width ""
#+BIND: org-latex-bib-compiler "biber"

#+LaTeX_CLASS: scrreprt
#+LaTeX_CLASS_OPTIONS: [a4paper, 10pt, DIV=12, parskip=full]
#+LaTeX_HEADER_EXTRA: \input{preamble.tex}
#+LATEX_HEADER_EXTRA: \addbibresource{ref.bib}

#+PROPERTY: header-args:matlab  :session *MATLAB*
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :results none
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :noweb yes
#+PROPERTY: header-args:matlab+ :tangle no
#+PROPERTY: header-args:matlab+ :mkdirp yes
#+PROPERTY: header-args:matlab+ :output-dir figs

#+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 file raw replace
#+PROPERTY: header-args:latex+ :buffer no
#+PROPERTY: header-args:latex+ :tangle no
#+PROPERTY: header-args:latex+ :eval no-export
#+PROPERTY: header-args:latex+ :exports results
#+PROPERTY: header-args:latex+ :mkdirp yes
#+PROPERTY: header-args:latex+ :output-dir figs
#+PROPERTY: header-args:latex+ :post pdf2svg(file=*this*, ext="png")
:END:

#+begin_export html
  <hr>
  <p>This report is also available as a <a href="./test-bench-apa300ml.pdf">pdf</a>.</p>
  <hr>
#+end_export

* Introduction                                                        :ignore:

The goal of this test bench is to extract all the important parameters of the Amplified Piezoelectric Actuator APA300ML.

This include:
- Stroke
- Stiffness
- Hysteresis
- Gain from the applied voltage $V_a$ to the generated Force $F_a$
- Gain from the sensor stack strain $\delta L$ to the generated voltage $V_s$
- Dynamical behavior

#+name: fig:apa300ML
#+caption: Picture of the APA300ML
#+attr_latex: :width 0.8\linewidth
[[file:figs/apa300ML.png]]

* Model of an Amplified Piezoelectric Actuator and Sensor

Consider a schematic of the Amplified Piezoelectric Actuator in Figure [[fig:apa_model_schematic]].

#+name: fig:apa_model_schematic
#+caption: Amplified Piezoelectric Actuator Schematic
[[file:figs/apa_model_schematic.png]]

A voltage $V_a$ applied to the actuator stacks will induce an actuator force $F_a$:
\begin{equation}
  F_a = g_a \cdot V_a
\end{equation}

A change of length $dl$ of the sensor stack will induce a voltage $V_s$:
\begin{equation}
  V_s = g_s \cdot dl
\end{equation}

We wish here to experimental measure $g_a$ and $g_s$.

The block-diagram model of the piezoelectric actuator is then as shown in Figure [[fig:apa-model-simscape-schematic]].

#+begin_src latex :file apa-model-simscape-schematic.pdf
  \begin{tikzpicture}
    \node[block={2.0cm}{2.0cm}, align=center] (model) at (0,0){Simscape\\Model};
    \node[block, left=1.0 of model] (ga){$g_a(s)$};
    \node[block, right=1.0 of model] (gs){$g_s(s)$};

    \draw[<-] (ga.west) -- node[midway, above]{$V_a$} node[midway, below]{$[V]$} ++(-1.0, 0);
    \draw[->] (ga.east) --node[midway, above]{$F_a$} node[midway, below]{$[N]$} (model.west);
    \draw[->] (model.east) --node[midway, above]{$dl$} node[midway, below]{$[m]$} (gs.west);
    \draw[->] (gs.east) -- node[midway, above]{$V_s$} node[midway, below]{$[V]$} ++(1.0, 0);
  \end{tikzpicture}
#+end_src

#+name: fig:apa-model-simscape-schematic
#+caption: Model of the APA with Simscape/Simulink
#+RESULTS:
[[file:figs/apa-model-simscape-schematic.png]]

* Geometrical Measurements
** Introduction                                                      :ignore:

The received APA are shown in Figure [[fig:received_apa]].

#+name: fig:received_apa
#+caption: Received APA
#+attr_latex: :width 0.6\linewidth
[[file:figs/IMG_20210224_143500.jpg]]

** 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/');
addpath('./matlab/');
#+end_src

#+begin_src matlab :eval no
addpath('./mat/');
#+end_src

** Measurement Setup

The flatness corresponding to the two interface planes are measured as shown in Figure [[fig:flatness_meas_setup]].

#+name: fig:flatness_meas_setup
#+caption: Measurement Setup
#+attr_latex: :width 0.6\linewidth
[[file:figs/IMG_20210224_143809.jpg]]

** Measurement Results

The height (Z) measurements at the 8 locations (4 points by plane) are defined below.
#+begin_src matlab
apa1  = 1e-6*[0, -0.5 , 3.5  , 3.5  , 42   , 45.5, 52.5 , 46];
apa2  = 1e-6*[0, -2.5 , -3   , 0    , -1.5 , 1   , -2   , -4];
apa3  = 1e-6*[0, -1.5 , 15   , 17.5 , 6.5  , 6.5 , 21   , 23];
apa4  = 1e-6*[0, 6.5  , 14.5 , 9    , 16   , 22  , 29.5 , 21];
apa5  = 1e-6*[0, -12.5, 16.5 , 28.5 , -43  , -52 , -22.5, -13.5];
apa6  = 1e-6*[0, -8   , -2   , 5    , -57.5, -62 , -55.5, -52.5];
apa7  = 1e-6*[0, 19.5 , -8   , -29.5, 75   , 97.5, 70   , 48];
apa7b = 1e-6*[0, 9    , -18.5, -30  , 31   , 46.5, 16.5 , 7.5];
apa = {apa1, apa2, apa3, apa4, apa5, apa6, apa7b};
#+end_src

The X/Y Positions of the 8 measurement points are defined below.
#+begin_src matlab
W = 20e-3;   % Width [m]
L = 61e-3;   % Length [m]
d = 1e-3;    % Distance from border [m]
l = 15.5e-3; % [m]

pos = [[-L/2 + d; W/2 - d], [-L/2 + l - d; W/2 - d], [-L/2 + l - d; -W/2 + d], [-L/2 + d; -W/2 + d], [L/2 - l + d; W/2 - d], [L/2 - d; W/2 - d], [L/2 - d; -W/2 + d], [L/2 - l + d; -W/2 + d]];
#+end_src

Finally, the flatness is estimated by fitting a plane through the 8 points using the =fminsearch= command.
#+begin_src matlab
apa_d = zeros(1, 7);
for i = 1:7
    fun = @(x)max(abs(([pos; apa{i}]-[0;0;x(1)])'*([x(2:3);1]/norm([x(2:3);1]))));
    x0 = [0;0;0];
    [x, min_d] = fminsearch(fun,x0);
    apa_d(i) = min_d;
end
#+end_src

The obtained flatness are shown in Table [[tab:flatness_meas]].

#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
data2orgtable(1e6*apa_d', {'APA 1', 'APA 2', 'APA 3', 'APA 4', 'APA 5', 'APA 6', 'APA 7'}, {'*Flatness* $[\mu m]$'}, ' %.1f ');
#+end_src

#+name: tab:flatness_meas
#+caption: Estimated flatness
#+attr_latex: :environment tabularx :width 0.25\linewidth :align lc
#+attr_latex: :center t :booktabs t :float t
#+RESULTS:
|       | *Flatness* $[\mu m]$ |
|-------+----------------------|
| APA 1 |                  8.9 |
| APA 2 |                  3.1 |
| APA 3 |                  9.1 |
| APA 4 |                  3.0 |
| APA 5 |                  1.9 |
| APA 6 |                  7.1 |
| APA 7 |                 18.7 |

* Electrical Measurements

#+begin_note
The capacitance of the stacks is measure with the [[https://www.gwinstek.com/en-global/products/detail/LCR-800][LCR-800 Meter]] ([[file:doc/DS_LCR-800_Series_V2_E.pdf][doc]])
#+end_note

#+name: fig:LCR_meter
#+caption: LCR Meter used for the measurements
#+attr_latex: :width 0.9\linewidth
[[file:figs/IMG_20210312_120337.jpg]]

The excitation frequency is set to be 1kHz.

#+name: tab:apa300ml_capacitance
#+caption: Capacitance measured with the LCR meter. The excitation signal is a sinus at 1kHz
#+attr_latex: :environment tabularx :width 0.5\linewidth :align lcc
#+attr_latex: :center t :booktabs t :float t
|       | *Sensor Stack* | *Actuator Stacks* |
|-------+----------------+-------------------|
| APA 1 |           5.10 |             10.03 |
| APA 2 |           4.99 |              9.85 |
| APA 3 |           1.72 |              5.18 |
| APA 4 |           4.94 |              9.82 |
| APA 5 |           4.90 |              9.66 |
| APA 6 |           4.99 |              9.91 |
| APA 7 |           4.85 |              9.85 |

#+begin_warning
There is clearly a problem with APA300ML number 3
#+end_warning

* Stroke measurement
** Introduction                                                      :ignore:

We here wish to estimate the stroke of the APA.

To do so, one side of the APA is fixed, and a displacement probe is located on the other side as shown in Figure [[fig:stroke_test_bench]].

Then, a voltage is applied on either one or two stacks using a DAC and a voltage amplifier.

#+begin_note
Here are the documentation of the equipment used for this test bench:
- *Voltage Amplifier*: [[file:doc/PD200-V7-R1.pdf][PD200]] with a gain of 20
- *16bits DAC*: [[file:doc/IO131-OEM-Datasheet.pdf][IO313 Speedgoat card]]
- *Displacement Probe*: [[file:doc/Millimar--3723046--BA--C1208-C1216-C1240--FR--2016-11-08.pdf][Millimar C1216 electronics]] and [[file:doc/tmp3m0cvmue_7888038c-cdc8-48d8-a837-35de02760685.pdf][Millimar 1318 probe]]
#+end_note

#+name: fig:stroke_test_bench
#+caption: Bench to measured the APA stroke
#+attr_latex: :width 0.9\linewidth
[[file:figs/CE0EF55E-07B7-461B-8CDB-98590F68D15B.jpeg]]

** 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/');
addpath('./matlab/');
#+end_src

#+begin_src matlab :eval no
addpath('./mat/');
#+end_src

** Voltage applied on one stack

Let's first look at the relation between the voltage applied to *one* stack to the displacement of the APA as measured by the displacement probe.

#+begin_src matlab :exports none
apa300ml_1s = {};
for i = 1:7
    apa300ml_1s(i) = {load(['mat/stroke_apa_1stacks_' num2str(i) '.mat'], 't', 'V', 'd')};
end
#+end_src

#+begin_src matlab :exports none
for i = 1:7
    t = apa300ml_1s{i}.t;
    apa300ml_1s{i}.d = apa300ml_1s{i}.d - mean(apa300ml_1s{i}.d(t > 1.9 & t < 2.0));
    apa300ml_1s{i}.d = apa300ml_1s{i}.d(t > 2.0 & t < 10.0);
    apa300ml_1s{i}.V = apa300ml_1s{i}.V(t > 2.0 & t < 10.0);
    apa300ml_1s{i}.t = apa300ml_1s{i}.t(t > 2.0 & t < 10.0);
end
#+end_src

The applied voltage is shown in Figure [[fig:apa_stroke_voltage_time]].

#+begin_src matlab :exports none
figure;
plot(apa300ml_1s{1}.t, 20*apa300ml_1s{1}.V)
xlabel('Time [s]'); ylabel('Voltage [V]');
ylim([-20,160]); yticks([-20     0    20    40    60    80   100   120   140   160]);
#+end_src

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

#+name: fig:apa_stroke_voltage_time
#+caption: Applied voltage as a function of time
#+RESULTS:
[[file:figs/apa_stroke_voltage_time.png]]

The obtained displacement is shown in Figure [[fig:apa_stroke_time_1s]].
The displacement is set to zero at initial time when the voltage applied is -20V.

#+begin_src matlab :exports none
figure;
hold on;
for i = 1:7
    plot(apa300ml_1s{i}.t, 1e6*apa300ml_1s{i}.d, 'DisplayName', sprintf('APA %i', i))
end
hold off;
xlabel('Time [s]'); ylabel('Displacement [$\mu m$]')
legend('location', 'southeast', 'FontSize', 8)
#+end_src

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

#+name: fig:apa_stroke_time_1s
#+caption: Displacement as a function of time for all the APA300ML
#+RESULTS:
[[file:figs/apa_stroke_time_1s.png]]

Finally, the displacement is shown as a function of the applied voltage in Figure [[fig:apa_d_vs_V_1s]].
We can clearly see that there is a problem with the APA 3.
Also, there is a large hysteresis.

#+begin_src matlab :exports none
figure;
hold on;
for i = 1:7
    plot(20*apa300ml_1s{i}.V, 1e6*apa300ml_1s{i}.d, 'DisplayName', sprintf('APA %i', i))
end
hold off;
xlabel('Voltage [V]'); ylabel('Displacement [$\mu m$]')
legend('location', 'southwest', 'FontSize', 8)
xlim([-20, 160]); ylim([-140, 0]);
#+end_src

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

#+name: fig:apa_d_vs_V_1s
#+caption: Displacement as a function of the applied voltage
#+RESULTS:
[[file:figs/apa_d_vs_V_1s.png]]

#+begin_important
We can clearly see from Figure [[fig:apa_d_vs_V_1s]] that there is a problem with the APA number 3.
#+end_important

** Voltage applied on two stacks

Now look at the relation between the voltage applied to the *two* other stacks to the displacement of the APA as measured by the displacement probe.

#+begin_src matlab :exports none
apa300ml_2s = {};
for i = 1:7
    apa300ml_2s(i) = {load(['mat/stroke_apa_2stacks_' num2str(i) '.mat'], 't', 'V', 'd')};
end
#+end_src

#+begin_src matlab :exports none
for i = 1:7
    t = apa300ml_2s{i}.t;
    apa300ml_2s{i}.d = apa300ml_2s{i}.d - mean(apa300ml_2s{i}.d(t > 1.9 & t < 2.0));
    apa300ml_2s{i}.d = apa300ml_2s{i}.d(t > 2.0 & t < 10.0);
    apa300ml_2s{i}.V = apa300ml_2s{i}.V(t > 2.0 & t < 10.0);
    apa300ml_2s{i}.t = apa300ml_2s{i}.t(t > 2.0 & t < 10.0);
end
#+end_src

The obtained displacement is shown in Figure [[fig:apa_stroke_time_2s]].
The displacement is set to zero at initial time when the voltage applied is -20V.

#+begin_src matlab :exports none
figure;
hold on;
for i = 1:7
    plot(apa300ml_2s{i}.t, 1e6*apa300ml_2s{i}.d, 'DisplayName', sprintf('APA %i', i))
end
hold off;
xlabel('Time [s]'); ylabel('Displacement [$\mu m$]')
legend('location', 'southeast', 'FontSize', 8)
ylim([-250, 0]);
#+end_src

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

#+name: fig:apa_stroke_time_2s
#+caption: Displacement as a function of time for all the APA300ML
#+RESULTS:
[[file:figs/apa_stroke_time_2s.png]]

Finally, the displacement is shown as a function of the applied voltage in Figure [[fig:apa_d_vs_V_2s]].
We can clearly see that there is a problem with the APA 3.
Also, there is a large hysteresis.

#+begin_src matlab :exports none
figure;
hold on;
for i = 1:7
    plot(20*apa300ml_2s{i}.V, 1e6*apa300ml_2s{i}.d, 'DisplayName', sprintf('APA %i', i))
end
hold off;
xlabel('Voltage [V]'); ylabel('Displacement [$\mu m$]')
legend('location', 'southwest', 'FontSize', 8)
xlim([-20, 160]); ylim([-250, 0]);
#+end_src

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

#+name: fig:apa_d_vs_V_2s
#+caption: Displacement as a function of the applied voltage
#+RESULTS:
[[file:figs/apa_d_vs_V_2s.png]]

** Voltage applied on all three stacks

Finally, we can combine the two measurements to estimate the relation between the displacement and the voltage applied to the *three* stacks (Figure [[fig:apa_d_vs_V_3s]]).

#+begin_src matlab :exports none
apa300ml_3s = {};
for i = 1:7
    apa300ml_3s(i) = apa300ml_1s(i);
    apa300ml_3s{i}.d = apa300ml_1s{i}.d + apa300ml_2s{i}.d;
end
#+end_src

#+begin_src matlab :exports none
figure;
hold on;
for i = 1:7
    plot(20*apa300ml_3s{i}.V, 1e6*apa300ml_3s{i}.d, 'DisplayName', sprintf('APA %i', i))
end
hold off;
xlabel('Voltage [V]'); ylabel('Displacement [$\mu m$]')
legend('location', 'southwest', 'FontSize', 8)
xlim([-20, 160]); ylim([-400, 0]);
#+end_src

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

#+name: fig:apa_d_vs_V_3s
#+caption: Displacement as a function of the applied voltage
#+RESULTS:
[[file:figs/apa_d_vs_V_3s.png]]

The obtained maximum stroke for all the APA are summarized in Table [[tab:apa_measured_stroke]].

#+begin_src matlab :exports none
apa300ml_stroke = zeros(1, 7);
for i = 1:7
    apa300ml_stroke(i) = max(apa300ml_3s{i}.d) - min(apa300ml_3s{i}.d);
end
#+end_src

#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
data2orgtable(1e6*apa300ml_stroke', {'APA 1', 'APA 2', 'APA 3', 'APA 4', 'APA 5', 'APA 6', 'APA 7'}, {'*Stroke* $[\mu m]$'}, ' %.1f ');
#+end_src

#+name: tab:apa_measured_stroke
#+caption: Measured maximum stroke
#+attr_latex: :environment tabularx :width 0.25\linewidth :align lc
#+attr_latex: :center t :booktabs t :float t
#+RESULTS:
|       | *Stroke* $[\mu m]$ |
|-------+--------------------|
| APA 1 |              373.2 |
| APA 2 |              365.5 |
| APA 3 |              181.7 |
| APA 4 |              359.7 |
| APA 5 |              361.5 |
| APA 6 |              363.9 |
| APA 7 |              358.4 |

* TODO Stiffness measurement
** 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/');
addpath('./matlab/');
#+end_src

#+begin_src matlab :eval no
addpath('./mat/');
#+end_src

** APA test
#+begin_src matlab
load('meas_stiff_apa_1_x.mat', 't', 'F', 'd');
#+end_src

#+begin_src matlab
figure;
plot(t, F)
#+end_src

#+begin_src matlab
%% Automatic Zero of the force
F = F - mean(F(t > 0.1 & t < 0.3));

%% Start measurement at t = 0.2 s
d = d(t > 0.2);
F = F(t > 0.2);
t = t(t > 0.2); t = t - t(1);
#+end_src

#+begin_src matlab
i_l_start = find(F > 0.3, 1, 'first');
[~, i_l_stop] = max(F);
#+end_src

#+begin_src matlab
F_l = F(i_l_start:i_l_stop);
d_l = d(i_l_start:i_l_stop);
#+end_src

#+begin_src matlab
fit_l = polyfit(F_l, d_l, 1);

% %% Reset displacement based on fit
% d = d - fit_l(2);
% fit_s(2) = fit_s(2) - fit_l(2);
% fit_l(2) = 0;

% %% Estimated Stroke
% F_max = fit_s(2)/(fit_l(1) - fit_s(1));
% d_max = fit_l(1)*F_max;
#+end_src

#+begin_src matlab
h^2/fit_l(1)
#+end_src

#+begin_src matlab
figure;
hold on;
plot(F,d,'k')
plot(F_l, d_l)
plot(F_l, F_l*fit_l(1) + fit_l(2), '--')
#+end_src
* Test-Bench Description

#+begin_note
Here are the documentation of the equipment used for this test bench:
- Voltage Amplifier: [[file:doc/PD200-V7-R1.pdf][PD200]]
- Amplified Piezoelectric Actuator: [[file:doc/APA300ML.pdf][APA300ML]]
- DAC/ADC: Speedgoat [[file:doc/IO131-OEM-Datasheet.pdf][IO313]]
- Encoder: [[file:doc/L-9517-9678-05-A_Data_sheet_VIONiC_series_en.pdf][Renishaw Vionic]] and used [[file:doc/L-9517-9862-01-C_Data_sheet_RKLC_EN.pdf][Ruler]]
- Interferometer: [[https://www.attocube.com/en/products/laser-displacement-sensor/displacement-measuring-interferometer][Attocube IDS3010]]
#+end_note

#+name: fig:test_bench_apa_alone
#+caption: Schematic of the Test Bench
[[file:figs/test_bench_apa_alone.png]]

* Measurement Procedure
<<sec:measurement_procedure>>
** Introduction                                                      :ignore:

** Stroke Measurement

Using the PD200 amplifier, output a voltage:
\[ V_a = 65 + 85 \sin(2\pi \cdot t) \]
To have a quasi-static excitation between -20 and 150V.

As the gain of the PD200 amplifier is 20, the DAC output voltage should be:
\[ V_{dac}(t) = 3.25 + 4.25\sin(2\pi \cdot t) \]

Verify that the voltage offset of the PD200 is zero!

Measure the output vertical displacement $d$ using the interferometer.

Then, plot $d$ as a function of $V_a$, and perform a linear regression.
Conclude on the obtained stroke.

** Stiffness Measurement

Add some (known) weight $\delta m g$ on the suspended mass and measure the deflection $\delta d$.
This can be tested when the piezoelectric stacks are open-circuit.

As the stiffness will be around $k \approx 10^6 N/m$, an added mass of $m \approx 100g$ will induce a static deflection of $\approx 1\mu m$ which should be large enough for a precise measurement using the interferometer.

Then the obtained stiffness is:
\begin{equation}
  k = \frac{\delta m g}{\delta d}
\end{equation}

** Hysteresis measurement

Supply a quasi static sinusoidal excitation $V_a$ at different voltages.

The offset should be 65V, and the sin amplitude can range from 1V up to 85V.

For each excitation amplitude, the vertical displacement $d$ of the mass is measured.

Then, $d$ is plotted as a function of $V_a$ for all the amplitudes.

#+name: fig:expected_hysteresis
#+caption: Expected Hysteresis cite:poel10_explor_activ_hard_mount_vibrat
#+attr_latex: :width 0.8\linewidth
[[file:figs/expected_hysteresis.png]]

** Piezoelectric Actuator Constant

Using the measurement test-bench, it is rather easy the determine the static gain between the applied voltage $V_a$ to the induced displacement $d$.
Use a quasi static (1Hz) excitation signal $V_a$ on the piezoelectric stack and measure the vertical displacement $d$.
Perform a linear regression to obtain:
\begin{equation}
  d = g_{d/V_a} \cdot V_a
\end{equation}

Using the Simscape model of the APA, it is possible to determine the static gain between the actuator force $F_a$ to the induced displacement $d$:
\begin{equation}
  d = g_{d/F_a} \cdot F_a
\end{equation}

From the two gains, it is then easy to determine $g_a$:
\begin{equation}
  g_a = \frac{F_a}{V_a} = \frac{F_a}{d} \cdot \frac{d}{V_a} = \frac{g_{d/V_a}}{g_{d/F_a}}
\end{equation}

** Piezoelectric Sensor Constant

From a quasi static excitation of the piezoelectric stack, measure the gain from $V_a$ to $V_s$:
\begin{equation}
  V_s = g_{V_s/V_a} V_a
\end{equation}

Note here that there is an high pass filter formed by the piezo capacitor and parallel resistor.
The excitation frequency should then be in between the cut-off frequency of this high pass filter and the first resonance.

Alternatively, the gain can be computed from the dynamical identification and taking the gain at the wanted frequency.

Using the simscape model, compute the static gain from the actuator force $F_a$ to the strain of the sensor stack $dl$:
\begin{equation}
  dl = g_{dl/F_a} F_a
\end{equation}

Then, the static gain from the sensor stack strain $dl$ to the general voltage $V_s$ is:
\begin{equation}
  g_s = \frac{V_s}{dl} = \frac{V_s}{V_a} \cdot \frac{V_a}{F_a} \cdot \frac{F_a}{dl} = \frac{g_{V_s/V_a}}{g_a \cdot g_{dl/F_a}}
\end{equation}

Alternatively, we could impose an external force to add strain in the APA that should be equally present in all the 3 stacks and equal to 1/5 of the vertical strain.
This external force can be some weight added, or a piezo in parallel.

** Capacitance Measurement

Measure the capacitance of the 3 stacks individually using a precise multi-meter.

** Dynamical Behavior

Perform a system identification from $V_a$ to the measured displacement $d$ by the interferometer and by the encoder, and to the generated voltage $V_s$.

This can be performed using different excitation signals.

This can also be performed with and without the encoder fixed to the APA.

** Compare the results obtained for all 7 APA300ML

Compare all the obtained parameters for all the test APA.

* FRF measurement
** Introduction                                                      :ignore:

- [ ] Schematic of the measurement

| Variable |                              | Unit | Hardware                      |
|----------+------------------------------+------+-------------------------------|
| =Va=     | Output DAC voltage           | [V]  | DAC - Ch. 1 => PD200 => APA   |
| =Vs=     | Measured stack voltage (ADC) | [V]  | APA => ADC - Ch. 1            |
| =de=     | Encoder Measurement          | [m]  | PEPU Ch. 1 - IO318(1) - Ch. 1 |
| =da=     | Attocube Measurement         | [m]  | PEPU Ch. 2 - IO318(1) - Ch. 2 |
| =t=      | Time                         | [s]  |                               |

** 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/src/');
addpath('./matlab/');
#+end_src

#+begin_src matlab :eval no
addpath('./src/');
#+end_src

** =frf_setup.m= - Measurement Setup
:PROPERTIES:
:header-args:matlab: :tangle matlab/frf_setup.m
:header-args:matlab+: :comments no
:END:
#+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 :eval no :exports none
addpath('./src/');
#+end_src

First is defined the sampling frequency:
#+begin_src matlab
%% Simulation configuration
Fs   = 10e3; % Sampling Frequency [Hz]
Ts   = 1/Fs; % Sampling Time [s]
#+end_src

#+begin_src matlab
%% Data record configuration
Trec_start = 5;  % Start time for Recording [s]
Trec_dur   = 100; % Recording Duration [s]
#+end_src

#+begin_src matlab
Tsim = 2*Trec_start + Trec_dur;  % Simulation Time [s]
#+end_src

The maximum excitation voltage at resonance is 9Vrms, therefore corresponding to 0.6V of output DAC voltage.
#+begin_src matlab
%% Sweep Sine
gc = 0.1;
xi = 0.5;
wn = 2*pi*94.3;

% Notch filter at the resonance of the APA
G_sweep = 0.2*(s^2 + 2*gc*xi*wn*s + wn^2)/(s^2 + 2*xi*wn*s + wn^2);

V_sweep = generateSweepExc('Ts',      Ts, ...
                           'f_start', 10, ...
                           'f_end',   1e3, ...
                           'V_mean',  3.25, ...
                           't_start', Trec_start, ...
                           'exc_duration', Trec_dur, ...
                           'sweep_type',   'log', ...
                           'V_exc', G_sweep*1/(1 + s/2/pi/500));
#+end_src

#+begin_src matlab :exports none :tangle no
figure;
tiledlayout(1, 2, 'TileSpacing', 'Normal', 'Padding', 'None');

ax1 = nexttile;
plot(V_sweep(1,:), V_sweep(2,:));
xlabel('Time [s]'); ylabel('Amplitude [V]');

ax2 = nexttile;
win = hanning(floor(length(V_sweep(2,:))/80));
[pxx, f] = pwelch(V_sweep(2,:), win, 0, [], Fs);
plot(f, pxx)
xlabel('Frequency [Hz]'); ylabel('Power Spectral Density [$V^2/Hz$]');
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlim([1, Fs/2]); ylim([1e-10, 1e0]);
#+end_src

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

#+name: fig:frf_meas_sweep_excitation
#+caption: Example of Sweep Sin excitation signal
#+RESULTS:
[[file:figs/frf_meas_sweep_excitation.png]]


A white noise excitation signal can be very useful in order to obtain a first idea of the plant FRF.
The gain can be gradually increased until satisfactory output is obtained.
#+begin_src matlab
%% Shaped Noise
V_noise = generateShapedNoise('Ts', 1/Fs, ...
                              'V_mean', 3.25, ...
                              't_start', Trec_start, ...
                              'exc_duration', Trec_dur, ...
                              'smooth_ends', true, ...
                              'V_exc', 0.05/(1 + s/2/pi/10));
#+end_src

#+begin_src matlab :exports none :tangle no
figure;
tiledlayout(1, 2, 'TileSpacing', 'Normal', 'Padding', 'None');

ax1 = nexttile;
plot(V_noise(1,:), V_noise(2,:));
xlabel('Time [s]'); ylabel('Amplitude [V]');

ax2 = nexttile;
win = hanning(floor(length(V_noise)/8));
[pxx, f] = pwelch(V_noise(2,:), win, 0, [], Fs);
plot(f, pxx)
xlabel('Frequency [Hz]'); ylabel('Power Spectral Density [$V^2/Hz$]');
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlim([1, Fs/2]); ylim([1e-10, 1e0]);
#+end_src

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

#+name: fig:frf_meas_noise_excitation
#+caption: Example of Shaped noise excitation signal
#+RESULTS:
[[file:figs/frf_meas_noise_excitation.png]]

#+begin_src matlab
%% Sinus excitation with increasing amplitude
V_sin = generateSinIncreasingAmpl('Ts', 1/Fs, ...
                                    'V_mean', 3.25, ...
                                    'sin_ampls', [0.1, 0.2, 0.4, 1, 2, 4], ...
                                    'sin_period', 1, ...
                                    'sin_num', 5, ...
                                    't_start', 10, ...
                                    'smooth_ends', true);
#+end_src

#+begin_src matlab :exports none :tangle no
figure;
plot(V_sin(1,:), V_sin(2,:));
xlabel('Time [s]'); ylabel('Amplitude [V]');
#+end_src

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

#+name: fig:frf_meas_sin_excitation
#+caption: Example of Shaped noise excitation signal
#+RESULTS:
[[file:figs/frf_meas_sin_excitation.png]]

Then, we select the wanted excitation signal.
#+begin_src matlab
%% Select the excitation signal
V_exc = timeseries(V_noise(2,:), V_noise(1,:));
#+end_src

#+begin_src matlab :exports none :eval no
figure;
tiledlayout(1, 2, 'TileSpacing', 'Normal', 'Padding', 'None');

ax1 = nexttile;
plot(V_exc(1,:), V_exc(2,:));
xlabel('Time [s]'); ylabel('Amplitude [V]');

ax2 = nexttile;
win = hanning(floor(length(V_exc)/8));
[pxx, f] = pwelch(V_exc(2,:), win, 0, [], Fs);
plot(f, pxx)
xlabel('Frequency [Hz]'); ylabel('Power Spectral Density [$V^2/Hz$]');
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlim([1, Fs/2]); ylim([1e-10, 1e0]);
#+end_src

#+begin_src matlab
%% Save data that will be loaded in the Simulink file
save('./frf_data.mat', 'Fs', 'Ts', 'Tsim', 'Trec_start', 'Trec_dur', 'V_exc');
#+end_src

** =frf_save.m= - Save Data
:PROPERTIES:
:header-args: :tangle matlab/frf_save.m
:header-args:matlab+: :comments no
:END:

First, we get data from the Speedgoat:
#+begin_src matlab
tg = slrt;

f = SimulinkRealTime.openFTP(tg);
mget(f, 'data/data.dat');
close(f);
#+end_src

And we load the data on the Workspace:
#+begin_src matlab
data = SimulinkRealTime.utils.getFileScopeData('data/data.dat').data;

da = data(:, 1); % Excitation Voltage (input of PD200) [V]
de = data(:, 2); % Measured voltage (force sensor) [V]
Vs = data(:, 3); % Measurment displacement (encoder) [m]
Va = data(:, 4); % Measurement displacement (attocube) [m]
t  = data(:, end); % Time [s]
#+end_src

And we save this to a =mat= file:
#+begin_src matlab
apa_number = 1;

save(sprintf('mat/frf_data_%i_huddle.mat', apa_number), 't', 'Va', 'Vs', 'de', 'da');
#+end_src

** Measurements on APA 1
*** Introduction                                                    :ignore:
Measurements are first performed on the APA number 1.

*** 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/');
addpath('./matlab/src/');
addpath('./matlab/');
#+end_src

#+begin_src matlab :eval no
addpath('./mat/');
addpath('./src/');
#+end_src

*** Huddle Test
#+begin_src matlab
load(sprintf('frf_data_%i_huddle.mat', 1), 't', 'Va', 'Vs', 'de', 'da');
#+end_src

#+begin_src matlab :exports none
figure;
tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(t, da, 'DisplayName', '$d_a$ - Attocube')
plot(t, de, 'DisplayName', '$d_e$ - Encoder')
hold off;
xlabel('Time [s]');
ylabel('Displacement [m]');
legend('location', 'northeast')

ax2 = nexttile;
hold on;
plot(t, Vs, 'DisplayName', '$V_s$ - Force Sensor')
hold off;
xlabel('Time [s]');
ylabel('Voltage [V]');
legend('location', 'northeast')
#+end_src

#+begin_src matlab
%% Parameter for Spectral analysis
Ts = (t(end) - t(1))/(length(t)-1);
Fs = 1/Ts;

win = hanning(ceil(10*Fs)); % Hannning Windows
#+end_src

#+begin_src matlab
[phh_da, f] = pwelch(da - mean(da), win, [], [], 1/Ts);
[phh_de, ~] = pwelch(de - mean(de), win, [], [], 1/Ts);
[phh_Vs, ~] = pwelch(Vs - mean(Vs), win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
figure;
tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f, phh_da, 'DisplayName', '$d_a$ - Attocube')
plot(f, phh_de, 'DisplayName', '$d_e$ - Encoder')
hold off;
xlabel('Time [s]'); ylabel('ASD [$m/\sqrt{Hz}$]');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
legend('location', 'northeast')
xlim([5e-1, 5e3]);

ax2 = nexttile;
hold on;
plot(f, phh_Vs, 'DisplayName', '$V_s$ - Force Sensor')
hold off;
xlabel('Time [s]'); ylabel('ASD [$V/\sqrt{Hz}$]');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
legend('location', 'northeast')
xlim([5e-1, 5e3]);
#+end_src

*** First identification with Noise
#+begin_src matlab
load(sprintf('mat/frf_data_%i_noise.mat', 1), 't', 'Va', 'Vs', 'da', 'de')
#+end_src

#+begin_src matlab
[pxx_da, f] = pwelch(da, win, [], [], 1/Ts);
[pxx_de, ~] = pwelch(de, win, [], [], 1/Ts);
[pxx_Vs, ~] = pwelch(Vs, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
figure;
tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f, sqrt(phh_da), 'DisplayName', 'Huddle')
plot(f, sqrt(pxx_da), 'DisplayName', 'Noise')
hold off;
xlabel('Time [s]'); ylabel('ASD Attocube [$m/\sqrt{Hz}$]');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
legend('location', 'northeast')

ax2 = nexttile;
hold on;
plot(f, sqrt(phh_Vs), 'DisplayName', 'Huddle')
plot(f, sqrt(pxx_Vs), 'DisplayName', 'Noise')
hold off;
xlabel('Time [s]'); ylabel('ASD [$V/\sqrt{Hz}$]');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
legend('location', 'northeast')

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

#+begin_src matlab
[G_dvf, f] = tfestimate(Va, de, win, [], [], 1/Ts);
[G_d,   ~] = tfestimate(Va, da, win, [], [], 1/Ts);
[G_iff, ~] = tfestimate(Va, Vs, win, [], [], 1/Ts);

[coh_dvf, ~] = mscohere(Va, de, win, [], [], 1/Ts);
[coh_d,   ~] = mscohere(Va, da, win, [], [], 1/Ts);
[coh_iff, ~] = mscohere(Va, Vs, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
%%
figure;
hold on;
plot(f, coh_dvf);
plot(f, coh_d);
plot(f, coh_iff);
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlim([1, 5e3]); ylim([0, 1]);
#+end_src

#+begin_src matlab :exports none
%%
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f, abs(G_d),   'DisplayName', 'Attocube');
plot(f, abs(G_dvf), 'DisplayName', 'Encoder');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d/V_a$ [m/V]'); set(gca, 'XTickLabel',[]);
hold off;
ylim([1e-8, 1e-3]);
legend('location', 'northeast');

ax2 = nexttile;
hold on;
plot(f, 180/pi*angle(G_d));
plot(f, 180/pi*angle(G_dvf));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

#+begin_src matlab :exports none
%%
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
plot(f, abs(G_iff));
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]);
hold off;

ax2 = nexttile;
plot(f, 180/pi*angle(G_iff));
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

*** Second identification with Sweep and high frequency noise
#+begin_src matlab
load(sprintf('mat/frf_data_%i_sweep.mat', 1), 't', 'Va', 'Vs', 'da', 'de')
#+end_src

#+begin_src matlab
[pxx_da, f] = pwelch(da, win, [], [], 1/Ts);
[pxx_de, ~] = pwelch(de, win, [], [], 1/Ts);
[pxx_Vs, ~] = pwelch(Vs, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
figure;
tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f, sqrt(phh_da), 'DisplayName', 'Huddle')
plot(f, sqrt(pxx_da), 'DisplayName', 'Noise')
hold off;
xlabel('Time [s]'); ylabel('ASD Attocube [$m/\sqrt{Hz}$]');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
legend('location', 'northeast')

ax2 = nexttile;
hold on;
plot(f, sqrt(phh_Vs), 'DisplayName', 'Huddle')
plot(f, sqrt(pxx_Vs), 'DisplayName', 'Noise')
hold off;
xlabel('Time [s]'); ylabel('ASD [$V/\sqrt{Hz}$]');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
legend('location', 'northeast')

linkaxes([ax1,ax2],'x');
xlim([1e1, 2e3]);
#+end_src

#+begin_src matlab
[G_dvf, f] = tfestimate(Va, de, win, [], [], 1/Ts);
[G_d,   ~] = tfestimate(Va, da, win, [], [], 1/Ts);
[G_iff, ~] = tfestimate(Va, Vs, win, [], [], 1/Ts);

[coh_dvf, ~] = mscohere(Va, de, win, [], [], 1/Ts);
[coh_d,   ~] = mscohere(Va, da, win, [], [], 1/Ts);
[coh_iff, ~] = mscohere(Va, Vs, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab
load(sprintf('mat/frf_data_%i_noise_high_freq.mat', 1), 't', 'Va', 'Vs', 'da', 'de')
#+end_src

#+begin_src matlab
[phf_da, f] = pwelch(da, win, [], [], 1/Ts);
[phf_de, ~] = pwelch(de, win, [], [], 1/Ts);
[phf_Vs, ~] = pwelch(Vs, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab
[cohhf_dvf, ~] = mscohere(Va, de, win, [], [], 1/Ts);
[cohhf_d,   ~] = mscohere(Va, da, win, [], [], 1/Ts);
[cohhf_iff, ~] = mscohere(Va, Vs, win, [], [], 1/Ts);

[Ghf_dvf, f] = tfestimate(Va, de, win, [], [], 1/Ts);
[Ghf_d,   ~] = tfestimate(Va, da, win, [], [], 1/Ts);
[Ghf_iff, ~] = tfestimate(Va, Vs, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
%%
figure;
hold on;
plot(f, coh_dvf);
plot(f, cohhf_dvf);
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlim([10, 2e3]); ylim([0, 1]);
#+end_src


#+begin_src matlab :exports none
%%
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f(f<350), abs(G_dvf(  f<350)), 'DisplayName', 'Attocube');
plot(f(f>350), abs(Ghf_dvf(f>350)), 'DisplayName', 'Encoder');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d/V_a$ [m/V]'); set(gca, 'XTickLabel',[]);
hold off;
ylim([1e-8, 1e-3]);
legend('location', 'northeast');

ax2 = nexttile;
hold on;
plot(f(f<350), 180/pi*angle(G_dvf(  f<350)));
plot(f(f>350), 180/pi*angle(Ghf_dvf(f>350)));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

#+begin_src matlab :exports none
%%
figure;
hold on;
plot(f, coh_dvf);
plot(f, coh_d);
plot(f, coh_iff);
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlim([10, 2e3]); ylim([0, 1]);
#+end_src

#+begin_src matlab :exports none
%%
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f, abs(G_d),   'DisplayName', 'Attocube');
plot(f, abs(G_dvf), 'DisplayName', 'Encoder');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d/V_a$ [m/V]'); set(gca, 'XTickLabel',[]);
hold off;
ylim([1e-8, 1e-3]);
legend('location', 'northeast');

ax2 = nexttile;
hold on;
plot(f, 180/pi*angle(G_d));
plot(f, 180/pi*angle(G_dvf));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

#+begin_src matlab :exports none
%%
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
plot(f, abs(G_iff));
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $V_s/V_a$ [V/V]'); set(gca, 'XTickLabel',[]);
hold off;

ax2 = nexttile;
plot(f, 180/pi*angle(G_iff));
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

*** Extract Parameters (Actuator/Sensor constants)

Quasi static gain between $d$ and $V_a$:
#+begin_src matlab
g_d_Va = mean(abs(G_dvf(f > 10 & f < 15)));
#+end_src

#+begin_src matlab :results value replace :exports results
sprintf('g_d_Va = %.1e [m/V]', g_d_Va)
#+end_src

#+RESULTS:
: g_d_Va = 1.7e-05 [m/V]

Quasi static gain between $V_s$ and $V_a$:
#+begin_src matlab
g_Vs_Va = mean(abs(G_iff(f > 10 & f < 15)));
#+end_src

#+begin_src matlab :results value replace :exports results
sprintf('g_Vs_Va = %.1e [V/V]', g_Vs_Va)
#+end_src

#+RESULTS:
: g_Vs_Va = 5.7e-01 [V/V]

** Comparison of all the APA
*** 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/');
addpath('./matlab/src/');
addpath('./matlab/');
#+end_src

#+begin_src matlab :eval no
addpath('./mat/');
addpath('./src/');
#+end_src

*** Stiffness - Comparison of the APA
In order to estimate the stiffness of the APA, a weight with known mass $m_a$ is added on top of the suspended granite and the deflection $d_e$ is measured using an encoder.
The APA stiffness is then:
\begin{equation}
  k_{\text{apa}} = \frac{m_a g}{d}
\end{equation}

Here are the number of the APA that have been measured:
#+begin_src matlab
apa_nums = [1 2 4 5 6 7 8];
#+end_src

The data are loaded.
#+begin_src matlab
apa_mass = {};
for i = 1:length(apa_nums)
    apa_mass(i) = {load(sprintf('frf_data_%i_add_mass_closed_circuit.mat', apa_nums(i)), 't', 'de')};
    % The initial displacement is set to zero
    apa_mass{i}.de = apa_mass{i}.de - mean(apa_mass{i}.de(apa_mass{i}.t<11));
end
#+end_src

The raw measurements are shown in Figure [[fig:apa_meas_k_time]].
All the APA seems to have similar stiffness except the APA 7 which should have an higher stiffness.

#+begin_question
It is however strange that the displacement $d_e$ when the mass is removed is higher for the APA 7 than for the other APA.
What could cause that?
#+end_question

#+begin_src matlab :exports none
figure;
hold on;
for i = 1:length(apa_nums)
  plot(apa_mass{i}.t, apa_mass{i}.de, 'DisplayName', sprintf('APA %i', apa_nums(i)));
end
hold off;
xlabel('Time [s]'); ylabel('Displacement $d_e$ [m]');
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2);
#+end_src

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

#+name: fig:apa_meas_k_time
#+caption: Raw measurements for all the APA. A mass of 6.4kg is added at arround 15s and removed at arround 22s
#+RESULTS:
[[file:figs/apa_meas_k_time.png]]

#+begin_src matlab
added_mass = 6.4; % Added mass [kg]
#+end_src

The stiffnesses are computed for all the APA and are summarized in Table [[tab:APA_measured_k]].

#+begin_src matlab :exports none
apa_k = zeros(length(apa_nums), 1);
for i = 1:length(apa_nums)
    apa_k(i) = 9.8 * added_mass / (mean(apa_mass{i}.de(apa_mass{i}.t > 12 & apa_mass{i}.t < 12.5)) - mean(apa_mass{i}.de(apa_mass{i}.t > 20 & apa_mass{i}.t < 20.5)));
end
#+end_src

#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
  data2orgtable(1e-6*apa_k, cellstr(num2str(apa_nums')), {'APA Num', '$k [N/\mu m]$'}, ' %.2f ');
#+end_src

#+name: tab:APA_measured_k
#+caption: Measured stiffnesses
#+attr_latex: :environment tabularx :width 0.3\linewidth :align cc
#+attr_latex: :center t :booktabs t :float t
#+RESULTS:
| APA Num | $k [N/\mu m]$ |
|---------+---------------|
|       1 |          1.68 |
|       2 |          1.69 |
|       4 |           1.7 |
|       5 |           1.7 |
|       6 |           1.7 |
|       7 |          1.93 |
|       8 |          1.73 |

#+begin_important
The APA300ML manual specifies the nominal stiffness to be $1.8\,[N/\mu m]$ which is very close to what have been measured.
Only the APA number 7 is a little bit off.
#+end_important

*** Stiffness - Effect of connecting the actuator stack to the amplifier and the sensor stack to the ADC
We wish here to see if the stiffness changes when the actuator stacks are not connected to the amplifier and the sensor stacks are not connected to the ADC.

Note here that the resistor in parallel to the sensor stack is present in both cases.

First, the data are loaded.
#+begin_src matlab
add_mass_oc = load(sprintf('frf_data_%i_add_mass_open_circuit.mat',   1), 't', 'de');
add_mass_cc = load(sprintf('frf_data_%i_add_mass_closed_circuit.mat', 1), 't', 'de');
#+end_src

And the initial displacement is set to zero.
#+begin_src matlab
add_mass_oc.de = add_mass_oc.de - mean(add_mass_oc.de(add_mass_oc.t<11));
add_mass_cc.de = add_mass_cc.de - mean(add_mass_cc.de(add_mass_cc.t<11));
#+end_src

#+begin_src matlab :exports none
figure;
hold on;
plot(add_mass_oc.t, add_mass_oc.de, 'DisplayName', 'Not connected');
plot(add_mass_cc.t, add_mass_cc.de, 'DisplayName', 'Connected');
hold off;
xlabel('Time [s]'); ylabel('Displacement $d_e$ [m]');
legend('location', 'northeast');
#+end_src

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

#+name: fig:apa_meas_k_time_oc_cc
#+caption: Measured displacement
#+RESULTS:
[[file:figs/apa_meas_k_time_oc_cc.png]]


And the stiffness is estimated in both case.
The results are shown in Table [[tab:APA_measured_k_oc_cc]].

#+begin_src matlab
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_cc = 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)));
#+end_src

#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
data2orgtable(1e-6*[apa_k_oc; apa_k_cc], {'Not connected', 'Connected'}, {'$k [N/\mu m]$'}, ' %.1f ');
#+end_src

#+name: tab:APA_measured_k_oc_cc
#+caption: Measured stiffnesses on "open" and "closed" circuits
#+attr_latex: :environment tabularx :width 0.3\linewidth :align cc
#+attr_latex: :center t :booktabs t :float t
#+RESULTS:
|               | $k [N/\mu m]$ |
|---------------+---------------|
| Not connected |           2.3 |
| Connected     |           1.7 |

#+begin_important
Clearly, connecting the actuator stacks to the amplified (basically equivalent as to short circuiting them) lowers the stiffness.
#+end_important

*** Hysteresis
We here wish to visually see the amount of hysteresis present in the APA.

The data is loaded.
#+begin_src matlab
apa_hyst = load('frf_data_1_hysteresis.mat', 't', 'Va', 'de');
% Initial time set to zero
apa_hyst.t = apa_hyst.t - apa_hyst.t(1);
#+end_src

The excitation voltage amplitudes are:
#+begin_src matlab
ampls = [0.1, 0.2, 0.4, 1, 2, 4]; % Excitation voltage amplitudes
#+end_src

The excitation voltage and the measured displacement are shown in Figure [[fig:hyst_exc_signal_time]].
#+begin_src matlab :exports none
figure;
tiledlayout(1, 2, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
plot(apa_hyst.t, apa_hyst.Va)
xlabel('Time [s]'); ylabel('Output Voltage [V]');

ax2 = nexttile;
plot(apa_hyst.t, apa_hyst.de)
xlabel('Time [s]'); ylabel('Measured Displacement [m]');
#+end_src

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

#+name: fig:hyst_exc_signal_time
#+caption: Excitation voltage and measured displacement
#+RESULTS:
[[file:figs/hyst_exc_signal_time.png]]

For each amplitude, we only take the last sinus in order to reduce possible transients.
Also, it is centered on zero.

The measured displacement at a function of the output voltage are shown in Figure [[fig:hyst_results_multi_ampl]].
#+begin_src matlab :exports none
figure;
tiledlayout(1, 3, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile([1,2]);
hold on;
for i = flip(1:6)
    i_lim = apa_hyst.t > i*5-1 & apa_hyst.t < i*5;
    plot(apa_hyst.Va(i_lim) - mean(apa_hyst.Va(i_lim)), apa_hyst.de(i_lim) - mean(apa_hyst.de(i_lim)), ...
         'DisplayName', sprintf('$V_a = %.1f [V]$', ampls(i)))
end
hold off;
xlabel('Output Voltage [V]'); ylabel('Measured Displacement [m]');
legend('location', 'northeast');
xlim([-4, 4]); ylim([-1.2e-4, 1.2e-4]);

ax2 = nexttile;
hold on;
for i = flip(1:6)
    i_lim = apa_hyst.t > i*5-1 & apa_hyst.t < i*5;
    plot(apa_hyst.Va(i_lim) - mean(apa_hyst.Va(i_lim)), apa_hyst.de(i_lim) - mean(apa_hyst.de(i_lim)))
end
hold off;
xlim([-0.4, 0.4]); ylim([-0.8e-5, 0.8e-5]);
#+end_src

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

#+name: fig:hyst_results_multi_ampl
#+caption: Obtained hysteresis for multiple excitation amplitudes
#+RESULTS:
[[file:figs/hyst_results_multi_ampl.png]]

#+begin_important
It is quite clear that hysteresis is increasing with the excitation amplitude.

Also, no hysteresis is found on the sensor stack voltage.
#+end_important

*** FRF Identification - Setup
The identification is performed in three steps:
1. White noise excitation with small amplitude.
   This is used to determine the main resonance of the system.
2. Sweep sine excitation with the amplitude lowered around the resonance.
   The sweep sine is from 10Hz to 400Hz.
3. High frequency noise.
   The noise is band-passed between 300Hz and 2kHz.

Then, the result of the second identification is used between 10Hz and 350Hz and the result of the third identification if used between 350Hz and 2kHz.

Here are the APA numbers that have been measured.
#+begin_src matlab
apa_nums = [1 2 4 5 6 7 8];
#+end_src

The data are loaded for both the second and third identification:
#+begin_src matlab
%% Second identification
apa_sweep = {};
for i = 1:length(apa_nums)
    apa_sweep(i) = {load(sprintf('frf_data_%i_sweep.mat', apa_nums(i)), 't', 'Va', 'Vs', 'de', 'da')};
end

%% Third identification
apa_noise_hf = {};
for i = 1:length(apa_nums)
    apa_noise_hf(i) = {load(sprintf('frf_data_%i_noise_hf.mat', apa_nums(i)), 't', 'Va', 'Vs', 'de', 'da')};
end
#+end_src

The time is the same for all measurements.
#+begin_src matlab
%% Time vector
t = apa_sweep{1}.t - apa_sweep{1}.t(1) ; % Time vector [s]

%% Sampling
Ts = (t(end) - t(1))/(length(t)-1); % Sampling Time [s]
Fs = 1/Ts; % Sampling Frequency [Hz]
#+end_src

Then we defined a "Hanning" windows that will be used for the spectral analysis:
#+begin_src matlab
win = hanning(ceil(0.5*Fs)); % Hannning Windows
#+end_src

We get the frequency vector that will be the same for all the frequency domain analysis.
#+begin_src matlab
% Only used to have the frequency vector "f"
[~, f] = tfestimate(apa_sweep{1}.Va, apa_sweep{1}.de, win, [], [], 1/Ts);
#+end_src

*** FRF Identification - DVF
In this section, the dynamics from $V_a$ to $d_e$ is identified.

We compute the coherence for 2nd and 3rd identification:
#+begin_src matlab
%% Coherence computation
coh_sweep = zeros(length(f), length(apa_nums));
for i = 1:length(apa_nums)
    [coh, ~] = mscohere(apa_sweep{i}.Va, apa_sweep{i}.de, win, [], [], 1/Ts);
    coh_sweep(:, i) = coh;
end

coh_noise_hf = zeros(length(f), length(apa_nums));
for i = 1:length(apa_nums)
    [coh, ~] = mscohere(apa_noise_hf{i}.Va, apa_noise_hf{i}.de, win, [], [], 1/Ts);
    coh_noise_hf(:, i) = coh;
end
#+end_src

The coherence is shown in Figure [[fig:frf_dvf_plant_coh]].
It is clear that the Sweep sine gives good coherence up to 400Hz and that the high frequency noise excitation signal helps increasing a little bit the coherence at high frequency.

#+begin_src matlab :exports none
colors = get(gca,'colororder');

figure;
hold on;
plot(f, coh_noise_hf(:, 1), 'color', [colors(1, :), 0.5], ...
     'DisplayName', 'HF Noise');
plot(f, coh_sweep(:, 1), 'color', [colors(2, :), 0.5], ...
     'DisplayName', 'Sweep');
for i = 2:length(apa_nums)
    plot(f, coh_noise_hf(:, i), 'color', [colors(1, :), 0.5], ...
         'HandleVisibility', 'off');
    plot(f, coh_sweep(:, i), 'color', [colors(2, :), 0.5], ...
         'HandleVisibility', 'off');
end;
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Coherence [-]');
xlim([5, 5e3]); ylim([0, 1]);
legend('location', 'southeast');
#+end_src

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

#+name: fig:frf_dvf_plant_coh
#+caption: Obtained coherence for the plant from $V_a$ to $d_e$
#+RESULTS:
[[file:figs/frf_dvf_plant_coh.png]]


Then, the transfer function from the DAC output voltage $V_a$ to the measured displacement by the encoders is computed:
#+begin_src matlab
%% Transfer function estimation
dvf_sweep = zeros(length(f), length(apa_nums));
for i = 1:length(apa_nums)
    [frf, ~] = tfestimate(apa_sweep{i}.Va, apa_sweep{i}.de, win, [], [], 1/Ts);
    dvf_sweep(:, i) = frf;
end

dvf_noise_hf = zeros(length(f), length(apa_nums));
for i = 1:length(apa_nums)
    [frf, ~] = tfestimate(apa_noise_hf{i}.Va, apa_noise_hf{i}.de, win, [], [], 1/Ts);
    dvf_noise_hf(:, i) = frf;
end
#+end_src

The obtained transfer functions are shown in Figure [[fig:frf_dvf_plant_tf]].

They are all superimposed except for the APA7.

#+begin_question
Why is the APA7 off?
We could think that the APA7 is stiffer, but also the mass line is off.

It seems that there is a "gain" problem.
The encoder seems fine (it measured the same as the Interferometer).
Maybe it could be due to the amplifier?
#+end_question

#+begin_question
Why is there a double resonance at around 94Hz?
#+end_question

#+begin_src matlab :exports none
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
for i = 1:length(apa_nums)
    plot(f(f> 350), abs(dvf_noise_hf(f> 350, i)), 'color', colors(i, :), ...
         'DisplayName', sprintf('APA %i', apa_nums(i)));
    plot(f(f<=350), abs(dvf_sweep(   f<=350, i)), 'color', colors(i, :), ...
         'HandleVisibility', 'off');
end
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d_e/V_a$ [m/V]'); set(gca, 'XTickLabel',[]);
hold off;
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2);
ylim([1e-9, 1e-3]);

ax2 = nexttile;
hold on;
for i = 1:length(apa_nums)
    plot(f(f> 350), 180/pi*angle(dvf_noise_hf(f> 350, i)), 'color', colors(i, :));
    plot(f(f<=350), 180/pi*angle(dvf_sweep(   f<=350, i)), 'color', colors(i, :));
end
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

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

#+name: fig:frf_dvf_plant_tf
#+caption: Estimated FRF for the DVF plant (transfer function from $V_a$ to the encoder $d_e$)
#+RESULTS:
[[file:figs/frf_dvf_plant_tf.png]]


*** FRF Identification - IFF
In this section, the dynamics from $V_a$ to $V_s$ is identified.

First the coherence is computed and shown in Figure [[fig:frf_iff_plant_coh]].
The coherence is very nice from 10Hz to 2kHz.
It is only dropping near a zeros at 40Hz, and near the resonance at 95Hz (the excitation amplitude being lowered).

#+begin_src matlab
%% Coherence
coh_sweep = zeros(length(f), length(apa_nums));
for i = 1:length(apa_nums)
    [coh, ~] = mscohere(apa_sweep{i}.Va, apa_sweep{i}.Vs, win, [], [], 1/Ts);
    coh_sweep(:, i) = coh;
end

coh_noise_hf = zeros(length(f), length(apa_nums));
for i = 1:length(apa_nums)
    [coh, ~] = mscohere(apa_noise_hf{i}.Va, apa_noise_hf{i}.Vs, win, [], [], 1/Ts);
    coh_noise_hf(:, i) = coh;
end
#+end_src

#+begin_src matlab :exports none
colors = get(gca,'colororder');

figure;
hold on;
plot(f, coh_noise_hf(:, 1), 'color', [colors(1, :), 0.5], 'DisplayName', 'HF Noise');
plot(f, coh_sweep(   :, 1), 'color', [colors(2, :), 0.5], 'DisplayName', 'Sweep');
for i = 2:length(apa_nums)
    plot(f, coh_noise_hf(:, i), 'color', [colors(1, :), 0.5], ...
         'HandleVisibility', 'off');
    plot(f, coh_sweep(   :, i), 'color', [colors(2, :), 0.5], ...
         'HandleVisibility', 'off');
end;
hold off;
xlabel('Frequency [Hz]'); ylabel('Coherence [-]');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlim([5, 5e3]); ylim([0, 1]);
legend('location', 'southeast');
#+end_src

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

#+name: fig:frf_iff_plant_coh
#+caption: Obtained coherence for the IFF plant
#+RESULTS:
[[file:figs/frf_iff_plant_coh.png]]

Then the FRF are estimated and shown in Figure [[fig:frf_iff_plant_tf]]
#+begin_src matlab
%% FRF estimation of the transfer function from Va to Vs
iff_sweep = zeros(length(f), length(apa_nums));
for i = 1:length(apa_nums)
    [frf, ~] = tfestimate(apa_sweep{i}.Va, apa_sweep{i}.Vs, win, [], [], 1/Ts);
    iff_sweep(:, i) = frf;
end

iff_noise_hf = zeros(length(f), length(apa_nums));
for i = 1:length(apa_nums)
    [frf, ~] = tfestimate(apa_noise_hf{i}.Va, apa_noise_hf{i}.Vs, win, [], [], 1/Ts);
    iff_noise_hf(:, i) = frf;
end
#+end_src

#+begin_src matlab :exports none
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
for i = 1:length(apa_nums)
    plot(f(f> 350), abs(iff_noise_hf(f> 350, i)), 'color', colors(i, :), ...
         'DisplayName', sprintf('APA %i', apa_nums(i)));
    plot(f(f<=350), abs(iff_sweep(   f<=350, i)), 'color', colors(i, :), ...
         'HandleVisibility', 'off');
end
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', 'FontSize', 8, 'NumColumns', 2);

ax2 = nexttile;
hold on;
for i = 1:length(apa_nums)
    plot(f(f> 350), 180/pi*angle(iff_noise_hf(f> 350, i)), 'color', colors(i, :));
    plot(ff<=350), 180/pi*angle(iff_sweep(   f<=350, i)), 'color', colors(i, :));
end
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

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

#+name: fig:frf_iff_plant_tf
#+caption:Identified IFF Plant
#+RESULTS:
[[file:figs/frf_iff_plant_tf.png]]

*** Effect of the resistor on the IFF Plant
A resistor is added in parallel with the sensor stack.
This has the effect to form a high pass filter with the capacitance of the stack.

We here measured the low frequency transfer function from $V_a$ to $V_s$ with and without this resistor.

#+begin_src matlab
% With the resistor
wi_k = load('frf_data_2_sweep_lf_with_R.mat', 't', 'Vs', 'Va');

% Without the resistor
wo_k = load('frf_data_2_sweep_lf.mat',        't', 'Vs', 'Va');
#+end_src

#+begin_src matlab
t = wo_k.t; % Time vector [s]

Ts = (t(end) - t(1))/(length(t)-1); % Sampling Time [s]
Fs = 1/Ts; % Sampling Frequency [Hz]
#+end_src

We use a very long "Hanning" window for the spectral analysis in order to estimate the low frequency behavior.
#+begin_src matlab
win = hanning(ceil(50*Fs)); % Hannning Windows
#+end_src

And we estimate the transfer function from $V_a$ to $V_s$ in both cases:
#+begin_src matlab
[frf_wo_k, f] = tfestimate(wo_k.Va, wo_k.Vs, win, [], [], 1/Ts);
[frf_wi_k, ~] = tfestimate(wi_k.Va, wi_k.Vs, win, [], [], 1/Ts);
#+end_src


#+begin_src matlab
f0 = 0.35;
G_hpf = 0.6*(s/2*pi*f0)/(1 + s/2*pi*f0);
#+end_src

#+begin_src matlab :exports none
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f, abs(frf_wo_k));
plot(f, abs(frf_wi_k));
plot(f, abs(squeeze(freqresp(G_hpf, f, 'Hz'))), 'k--');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $V_{out}/V_{in}$ [V/V]'); set(gca, 'XTickLabel',[]);
hold off;
ylim([1e-3, 1e1]);

ax2 = nexttile;
hold on;
plot(f, 180/pi*angle(frf_wo_k), 'DisplayName', 'Without $k$');
plot(f, 180/pi*angle(frf_wi_k), 'DisplayName', 'With $k$');
plot(f, 180/pi*angle(squeeze(freqresp(G_hpf, f, 'Hz'))), 'k--', 'DisplayName', sprintf('HPF $f_o = %.2f [Hz]$', f0));
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]);
legend('location', 'northeast')

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

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

#+name: fig:frf_iff_effect_R
#+caption: Transfer function from $V_a$ to $V_s$ with and without the resistor $k$
#+RESULTS:
[[file:figs/frf_iff_effect_R.png]]


** TODO Measurement on Strut 1
*** Introduction                                                    :ignore:
Measurements are first performed on the strut number 1 with:
- APA 1
- flex 1 and flex 2

*** 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/');
addpath('./matlab/src/');
addpath('./matlab/');
#+end_src

#+begin_src matlab :eval no
addpath('./mat/');
addpath('./src/');
#+end_src

*** Without Encoder
**** FRF Identification - Setup
The identification is performed in three steps:
1. White noise excitation with small amplitude.
   This is used to determine the main resonance of the system.
2. Sweep sine excitation with the amplitude lowered around the resonance.
   The sweep sine is from 10Hz to 400Hz.
3. High frequency noise.
   The noise is band-passed between 300Hz and 2kHz.

Then, the result of the second identification is used between 10Hz and 350Hz and the result of the third identification if used between 350Hz and 2kHz.

#+begin_src matlab
leg_sweep    = load(sprintf('frf_data_leg_%i_sweep.mat',    1), 't', 'Va', 'Vs', 'de', 'da');
leg_noise_hf = load(sprintf('frf_data_leg_%i_noise_hf.mat', 1), 't', 'Va', 'Vs', 'de', 'da');
#+end_src

The time is the same for all measurements.
#+begin_src matlab
%% Time vector
t = leg_sweep.t - leg_sweep.t(1) ; % Time vector [s]

%% Sampling
Ts = (t(end) - t(1))/(length(t)-1); % Sampling Time [s]
Fs = 1/Ts; % Sampling Frequency [Hz]
#+end_src

Then we defined a "Hanning" windows that will be used for the spectral analysis:
#+begin_src matlab
win = hanning(ceil(0.5*Fs)); % Hannning Windows
#+end_src

We get the frequency vector that will be the same for all the frequency domain analysis.
#+begin_src matlab
% Only used to have the frequency vector "f"
[~, f] = tfestimate(leg_sweep.Va, leg_sweep.de, win, [], [], 1/Ts);
#+end_src

**** TODO FRF Identification - DVF
In this section, the dynamics from $V_a$ to $d_e$ is identified.

We compute the coherence for 2nd and 3rd identification:
#+begin_src matlab
[coh_sweep, ~]    = mscohere(leg_sweep.Va,    leg_sweep.da,    win, [], [], 1/Ts);
[coh_noise_hf, ~] = mscohere(leg_noise_hf.Va, leg_noise_hf.da, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
colors = get(gca,'colororder');

figure;
hold on;
plot(f, coh_noise_hf(:, 1), 'color', [colors(1, :), 0.5], ...
     'DisplayName', 'HF Noise');
plot(f, coh_sweep(:, 1), 'color', [colors(2, :), 0.5], ...
     'DisplayName', 'Sweep');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Coherence [-]');
% xlim([5, 5e3]); ylim([0, 1]);
legend('location', 'southeast');
#+end_src

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

#+name: fig:frf_dvf_plant_coh
#+caption: Obtained coherence for the plant from $V_a$ to $d_e$
#+RESULTS:
[[file:figs/frf_dvf_plant_coh.png]]


#+begin_src matlab
[dvf_sweep, ~]    = tfestimate(leg_sweep.Va,    leg_sweep.da,    win, [], [], 1/Ts);
[dvf_noise_hf, ~] = tfestimate(leg_noise_hf.Va, leg_noise_hf.da, win, [], [], 1/Ts);
#+end_src

The obtained transfer functions are shown in Figure [[fig:frf_dvf_plant_tf]].

They are all superimposed except for the APA7.

#+begin_question
Why is the APA7 off?
We could think that the APA7 is stiffer, but also the mass line is off.

It seems that there is a "gain" problem.
The encoder seems fine (it measured the same as the Interferometer).
Maybe it could be due to the amplifier?
#+end_question

#+begin_question
Why is there a double resonance at around 94Hz?
#+end_question

#+begin_src matlab :exports none
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f(f> 350), abs(dvf_noise_hf(f> 350)), 'color', colors(1, :));
plot(f(f<=350), abs(dvf_sweep(   f<=350)), 'color', colors(1, :));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d_e/V_a$ [m/V]'); set(gca, 'XTickLabel',[]);
hold off;
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2);
ylim([1e-9, 1e-3]);

ax2 = nexttile;
hold on;
plot(f(f> 350), 180/pi*angle(dvf_noise_hf(f> 350)), 'color', colors(1, :));
plot(f(f<=350), 180/pi*angle(dvf_sweep(   f<=350)), 'color', colors(1, :));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

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

#+name: fig:frf_dvf_plant_tf
#+caption: Estimated FRF for the DVF plant (transfer function from $V_a$ to the encoder $d_e$)
#+RESULTS:
[[file:figs/frf_dvf_plant_tf.png]]

**** TODO FRF Identification - IFF
In this section, the dynamics from $V_a$ to $V_s$ is identified.

First the coherence is computed and shown in Figure [[fig:frf_iff_plant_coh]].
The coherence is very nice from 10Hz to 2kHz.
It is only dropping near a zeros at 40Hz, and near the resonance at 95Hz (the excitation amplitude being lowered).

#+begin_src matlab
[coh_sweep, ~] = mscohere(leg_sweep.Va, leg_sweep.Vs, win, [], [], 1/Ts);
[coh_noise_hf, ~] = mscohere(leg_noise_hf.Va, leg_noise_hf.Vs, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
colors = get(gca,'colororder');

figure;
hold on;
plot(f, coh_noise_hf, 'color', [colors(1, :), 0.5], 'DisplayName', 'HF Noise');
plot(f, coh_sweep,    'color', [colors(2, :), 0.5], 'DisplayName', 'Sweep');
hold off;
xlabel('Frequency [Hz]'); ylabel('Coherence [-]');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlim([5, 5e3]); ylim([0, 1]);
legend('location', 'southeast');
#+end_src

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

#+name: fig:frf_iff_plant_coh
#+caption: Obtained coherence for the IFF plant
#+RESULTS:
[[file:figs/frf_iff_plant_coh.png]]

Then the FRF are estimated and shown in Figure [[fig:frf_iff_plant_tf]]
#+begin_src matlab
[iff_sweep, ~] = tfestimate(leg_sweep.Va, leg_sweep.Vs, win, [], [], 1/Ts);
[iff_noise_hf, ~] = tfestimate(leg_noise_hf.Va, leg_noise_hf.Vs, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f(f> 350), abs(iff_noise_hf(f> 350)), 'color', colors(1, :));
plot(f(f<=350), abs(iff_sweep(   f<=350)), 'color', colors(1, :));
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', 'FontSize', 8, 'NumColumns', 2);

ax2 = nexttile;
hold on;
plot(f(f> 350), 180/pi*angle(iff_noise_hf(f> 350)), 'color', colors(1, :));
plot(f(f<=350), 180/pi*angle(iff_sweep(   f<=350)), 'color', colors(1, :));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

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

#+name: fig:frf_iff_plant_tf
#+caption:Identified IFF Plant
#+RESULTS:
[[file:figs/frf_iff_plant_tf.png]]

*** With Encoder
**** FRF Identification - Setup
The identification is performed in three steps:
1. White noise excitation with small amplitude.
   This is used to determine the main resonance of the system.
2. Sweep sine excitation with the amplitude lowered around the resonance.
   The sweep sine is from 10Hz to 400Hz.
3. High frequency noise.
   The noise is band-passed between 300Hz and 2kHz.

Then, the result of the second identification is used between 10Hz and 350Hz and the result of the third identification if used between 350Hz and 2kHz.

#+begin_src matlab
leg_enc_sweep    = load(sprintf('frf_data_leg_coder_%i_noise.mat',    1), 't', 'Va', 'Vs', 'de', 'da');
leg_enc_noise_hf = load(sprintf('frf_data_leg_coder_%i_noise_hf.mat', 1), 't', 'Va', 'Vs', 'de', 'da');
#+end_src

**** TODO FRF Identification - DVF
In this section, the dynamics from $V_a$ to $d_e$ is identified.

We compute the coherence for 2nd and 3rd identification:
#+begin_src matlab
[coh_enc_sweep, ~]    = mscohere(leg_enc_sweep.Va,    leg_enc_sweep.de,    win, [], [], 1/Ts);
[coh_enc_noise_hf, ~] = mscohere(leg_enc_noise_hf.Va, leg_enc_noise_hf.de, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
colors = get(gca,'colororder');

figure;
hold on;
plot(f, coh_enc_noise_hf(:, 1), 'color', [colors(1, :), 0.5], ...
     'DisplayName', 'HF Noise');
plot(f, coh_enc_sweep(:, 1), 'color', [colors(2, :), 0.5], ...
     'DisplayName', 'Sweep');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Coherence [-]');
% xlim([5, 5e3]); ylim([0, 1]);
legend('location', 'southeast');
#+end_src

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

#+name: fig:frf_dvf_plant_coh
#+caption: Obtained coherence for the plant from $V_a$ to $d_e$
#+RESULTS:
[[file:figs/frf_dvf_plant_coh.png]]


#+begin_src matlab
[dvf_enc_sweep, ~]    = tfestimate(leg_enc_sweep.Va,    leg_enc_sweep.de,    win, [], [], 1/Ts);
[dvf_enc_noise_hf, ~] = tfestimate(leg_enc_noise_hf.Va, leg_enc_noise_hf.de, win, [], [], 1/Ts);
#+end_src

The obtained transfer functions are shown in Figure [[fig:frf_dvf_plant_tf]].

They are all superimposed except for the APA7.

#+begin_question
Why is the APA7 off?
We could think that the APA7 is stiffer, but also the mass line is off.

It seems that there is a "gain" problem.
The encoder seems fine (it measured the same as the Interferometer).
Maybe it could be due to the amplifier?
#+end_question

#+begin_question
Why is there a double resonance at around 94Hz?
#+end_question

#+begin_src matlab :exports none
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f(f> 350), abs(dvf_enc_noise_hf(f> 350)), 'color', colors(1, :));
plot(f(f<=350), abs(dvf_enc_sweep(   f<=350)), 'color', colors(1, :));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d_e/V_a$ [m/V]'); set(gca, 'XTickLabel',[]);
hold off;
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2);
ylim([1e-9, 1e-3]);

ax2 = nexttile;
hold on;
plot(f(f> 350), 180/pi*angle(dvf_enc_noise_hf(f> 350)), 'color', colors(1, :));
plot(f(f<=350), 180/pi*angle(dvf_enc_sweep(   f<=350)), 'color', colors(1, :));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

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

#+name: fig:frf_dvf_plant_tf
#+caption: Estimated FRF for the DVF plant (transfer function from $V_a$ to the encoder $d_e$)
#+RESULTS:
[[file:figs/frf_dvf_plant_tf.png]]

**** Comparison with Interferometer

#+begin_src matlab
[dvf_int_sweep, ~]    = tfestimate(leg_enc_sweep.Va,    leg_enc_sweep.da,    win, [], [], 1/Ts);
[dvf_int_noise_hf, ~] = tfestimate(leg_enc_noise_hf.Va, leg_enc_noise_hf.da, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f(f> 350), abs(dvf_enc_noise_hf(f> 350)), 'color', colors(1, :), ...
     'DisplayName', 'Encoder');
plot(f(f<=350), abs(dvf_enc_sweep(   f<=350)), 'color', colors(1, :), ...
     'HandleVisibility', 'off');
plot(f(f> 350), abs(dvf_int_noise_hf(f> 350)), 'color', colors(2, :), ...
     'DisplayName', 'Interferometer');
plot(f(f<=350), abs(dvf_int_sweep(   f<=350)), 'color', colors(2, :), ...
     'HandleVisibility', 'off');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d_e/V_a$ [m/V]'); set(gca, 'XTickLabel',[]);
hold off;
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2);
ylim([1e-9, 1e-3]);

ax2 = nexttile;
hold on;
plot(f(f> 350), 180/pi*angle(dvf_enc_noise_hf(f> 350)), 'color', colors(1, :));
plot(f(f<=350), 180/pi*angle(dvf_enc_sweep(   f<=350)), 'color', colors(1, :));
plot(f(f> 350), 180/pi*angle(dvf_int_noise_hf(f> 350)), 'color', colors(1, :));
plot(f(f<=350), 180/pi*angle(dvf_int_sweep(   f<=350)), 'color', colors(1, :));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

#+begin_important
Clearly using the encoder like this will not be useful.

Probably, the encoders will have to be fixed on the plates so that the resonances of the APA are not a problem anymore.
#+end_important

**** APA Resonances Frequency

#+begin_src matlab :exports none
figure;
hold on;
plot(f(f> 350), abs(dvf_enc_noise_hf(f> 350)), 'k-');
plot(f(f<=350), abs(dvf_enc_sweep(   f<=350)), 'k-');

text(93,  4e-4,  {'93Hz'}, 'VerticalAlignment','bottom','HorizontalAlignment','center')
text(197, 1.3e-4,{'197Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
text(290, 4e-6,  {'290Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
text(376, 1.4e-6,{'376Hz'},'VerticalAlignment','bottom','HorizontalAlignment','center')
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-7, 1e-3]); xlim([10, 2e3]);
#+end_src

This is very close to what was estimated using the FEM.

#+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]]

#+begin_important
The resonances are indeed due to limited stiffness of the APA.
#+end_important

**** FRF Identification - IFF
In this section, the dynamics from $V_a$ to $V_s$ is identified.

First the coherence is computed and shown in Figure [[fig:frf_iff_plant_coh]].
The coherence is very nice from 10Hz to 2kHz.
It is only dropping near a zeros at 40Hz, and near the resonance at 95Hz (the excitation amplitude being lowered).

#+begin_src matlab
[coh_enc_sweep,    ~] = mscohere(leg_enc_sweep.Va,    leg_enc_sweep.Vs,    win, [], [], 1/Ts);
[coh_enc_noise_hf, ~] = mscohere(leg_enc_noise_hf.Va, leg_enc_noise_hf.Vs, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
colors = get(gca,'colororder');

figure;
hold on;
plot(f, coh_enc_noise_hf, 'color', [colors(1, :), 0.5], 'DisplayName', 'HF Noise');
plot(f, coh_enc_sweep,    'color', [colors(2, :), 0.5], 'DisplayName', 'Sweep');
hold off;
xlabel('Frequency [Hz]'); ylabel('Coherence [-]');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlim([5, 5e3]); ylim([0, 1]);
legend('location', 'southeast');
#+end_src

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

#+name: fig:frf_iff_plant_coh
#+caption: Obtained coherence for the IFF plant
#+RESULTS:
[[file:figs/frf_iff_plant_coh.png]]

Then the FRF are estimated and shown in Figure [[fig:frf_iff_plant_tf]]
#+begin_src matlab
[iff_enc_sweep,    ~] = tfestimate(leg_enc_sweep.Va,    leg_enc_sweep.Vs,    win, [], [], 1/Ts);
[iff_enc_noise_hf, ~] = tfestimate(leg_enc_noise_hf.Va, leg_enc_noise_hf.Vs, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f(f> 350), abs(iff_enc_noise_hf(f> 350)), 'color', colors(1, :));
plot(f(f<=350), abs(iff_enc_sweep(   f<=350)), 'color', colors(1, :));
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', 'FontSize', 8, 'NumColumns', 2);

ax2 = nexttile;
hold on;
plot(f(f> 350), 180/pi*angle(iff_enc_noise_hf(f> 350)), 'color', colors(1, :));
plot(f(f<=350), 180/pi*angle(iff_enc_sweep(   f<=350)), 'color', colors(1, :));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

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

#+name: fig:frf_iff_plant_tf
#+caption:Identified IFF Plant
#+RESULTS:
[[file:figs/frf_iff_plant_tf.png]]

**** Comparison to when the encoder is not fixed

#+begin_src matlab
[iff_sweep, ~]    = tfestimate(leg_sweep.Va,    leg_sweep.Vs,    win, [], [], 1/Ts);
[iff_noise_hf, ~] = tfestimate(leg_noise_hf.Va, leg_noise_hf.Vs, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(f(f> 350), abs(iff_enc_noise_hf(f> 350)), 'color', colors(1, :), ...
     'DisplayName', 'Encoder');
plot(f(f<=350), abs(iff_enc_sweep(   f<=350)), 'color', colors(1, :), ...
     'HandleVisibility', 'off');
plot(f(f> 350), abs(iff_noise_hf(f> 350)), 'color', colors(2, :), ...
     'DisplayName', 'Withotu Encoder');
plot(f(f<=350), abs(iff_sweep(   f<=350)), 'color', colors(2, :), ...
     'HandleVisibility', 'off');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d_e/V_a$ [m/V]'); set(gca, 'XTickLabel',[]);
hold off;
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2);
ylim([1e-2, 1e2]);

ax2 = nexttile;
hold on;
plot(f(f> 350), 180/pi*angle(iff_enc_noise_hf(f> 350)), 'color', colors(1, :));
plot(f(f<=350), 180/pi*angle(iff_enc_sweep(   f<=350)), 'color', colors(1, :));
plot(f(f> 350), 180/pi*angle(iff_noise_hf(f> 350)), 'color', colors(2, :));
plot(f(f<=350), 180/pi*angle(iff_sweep(   f<=350)), 'color', colors(2, :));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);

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

#+begin_important
We can see that the IFF does not change whether of not the encoder are fixed to the struts.
#+end_important

* Measurement Results

* TODO Compare with the FEM/Simscape Model                          :noexport:
:PROPERTIES:
:header-args:matlab+: :tangle matlab/APA300ML.m
:END:

** Introduction                                                      :ignore:
In this section, the Amplified Piezoelectric Actuator APA300ML ([[file:doc/APA300ML.pdf][doc]]) is modeled using a Finite Element Software.
Then a /super element/ is exported and imported in Simscape where its dynamic is studied.

A 3D view of the Amplified Piezoelectric Actuator (APA300ML) is shown in Figure [[fig:apa300ml_ansys]].
The remote point used are also shown in this figure.

#+name: fig:apa300ml_ansys
#+caption: Ansys FEM of the APA300ML
[[file:figs/apa300ml_ansys.jpg]]

** 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/APA300ML/');
#+end_src

#+begin_src matlab :eval no
  addpath('APA300ML/');
#+end_src

#+begin_src matlab
  open('APA300ML.slx');
#+end_src

** Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates
We first extract the stiffness and mass matrices.
#+begin_src matlab
  K = readmatrix('APA300ML_mat_K.CSV');
  M = readmatrix('APA300ML_mat_M.CSV');
#+end_src

#+begin_src matlab :exports results :results value table replace :tangle no
  data2orgtable(K(1:10, 1:10), {}, {}, ' %.1g ');
#+end_src

#+caption: First 10x10 elements of the Stiffness matrix
#+RESULTS:
| 200000000.0 |    30000.0 |  -20000.0 |     -70.0 | 300000.0 |    40.0 |  10000000.0 |    10000.0 |   -6000.0 |     30.0 |
|     30000.0 | 30000000.0 |    2000.0 | -200000.0 |     60.0 |   -10.0 |      4000.0 |  2000000.0 |    -500.0 |   9000.0 |
|    -20000.0 |     2000.0 | 7000000.0 |     -10.0 |    -30.0 |    10.0 |      6000.0 |      900.0 | -500000.0 |        3 |
|       -70.0 |  -200000.0 |     -10.0 |    1000.0 |     -0.1 |    0.08 |       -20.0 |    -9000.0 |         3 |    -30.0 |
|    300000.0 |       60.0 |     -30.0 |      -0.1 |    900.0 |     0.1 |     30000.0 |       20.0 |     -10.0 |     0.06 |
|        40.0 |      -10.0 |      10.0 |      0.08 |      0.1 | 10000.0 |        20.0 |          9 |        -5 |     0.03 |
|  10000000.0 |     4000.0 |    6000.0 |     -20.0 |  30000.0 |    20.0 | 200000000.0 |    10000.0 |    9000.0 |     50.0 |
|     10000.0 |  2000000.0 |     900.0 |   -9000.0 |     20.0 |       9 |     10000.0 | 30000000.0 |    -500.0 | 200000.0 |
|     -6000.0 |     -500.0 | -500000.0 |         3 |    -10.0 |      -5 |      9000.0 |     -500.0 | 7000000.0 |       -2 |
|        30.0 |     9000.0 |         3 |     -30.0 |     0.06 |    0.03 |        50.0 |   200000.0 |        -2 |   1000.0 |


#+begin_src matlab :exports results :results value table replace :tangle no
  data2orgtable(M(1:10, 1:10), {}, {}, ' %.1g ');
#+end_src

#+caption: First 10x10 elements of the Mass matrix
#+RESULTS:
|    0.01 |  -2e-06 |  1e-06 |  6e-09 |  5e-05 | -5e-09 | -0.0005 |  -7e-07 |  6e-07 | -3e-09 |
|  -2e-06 |    0.01 |  8e-07 | -2e-05 | -8e-09 |  2e-09 |  -9e-07 | -0.0002 |  1e-08 | -9e-07 |
|   1e-06 |   8e-07 |  0.009 |  5e-10 |  1e-09 | -1e-09 |  -5e-07 |   3e-08 |  6e-05 |  1e-10 |
|   6e-09 |  -2e-05 |  5e-10 |  3e-07 |  2e-11 | -3e-12 |   3e-09 |   9e-07 | -4e-10 |  3e-09 |
|   5e-05 |  -8e-09 |  1e-09 |  2e-11 |  6e-07 | -4e-11 |  -1e-06 |  -2e-09 |  1e-09 | -8e-12 |
|  -5e-09 |   2e-09 | -1e-09 | -3e-12 | -4e-11 |  1e-07 |  -2e-09 |  -1e-09 | -4e-10 | -5e-12 |
| -0.0005 |  -9e-07 | -5e-07 |  3e-09 | -1e-06 | -2e-09 |    0.01 |   1e-07 | -3e-07 | -2e-08 |
|  -7e-07 | -0.0002 |  3e-08 |  9e-07 | -2e-09 | -1e-09 |   1e-07 |    0.01 | -4e-07 |  2e-05 |
|   6e-07 |   1e-08 |  6e-05 | -4e-10 |  1e-09 | -4e-10 |  -3e-07 |  -4e-07 |  0.009 | -2e-10 |
|  -3e-09 |  -9e-07 |  1e-10 |  3e-09 | -8e-12 | -5e-12 |  -2e-08 |   2e-05 | -2e-10 |  3e-07 |


Then, we extract the coordinates of the interface nodes.
#+begin_src matlab
  [int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('APA300ML_out_nodes_3D.txt');
#+end_src

#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
  data2orgtable([[1:length(int_i)]', int_i, int_xyz], {}, {'Node i', 'Node Number', 'x [m]', 'y [m]', 'z [m]'}, ' %f ');
#+end_src

#+caption: Coordinates of the interface nodes
#+RESULTS:
| Node i | Node Number |   x [m] | y [m] |  z [m] |
|--------+-------------+---------+-------+--------|
|    1.0 |    697783.0 |     0.0 |   0.0 | -0.015 |
|    2.0 |    697784.0 |     0.0 |   0.0 |  0.015 |
|    3.0 |    697785.0 | -0.0325 |   0.0 |    0.0 |
|    4.0 |    697786.0 | -0.0125 |   0.0 |    0.0 |
|    5.0 |    697787.0 | -0.0075 |   0.0 |    0.0 |
|    6.0 |    697788.0 |  0.0125 |   0.0 |    0.0 |
|    7.0 |    697789.0 |  0.0325 |   0.0 |    0.0 |

#+begin_src matlab :exports results :results value table replace :tangle no
  data2orgtable([length(n_i); length(int_i); size(M,1) - 6*length(int_i); size(M,1)], {'Total number of Nodes', 'Number of interface Nodes', 'Number of Modes', 'Size of M and K matrices'}, {}, ' %.0f ');
#+end_src

#+caption: Some extracted parameters of the FEM
#+RESULTS:
| Total number of Nodes     |   7 |
| Number of interface Nodes |   7 |
| Number of Modes           | 120 |
| Size of M and K matrices  | 162 |

Using =K=, =M= and =int_xyz=, we can now use the =Reduced Order Flexible Solid= simscape block.

** Piezoelectric parameters
#+begin_src matlab
  Ga = 1; % [N/V]
  Gs = 1; % [V/m]
#+end_src

#+begin_src matlab
  m = 0.1; % [kg]
#+end_src

** Simscape Model
The flexible element is imported using the =Reduced Order Flexible Solid= simscape block.

Let's say we use two stacks as a force sensor and one stack as an actuator:
- A =Relative Motion Sensor= block is added between the nodes A and C
- An =Internal Force= block is added between the remote points E and B

The interface nodes are shown in Figure [[fig:apa300ml_ansys]].

One mass is fixed at one end of the piezo-electric stack actuator (remove point F), the other end is fixed to the world frame (remote point G).

** Identification of the APA Characteristics
*** Stiffness
#+begin_src matlab :exports none
  m = 0.0001;
#+end_src

The transfer function from vertical external force to the relative vertical displacement is identified.

#+begin_src matlab :exports none
  %% Name of the Simulink File
  mdl = 'APA300ML';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1;
  io(io_i) = linio([mdl, '/z'],  1, 'openoutput'); io_i = io_i + 1;

  G = linearize(mdl, io);
#+end_src

The inverse of its DC gain is the axial stiffness of the APA:
#+begin_src matlab :results replace value
  1e-6/dcgain(G) % [N/um]
#+end_src

#+RESULTS:
: 1.753

The specified stiffness in the datasheet is $k = 1.8\, [N/\mu m]$.

*** Resonance Frequency
The resonance frequency is specified to be between 650Hz and 840Hz.
This is also the case for the FEM model (Figure [[fig:apa300ml_resonance]]).

#+begin_src matlab :exports none
  freqs = logspace(2, 4, 5000);

  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G, freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Amplitude');
  hold off;
#+end_src

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

#+name: fig:apa300ml_resonance
#+caption: First resonance is around 800Hz
#+RESULTS:
[[file:figs/apa300ml_resonance.png]]

*** Amplification factor
The amplification factor is the ratio of the vertical displacement to the stack displacement.

#+begin_src matlab :exports none
  %% Name of the Simulink File
  mdl = 'APA300ML';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
  io(io_i) = linio([mdl, '/z'], 1, 'openoutput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/d'], 1, 'openoutput'); io_i = io_i + 1;

  G = linearize(mdl, io);
#+end_src

The ratio of the two displacement is computed from the FEM model.
#+begin_src matlab :results replace value
  abs(dcgain(G(1,1))./dcgain(G(2,1)))
#+end_src

#+RESULTS:
: 5.0749

This is actually correct and approximately corresponds to the ratio of the piezo height and length:
#+begin_src matlab :results replace value
  75/15
#+end_src

#+RESULTS:
: 5

*** Stroke

Estimation of the actuator stroke:
\[ \Delta H = A n \Delta L \]
with:
- $\Delta H$ Axial Stroke of the APA
- $A$ Amplification factor (5 for the APA300ML)
- $n$ Number of stack used
- $\Delta L$ Stroke of the stack (0.1% of its length)

#+begin_src matlab :results replace value
  1e6 * 5 * 3 * 20e-3 * 0.1e-2
#+end_src

#+RESULTS:
: 300

This is exactly the specified stroke in the data-sheet.

*** TODO Stroke BIS
- [ ] Identified the stroke form the transfer function from V to z

#+begin_src matlab :exports none
  %% Name of the Simulink File
  mdl = 'APA300ML';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/V'], 1, 'openinput');  io_i = io_i + 1;
  io(io_i) = linio([mdl, '/d'], 1, 'openoutput'); io_i = io_i + 1;

  G = linearize(mdl, io);

  1e6*170*abs(dcgain(G))
#+end_src

** Identification of the Dynamics from actuator to replace displacement
We first set the mass to be approximately zero.
#+begin_src matlab :exports none
  m = 0.01;
#+end_src

The dynamics is identified from the applied force to the measured relative displacement.
#+begin_src matlab :exports none
  %% Name of the Simulink File
  mdl = 'APA300ML';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
  io(io_i) = linio([mdl, '/z'], 1, 'openoutput'); io_i = io_i + 1;

  Gh = -linearize(mdl, io);
#+end_src

The same dynamics is identified for a payload mass of 10Kg.
#+begin_src matlab
  m = 10;
#+end_src

#+begin_src matlab :exports none
  %% Name of the Simulink File
  mdl = 'APA300ML';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
  io(io_i) = linio([mdl, '/z'], 1, 'openoutput'); io_i = io_i + 1;

  Ghm = -linearize(mdl, io);
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 4, 5000);

  figure;
  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');

  ax1 = nexttile([2,1]);
  hold on;
  plot(freqs, abs(squeeze(freqresp(Gh, freqs, 'Hz'))), '-');
  plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz'))), '-');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
  hold off;

  ax2 = nexttile;
  hold on;
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gh, freqs, 'Hz')))), '-', ...
       'DisplayName', '$m = 0kg$');
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Ghm, freqs, 'Hz')))), '-', ...
       'DisplayName', '$m = 10kg$');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  yticks(-360:90:360);
  ylim([-360 0]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;
  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
  legend('location', 'southwest');
#+end_src

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

#+name: fig:apa300ml_plant_dynamics
#+caption: Transfer function from forces applied by the stack to the axial displacement of the APA
#+RESULTS:
[[file:figs/apa300ml_plant_dynamics.png]]

The root locus corresponding to Direct Velocity Feedback with a mass of 10kg is shown in Figure [[fig:apa300ml_dvf_root_locus]].
#+begin_src matlab :exports none
  figure;

  gains = logspace(0, 5, 500);

  hold on;
  plot(real(pole(Ghm)),  imag(pole(G)), 'kx');
  plot(real(tzero(Ghm)),  imag(tzero(G)), 'ko');
  for k = 1:length(gains)
      cl_poles = pole(feedback(Ghm, gains(k)*s));
      plot(real(cl_poles), imag(cl_poles), 'k.');
  end
  hold off;
  axis square;
  xlim([-500, 10]); ylim([0, 510]);

  xlabel('Real Part'); ylabel('Imaginary Part');
#+end_src

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

#+name: fig:apa300ml_dvf_root_locus
#+caption: Root Locus for Direct Velocity Feedback
#+RESULTS:
[[file:figs/apa300ml_dvf_root_locus.png]]

** Identification of the Dynamics from actuator to force sensor
Let's use 2 stacks as a force sensor and 1 stack as force actuator.

The transfer function from actuator voltage to sensor voltage is identified and shown in Figure [[fig:apa300ml_iff_plant]].
#+begin_src matlab :exports none
  m = 10;
#+end_src

#+begin_src matlab :exports none
  %% Name of the Simulink File
  mdl = 'APA300ML';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Va'], 1, 'openinput');  io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Vs'], 1, 'openoutput'); io_i = io_i + 1;

  Giff = -linearize(mdl, io);
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 4, 5000);

  figure;
  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');

  ax1 = nexttile([2,1]);
  hold on;
  plot(freqs, abs(squeeze(freqresp(Giff, freqs, 'Hz'))), '-');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
  hold off;

  ax2 = nexttile;
  hold on;
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Giff, freqs, 'Hz')))), '-');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  yticks(-360:90:360);
  ylim([-180 180]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;
  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

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

#+name: fig:apa300ml_iff_plant
#+caption: Transfer function from actuator to force sensor
#+RESULTS:
[[file:figs/apa300ml_iff_plant.png]]

For root locus corresponding to IFF is shown in Figure [[fig:apa300ml_iff_root_locus]].
#+begin_src matlab :exports none
  figure;

  gains = logspace(0, 5, 500);

  hold on;
  plot(real(pole(Giff)),  imag(pole(Giff)), 'kx');
  plot(real(tzero(Giff)),  imag(tzero(Giff)), 'ko');
  for k = 1:length(gains)
      cl_poles = pole(feedback(Giff, gains(k)/s));
      plot(real(cl_poles), imag(cl_poles), 'k.');
  end
  hold off;
  axis square;
  xlim([-500, 10]); ylim([0, 510]);

  xlabel('Real Part'); ylabel('Imaginary Part');
#+end_src

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

#+name: fig:apa300ml_iff_root_locus
#+caption: Root Locus for IFF
#+RESULTS:
[[file:figs/apa300ml_iff_root_locus.png]]

** Identification for a simpler model
The goal in this section is to identify the parameters of a simple APA model from the FEM.
This can be useful is a lower order model is to be used for simulations.

The presented model is based on cite:souleille18_concep_activ_mount_space_applic.

The model represents the Amplified Piezo Actuator (APA) from Cedrat-Technologies (Figure [[fig:souleille18_model_piezo]]).
The parameters are shown in the table below.

#+name: fig:souleille18_model_piezo
#+caption: Picture of an APA100M from Cedrat Technologies. Simplified model of a one DoF payload mounted on such isolator
[[file:./figs/souleille18_model_piezo.png]]

#+caption:Parameters used for the model of the APA 100M
|       | Meaning                                                        |
|-------+----------------------------------------------------------------|
| $k_e$ | Stiffness used to adjust the pole of the isolator              |
| $k_1$ | Stiffness of the metallic suspension when the stack is removed |
| $k_a$ | Stiffness of the actuator                                      |
| $c_1$ | Added viscous damping                                          |

The goal is to determine $k_e$, $k_a$ and $k_1$ so that the simplified model fits the FEM model.

\[ \alpha = \frac{x_1}{f}(\omega=0) = \frac{\frac{k_e}{k_e + k_a}}{k_1 + \frac{k_e k_a}{k_e + k_a}} \]
\[ \beta = \frac{x_1}{F}(\omega=0) = \frac{1}{k_1 + \frac{k_e k_a}{k_e + k_a}} \]

If we can fix $k_a$, we can determine $k_e$ and $k_1$ with:
\[ k_e = \frac{k_a}{\frac{\beta}{\alpha} - 1} \]
\[ k_1 = \frac{1}{\beta} - \frac{k_e k_a}{k_e + k_a} \]

#+begin_src matlab :exports none
  m = 10;
#+end_src

#+begin_src matlab :exports none
  %% Name of the Simulink File
  mdl = 'APA300ML';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1; % External Vertical Force [N]
  io(io_i) = linio([mdl, '/w'],  1, 'openinput');  io_i = io_i + 1; % Base Motion [m]
  io(io_i) = linio([mdl, '/Fa'], 1, 'openinput');  io_i = io_i + 1; % Actuator Force [N]
  io(io_i) = linio([mdl, '/z'],  1, 'openoutput'); io_i = io_i + 1; % Vertical Displacement [m]
  io(io_i) = linio([mdl, '/Vs'], 1, 'openoutput'); io_i = io_i + 1; % Force Sensor [V]
  io(io_i) = linio([mdl, '/d'],  1, 'openoutput'); io_i = io_i + 1; % Stack Displacement [m]

  G = linearize(mdl, io);

  G.InputName = {'Fd', 'w', 'Fa'};
  G.OutputName = {'y', 'Fs', 'd'};
#+end_src

From the identified dynamics, compute $\alpha$ and $\beta$
#+begin_src matlab
  alpha = abs(dcgain(G('y', 'Fa')));
  beta  = abs(dcgain(G('y', 'Fd')));
#+end_src

$k_a$ is estimated using the following formula:
#+begin_src matlab
  ka = 0.8/abs(dcgain(G('y', 'Fa')));
#+end_src
The factor can be adjusted to better match the curves.

Then $k_e$ and $k_1$ are computed.
#+begin_src matlab
  ke = ka/(beta/alpha - 1);
  k1 = 1/beta - ke*ka/(ke + ka);
#+end_src

#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
  data2orgtable(1e-6*[ka; ke; k1], {'ka', 'ke', 'k1'}, {'Value [N/um]'}, ' %.1f ');
#+end_src

#+RESULTS:
|    | Value [N/um] |
|----+--------------|
| ka |         40.5 |
| ke |          1.5 |
| k1 |          0.4 |

The damping in the system is adjusted to match the FEM model if necessary.
#+begin_src matlab
  c1 = 1e2;
#+end_src

The analytical model of the simpler system is defined below:
#+begin_src matlab
  Ga = 1/(m*s^2 + k1 + c1*s + ke*ka/(ke + ka)) * ...
       [ 1 ,              k1 + c1*s + ke*ka/(ke + ka)  , ke/(ke + ka) ;
        -ke*ka/(ke + ka), ke*ka/(ke + ka)*m*s^2 ,       -ke/(ke + ka)*(m*s^2 + c1*s + k1)];

  Ga.InputName = {'Fd', 'w', 'Fa'};
  Ga.OutputName = {'y', 'Fs'};
#+end_src

And the DC gain is adjusted for the force sensor:
#+begin_src matlab
  F_gain = dcgain(G('Fs', 'Fd'))/dcgain(Ga('Fs', 'Fd'));
#+end_src

The dynamics of the FEM model and the simpler model are compared in Figure [[fig:apa300ml_comp_simpler_model]].

#+begin_src matlab :exports none
  freqs = logspace(0, 5, 1000);

  figure;
  tiledlayout(2, 3, 'TileSpacing', 'None', 'Padding', 'None');

  ax1 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G( 'y', 'w'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Ga('y', 'w'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  ylabel('$x_1/w$ [m/m]');
  ylim([1e-6, 1e2]);

  ax2 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fa'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Ga('y', 'Fa'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  ylabel('$x_1/f$ [m/N]');
  ylim([1e-14, 1e-6]);

  ax3 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fd'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Ga('y', 'Fd'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  ylabel('$x_1/F$ [m/N]');
  ylim([1e-14, 1e-4]);

  ax4 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'w'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(F_gain*Ga('Fs', 'w'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]');
  ylabel('$F_s/w$ [m/m]');
  ylim([1e2, 1e8]);

  ax5 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fa'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(F_gain*Ga('Fs', 'Fa'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]');
  ylabel('$F_s/f$ [m/N]');
  ylim([1e-4, 1e1]);

  ax6 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fd'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(F_gain*Ga('Fs', 'Fd'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]');
  ylabel('$F_s/F$ [m/N]');
  ylim([1e-7, 1e2]);

  linkaxes([ax1,ax2,ax3,ax4,ax5,ax6],'x');
#+end_src

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

#+name: fig:apa300ml_comp_simpler_model
#+caption: Comparison of the Dynamics between the FEM model and the simplified one
#+RESULTS:
[[file:figs/apa300ml_comp_simpler_model.png]]

The simplified model has also been implemented in Simscape.

The dynamics of the Simscape simplified model is identified and compared with the FEM one in Figure [[fig:apa300ml_comp_simpler_simscape]].
#+begin_src matlab :exports none
  %% Name of the Simulink File
  mdl = 'APA300ML_simplified';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1; % External Vertical Force [N]
  io(io_i) = linio([mdl, '/w'],  1, 'openinput');  io_i = io_i + 1; % Base Motion [m]
  io(io_i) = linio([mdl, '/Fa'], 1, 'openinput');  io_i = io_i + 1; % Actuator Force [N]
  io(io_i) = linio([mdl, '/y'],  1, 'openoutput'); io_i = io_i + 1; % Vertical Displacement [m]
  io(io_i) = linio([mdl, '/Fs'], 1, 'openoutput'); io_i = io_i + 1; % Force Sensor [V]

  Gs = linearize(mdl, io);

  Gs.InputName = {'Fd', 'w', 'Fa'};
  Gs.OutputName = {'y', 'Fs'};
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 5, 1000);

  figure;
  tiledlayout(2, 3, 'TileSpacing', 'None', 'Padding', 'None');

  ax1 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G( 'y', 'w'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Gs('y', 'w'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  ylabel('$x_1/w$ [m/m]');
  ylim([1e-6, 1e2]);

  ax2 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fa'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Gs('y', 'Fa'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  ylabel('$x_1/f$ [m/N]');
  ylim([1e-14, 1e-6]);

  ax3 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fd'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Gs('y', 'Fd'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]);
  ylabel('$x_1/F$ [m/N]');
  ylim([1e-14, 1e-4]);

  ax4 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'w'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(F_gain*Gs('Fs', 'w'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]');
  ylabel('$F_s/w$ [m/m]');
  ylim([1e2, 1e8]);

  ax5 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fa'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(F_gain*Gs('Fs', 'Fa'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]');
  ylabel('$F_s/f$ [m/N]');
  ylim([1e-4, 1e1]);

  ax6 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fd'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(F_gain*Gs('Fs', 'Fd'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]');
  ylabel('$F_s/F$ [m/N]');
  ylim([1e-7, 1e2]);

  linkaxes([ax1,ax2,ax3,ax4,ax5,ax6],'x');
#+end_src

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

#+name: fig:apa300ml_comp_simpler_simscape
#+caption: Comparison of the Dynamics between the FEM model and the simplified simscape model
#+RESULTS:
[[file:figs/apa300ml_comp_simpler_simscape.png]]

** Integral Force Feedback
In this section, Integral Force Feedback control architecture is applied on the APA300ML.

First, the plant (dynamics from voltage actuator to voltage sensor is identified).
#+begin_src matlab :exports none
  Kiff = tf(0);
#+end_src

The payload mass is set to 10kg.
#+begin_src matlab
  m = 10;
#+end_src

#+begin_src matlab :exports none
  %% Name of the Simulink File
  mdl = 'APA300ML_IFF';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/w'],         1, 'openinput');  io_i = io_i + 1;
  io(io_i) = linio([mdl, '/F'],         1, 'openinput');  io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Fd'],        1, 'openinput');  io_i = io_i + 1;
  io(io_i) = linio([mdl, '/z'],         1, 'openoutput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/APA300ML'],  1, 'openoutput'); io_i = io_i + 1;

  G_ol = linearize(mdl, io);
  G_ol.InputName  = {'w', 'f', 'F'};
  G_ol.OutputName  = {'x1', 'Fs'};

  G = G_ol({'Fs'}, {'f'});
#+end_src

The obtained dynamics is shown in Figure [[fig:piezo_amplified_iff_plant]].

#+begin_src matlab :exports none
  freqs = logspace(1, 5, 1000);

  figure;
  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');

  ax1 = nexttile([2,1]);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G, freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
  hold off;

  ax2 = nexttile;
  hold on;
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G, freqs, 'Hz')))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  yticks(-360:90:360);
  ylim([-390 30]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;
  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

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

#+name: fig:piezo_amplified_iff_plant
#+caption: IFF Plant
#+RESULTS:
[[file:figs/piezo_amplified_iff_plant.png]]

The controller is defined below and the loop gain is shown in Figure [[fig:piezo_amplified_iff_loop_gain]].
#+begin_src matlab
  Kiff = -1e3/s;
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(1, 5, 1000);

  figure;
  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');

  ax1 = nexttile([2,1]);
  hold on;
  plot(freqs, abs(squeeze(freqresp(G*Kiff, freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
  hold off;

  ax2 = nexttile;
  hold on;
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G*Kiff, freqs, 'Hz')))));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
  yticks(-360:90:360);
  ylim([-180 180]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;
  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

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

#+name: fig:piezo_amplified_iff_loop_gain
#+caption: IFF Loop Gain
#+RESULTS:
[[file:figs/piezo_amplified_iff_loop_gain.png]]

Now the closed-loop system is identified again and compare with the open loop system in Figure [[fig:piezo_amplified_iff_comp]].

It is the expected behavior as shown in the Figure [[fig:souleille18_results]] (from cite:souleille18_concep_activ_mount_space_applic).

#+begin_src matlab :exports none
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/w'],         1, 'openinput');  io_i = io_i + 1;
  io(io_i) = linio([mdl, '/F'],         1, 'openinput');  io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Fd'],        1, 'openinput');  io_i = io_i + 1;
  io(io_i) = linio([mdl, '/z'],         1, 'openoutput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/APA300ML'],  1, 'output');     io_i = io_i + 1;

  Giff = linearize(mdl, io);
  Giff.InputName  = {'w', 'f', 'F'};
  Giff.OutputName  = {'x1', 'Fs'};
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 3, 1000);

  figure;
  tiledlayout(2, 3, 'TileSpacing', 'None', 'Padding', 'None');

  ax1 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_ol('x1', 'w'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Giff('x1', 'w'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]); ylabel('$x_1/w$ [m/m]')

  ax2 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_ol('x1', 'f'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Giff('x1', 'f'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]); ylabel('$x_1/f$ [m/N]');

  ax3 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_ol('x1', 'F'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Giff('x1', 'F'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  set(gca, 'XTickLabel',[]); ylabel('$x_1/F$ [m/N]');

  ax4 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_ol('Fs', 'w'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Giff('Fs', 'w'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('$F_s/w$ [N/m]');

  ax5 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_ol('Fs', 'f'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Giff('Fs', 'f'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('$F_s/f$ [N/N]');

  ax6 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_ol('Fs', 'F'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Giff('Fs', 'F'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('$F_s/F$ [N/N]');

  linkaxes([ax1,ax2,ax3,ax4,ax5,ax6],'x');
#+end_src

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

#+name: fig:piezo_amplified_iff_comp
#+caption: OL and CL transfer functions
#+RESULTS:
[[file:figs/piezo_amplified_iff_comp.png]]

#+name: fig:souleille18_results
#+caption: Results obtained in cite:souleille18_concep_activ_mount_space_applic
[[file:figs/souleille18_results.png]]

* 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.jpg]]

** 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.jpg]]

The data is loaded.
#+begin_src matlab
bending_X = load('apa300ml_bending_X_top.mat');
#+end_src

The config for =tfestimate= is performed:
#+begin_src matlab
Ts = bending_X.Track1_X_Resolution; % Sampling frequency [Hz]
win = hann(ceil(1/Ts));
#+end_src

The transfer function from the input force to the output "rotation" (difference between the two measured distances).
#+begin_src matlab
[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.jpg]]

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.jpg]]

The data is loaded, and the transfer function computed.
#+begin_src matlab
torsion = load('apa300ml_torsion_left.mat');
[G_torsion, ~]  = tfestimate(torsion.Track1, torsion.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')
text(800, 6e-4,{'800Hz'}, '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]]

In order to verify that, the APA is excited on the top part such that the torsion mode should not be excited.
#+begin_src matlab
torsion = load('apa300ml_torsion_top.mat');
[G_torsion_top, ~]  = tfestimate(torsion.Track1, torsion.Track2, win, [], [], 1/Ts);
#+end_src

The two FRF are compared in Figure [[fig:apa300ml_meas_freq_torsion_z_comp]].
It is clear that the first two modes does not correspond to the torsional mode.
Maybe the resonance at 800Hz, or even higher resonances. It is difficult to conclude here.
#+begin_src matlab :exports none
figure;
hold on;
plot(f, abs(G_torsion), 'k-', 'DisplayName', 'Left excitation');
plot(f, abs(G_torsion_top), '-', 'DisplayName', 'Top excitation');
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')
text(800, 2e-3,{'800Hz'}, 'VerticalAlignment', 'bottom','HorizontalAlignment','center')
legend('location', 'northwest');
#+end_src

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

#+name: fig:apa300ml_meas_freq_torsion_z_comp
#+caption: Obtained FRF for the Z-torsion
#+RESULTS:
[[file:figs/apa300ml_meas_freq_torsion_z_comp.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

When two flexible joints are fixed at each ends of the APA, the APA is mostly in a free/free condition in terms of bending/torsion (the bending/torsional stiffness of the joints being very small).

In the current tests, the APA are in a fixed/free condition.
Therefore, it is quite obvious that we measured higher resonance frequencies than what is foreseen for the struts.
It is however quite interesting that there is a factor $\approx \sqrt{2}$ between the two (increased of the stiffness by a factor 2?).

#+name: tab:apa300ml_measured_modes_freq
#+caption: Measured frequency of the modes
#+attr_latex: :environment tabularx :width 0.6\linewidth :align ccc
#+attr_latex: :center t :booktabs t :float t
| Mode      | Strut Mode | Measured Frequency |
|-----------+------------+--------------------|
| X-Bending | 189Hz      | 280Hz              |
| Y-Bending | 285Hz      | 410Hz              |
| Z-Torsion | 400Hz      | ?                  |

* Function
** =generateSweepExc=: Generate sweep sinus excitation
:PROPERTIES:
:header-args:matlab+: :tangle ./matlab/src/generateSweepExc.m
:header-args:matlab+: :comments none :mkdirp yes :eval no
:END:
<<sec:generateSweepExc>>

*** Function description
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
function [U_exc] = generateSweepExc(args)
% generateSweepExc - Generate a Sweep Sine excitation signal
%
% Syntax: [U_exc] = generateSweepExc(args)
%
% Inputs:
%    - args - Optinal arguments:
%        - Ts              - Sampling Time                                              - [s]
%        - f_start         - Start frequency of the sweep                               - [Hz]
%        - f_end           - End frequency of the sweep                                 - [Hz]
%        - V_mean          - Mean value of the excitation voltage                       - [V]
%        - V_exc           - Excitation Amplitude for the Sweep, could be numeric or TF - [V]
%        - t_start         - Time at which the sweep begins                             - [s]
%        - exc_duration    - Duration of the sweep                                      - [s]
%        - sweep_type      - 'logarithmic' or 'linear'                                  - [-]
%        - smooth_ends     - 'true' or 'false': smooth transition between 0 and V_mean  - [-]
#+end_src

*** Optional Parameters
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
arguments
    args.Ts              (1,1) double  {mustBeNumeric, mustBePositive} = 1e-4
    args.f_start         (1,1) double  {mustBeNumeric, mustBePositive} = 1
    args.f_end           (1,1) double  {mustBeNumeric, mustBePositive} = 1e3
    args.V_mean          (1,1) double  {mustBeNumeric} = 0
    args.V_exc                                         = 1
    args.t_start         (1,1) double  {mustBeNumeric, mustBeNonnegative} = 5
    args.exc_duration    (1,1) double  {mustBeNumeric, mustBePositive} = 10
    args.sweep_type            char    {mustBeMember(args.sweep_type,{'log', 'lin'})} = 'lin'
    args.smooth_ends           logical {mustBeNumericOrLogical} = true
end
#+end_src

*** Sweep Sine part
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
t_sweep = 0:args.Ts:args.exc_duration;

if strcmp(args.sweep_type, 'log')
    V_exc = sin(2*pi*args.f_start * args.exc_duration/log(args.f_end/args.f_start) * (exp(log(args.f_end/args.f_start)*t_sweep/args.exc_duration) - 1));
elseif strcmp(args.sweep_type, 'lin')
    V_exc = sin(2*pi*(args.f_start + (args.f_end - args.f_start)/2/args.exc_duration*t_sweep).*t_sweep);
else
    error('sweep_type should either be equal to "log" or to "lin"');
end
#+end_src

#+begin_src matlab
if isnumeric(args.V_exc)
    V_sweep = args.V_mean + args.V_exc*V_exc;
elseif isct(args.V_exc)
    if strcmp(args.sweep_type, 'log')
        V_sweep = args.V_mean + abs(squeeze(freqresp(args.V_exc, args.f_start*(args.f_end/args.f_start).^(t_sweep/args.exc_duration), 'Hz')))'.*V_exc;
    elseif strcmp(args.sweep_type, 'lin')
        V_sweep = args.V_mean + abs(squeeze(freqresp(args.V_exc, args.f_start+(args.f_end-args.f_start)/args.exc_duration*t_sweep, 'Hz')))'.*V_exc;
    end
end
#+end_src

*** Smooth Ends
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
if args.t_start > 0
    t_smooth_start = args.Ts:args.Ts:args.t_start;

    V_smooth_start = zeros(size(t_smooth_start));
    V_smooth_end   = zeros(size(t_smooth_start));

    if args.smooth_ends
        Vd_max = args.V_mean/(0.7*args.t_start);

        V_d = zeros(size(t_smooth_start));
        V_d(t_smooth_start < 0.2*args.t_start) = t_smooth_start(t_smooth_start < 0.2*args.t_start)*Vd_max/(0.2*args.t_start);
        V_d(t_smooth_start > 0.2*args.t_start & t_smooth_start < 0.7*args.t_start) = Vd_max;
        V_d(t_smooth_start > 0.7*args.t_start & t_smooth_start < 0.9*args.t_start) = Vd_max - (t_smooth_start(t_smooth_start > 0.7*args.t_start & t_smooth_start < 0.9*args.t_start) - 0.7*args.t_start)*Vd_max/(0.2*args.t_start);

        V_smooth_start = cumtrapz(V_d)*args.Ts;

        V_smooth_end = args.V_mean - V_smooth_start;
    end
else
    V_smooth_start = [];
    V_smooth_end = [];
end
#+end_src

*** Combine Excitation signals
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
V_exc = [V_smooth_start, V_sweep, V_smooth_end];
t_exc = args.Ts*[0:1:length(V_exc)-1];
#+end_src

#+begin_src matlab
U_exc = [t_exc; V_exc];
#+end_src

** =generateShapedNoise=: Generate Shaped Noise excitation
:PROPERTIES:
:header-args:matlab+: :tangle ./matlab/src/generateShapedNoise.m
:header-args:matlab+: :comments none :mkdirp yes :eval no
:END:
<<sec:generateShapedNoise>>

*** Function description
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
function [U_exc] = generateShapedNoise(args)
% generateShapedNoise - Generate a Shaped Noise excitation signal
%
% Syntax: [U_exc] = generateShapedNoise(args)
%
% Inputs:
%    - args - Optinal arguments:
%        - Ts              - Sampling Time                                              - [s]
%        - V_mean          - Mean value of the excitation voltage                       - [V]
%        - V_exc           - Excitation Amplitude, could be numeric or TF               - [V rms]
%        - t_start         - Time at which the noise begins                             - [s]
%        - exc_duration    - Duration of the noise                                      - [s]
%        - smooth_ends     - 'true' or 'false': smooth transition between 0 and V_mean  - [-]
#+end_src

*** Optional Parameters
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
arguments
    args.Ts              (1,1) double  {mustBeNumeric, mustBePositive} = 1e-4
    args.V_mean          (1,1) double  {mustBeNumeric} = 0
    args.V_exc                                         = 1
    args.t_start         (1,1) double  {mustBeNumeric, mustBePositive} = 5
    args.exc_duration    (1,1) double  {mustBeNumeric, mustBePositive} = 10
    args.smooth_ends           logical {mustBeNumericOrLogical} = true
end
#+end_src

*** Shaped Noise
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
t_noise = 0:args.Ts:args.exc_duration;

#+end_src

#+begin_src matlab
if isnumeric(args.V_exc)
    V_noise = args.V_mean + args.V_exc*sqrt(1/args.Ts/2)*randn(length(t_noise), 1)';
elseif isct(args.V_exc)
    V_noise = args.V_mean + lsim(args.V_exc, sqrt(1/args.Ts/2)*randn(length(t_noise), 1), t_noise)';
end
#+end_src

*** Smooth Ends
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
t_smooth_start = args.Ts:args.Ts:args.t_start;

V_smooth_start = zeros(size(t_smooth_start));
V_smooth_end   = zeros(size(t_smooth_start));

if args.smooth_ends
    Vd_max = args.V_mean/(0.7*args.t_start);

    V_d = zeros(size(t_smooth_start));
    V_d(t_smooth_start < 0.2*args.t_start) = t_smooth_start(t_smooth_start < 0.2*args.t_start)*Vd_max/(0.2*args.t_start);
    V_d(t_smooth_start > 0.2*args.t_start & t_smooth_start < 0.7*args.t_start) = Vd_max;
    V_d(t_smooth_start > 0.7*args.t_start & t_smooth_start < 0.9*args.t_start) = Vd_max - (t_smooth_start(t_smooth_start > 0.7*args.t_start & t_smooth_start < 0.9*args.t_start) - 0.7*args.t_start)*Vd_max/(0.2*args.t_start);

    V_smooth_start = cumtrapz(V_d)*args.Ts;

    V_smooth_end = args.V_mean - V_smooth_start;
end
#+end_src

*** Combine Excitation signals
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
V_exc = [V_smooth_start, V_noise, V_smooth_end];
t_exc = args.Ts*[0:1:length(V_exc)-1];
#+end_src

#+begin_src matlab
U_exc = [t_exc; V_exc];
#+end_src

** =generateSinIncreasingAmpl=: Generate Sinus with increasing amplitude
:PROPERTIES:
:header-args:matlab+: :tangle ./matlab/src/generateSinIncreasingAmpl.m
:header-args:matlab+: :comments none :mkdirp yes :eval no
:END:
<<sec:generateSinIncreasingAmpl>>

*** Function description
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
function [U_exc] = generateSinIncreasingAmpl(args)
% generateSinIncreasingAmpl - Generate Sinus with increasing amplitude
%
% Syntax: [U_exc] = generateSinIncreasingAmpl(args)
%
% Inputs:
%    - args - Optinal arguments:
%        - Ts              - Sampling Time                                              - [s]
%        - V_mean          - Mean value of the excitation voltage                       - [V]
%        - sin_ampls       - Excitation Amplitudes                                      - [V]
%        - sin_freq        - Excitation Frequency                                       - [Hz]
%        - sin_num         - Number of period for each amplitude                        - [-]
%        - t_start         - Time at which the excitation begins                        - [s]
%        - smooth_ends     - 'true' or 'false': smooth transition between 0 and V_mean  - [-]
#+end_src

*** Optional Parameters
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
arguments
    args.Ts              (1,1) double  {mustBeNumeric, mustBePositive} = 1e-4
    args.V_mean          (1,1) double  {mustBeNumeric} = 0
    args.sin_ampls             double  {mustBeNumeric, mustBePositive} = [0.1, 0.2, 0.3]
    args.sin_period      (1,1) double  {mustBeNumeric, mustBePositive} = 1
    args.sin_num         (1,1) double  {mustBeNumeric, mustBePositive, mustBeInteger} = 3
    args.t_start         (1,1) double  {mustBeNumeric, mustBePositive} = 5
    args.smooth_ends           logical {mustBeNumericOrLogical} = true
end
#+end_src

*** Sinus excitation
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
t_noise = 0:args.Ts:args.sin_period*args.sin_num;
sin_exc = [];
#+end_src

#+begin_src matlab
for sin_ampl = args.sin_ampls
    sin_exc = [sin_exc, args.V_mean + sin_ampl*sin(2*pi/args.sin_period*t_noise)];
end
#+end_src

*** Smooth Ends
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
t_smooth_start = args.Ts:args.Ts:args.t_start;

V_smooth_start = zeros(size(t_smooth_start));
V_smooth_end   = zeros(size(t_smooth_start));

if args.smooth_ends
    Vd_max = args.V_mean/(0.7*args.t_start);

    V_d = zeros(size(t_smooth_start));
    V_d(t_smooth_start < 0.2*args.t_start) = t_smooth_start(t_smooth_start < 0.2*args.t_start)*Vd_max/(0.2*args.t_start);
    V_d(t_smooth_start > 0.2*args.t_start & t_smooth_start < 0.7*args.t_start) = Vd_max;
    V_d(t_smooth_start > 0.7*args.t_start & t_smooth_start < 0.9*args.t_start) = Vd_max - (t_smooth_start(t_smooth_start > 0.7*args.t_start & t_smooth_start < 0.9*args.t_start) - 0.7*args.t_start)*Vd_max/(0.2*args.t_start);

    V_smooth_start = cumtrapz(V_d)*args.Ts;

    V_smooth_end = args.V_mean - V_smooth_start;
end
#+end_src

*** Combine Excitation signals
:PROPERTIES:
:UNNUMBERED: t
:END:

#+begin_src matlab
V_exc = [V_smooth_start, sin_exc, V_smooth_end];
t_exc = args.Ts*[0:1:length(V_exc)-1];
#+end_src

#+begin_src matlab
U_exc = [t_exc; V_exc];
#+end_src

* Bibliography                                                        :ignore:
#+latex: \printbibliography