#+TITLE: Amplified Piezoelectric Stack Actuator
#+SETUPFILE: ./setup/org-setup-file.org

* Introduction                                                        :ignore:
The presented model is based on cite:souleille18_concep_activ_mount_space_applic.

The model represents the amplified piezo APA100M 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
|       | Value             | Meaning                                                        |
|-------+-------------------+----------------------------------------------------------------|
| $m$   | $1\,[kg]$         | Payload mass                                                   |
| $k_e$ | $4.8\,[N/\mu m]$  | Stiffness used to adjust the pole of the isolator              |
| $k_1$ | $0.96\,[N/\mu m]$ | Stiffness of the metallic suspension when the stack is removed |
| $k_a$ | $65\,[N/\mu m]$   | Stiffness of the actuator                                      |
| $c_1$ | $10\,[N/(m/s)]$   | Added viscous damping                                          |

* Simplified Model
** Matlab Init                                              :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
  <<matlab-dir>>
#+end_src

#+begin_src matlab :exports none :results silent :noweb yes
  <<matlab-init>>
#+end_src

#+BEGIN_SRC matlab
  simulinkproject('../');
#+END_SRC

#+begin_src matlab
  open 'amplified_piezo_model.slx'
#+end_src

** Parameters
#+begin_src matlab
  m = 1; % [kg]

  ke = 4.8e6; % [N/m]
  ce = 5; % [N/(m/s)]
  me = 0.001; % [kg]

  k1 = 0.96e6; % [N/m]
  c1 = 10; % [N/(m/s)]

  ka = 65e6; % [N/m]
  ca = 5; % [N/(m/s)]
  ma = 0.001; % [kg]

  h = 0.2; % [m]
#+end_src

IFF Controller:
#+begin_src matlab
  Kiff = -8000/s;
#+end_src

** Identification
Identification in open-loop.
#+begin_src matlab
  %% Name of the Simulink File
  mdl = 'amplified_piezo_model';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/w'],  1, 'openinput');  io_i = io_i + 1; % Base Motion
  io(io_i) = linio([mdl, '/f'],  1, 'openinput');  io_i = io_i + 1; % Actuator Inputs
  io(io_i) = linio([mdl, '/F'],  1, 'openinput');  io_i = io_i + 1; % External Force

  io(io_i) = linio([mdl, '/Fs'], 3, 'openoutput'); io_i = io_i + 1; % Force Sensors
  io(io_i) = linio([mdl, '/x1'], 1, 'openoutput'); io_i = io_i + 1; % Mass displacement

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

Identification in closed-loop.
#+begin_src matlab
  %% Name of the Simulink File
  mdl = 'amplified_piezo_model';

  %% Input/Output definition
  clear io; io_i = 1;
  io(io_i) = linio([mdl, '/w'],  1, 'input');  io_i = io_i + 1; % Base Motion
  io(io_i) = linio([mdl, '/f'],  1, 'input');  io_i = io_i + 1; % Actuator Inputs
  io(io_i) = linio([mdl, '/F'],  1, 'input');  io_i = io_i + 1; % External Force

  io(io_i) = linio([mdl, '/Fs'], 3, 'output'); io_i = io_i + 1; % Force Sensors
  io(io_i) = linio([mdl, '/x1'], 1, 'output'); io_i = io_i + 1; % Mass displacement

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

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

  figure;

  ax1 = subplot(2, 3, 1);
  title('$\displaystyle \frac{x_1}{w}$')
  hold on;
  plot(freqs, abs(squeeze(freqresp(G('x1', 'w'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Giff('x1', 'w'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]');xlabel('Frequency [Hz]');

  ax2 = subplot(2, 3, 2);
  title('$\displaystyle \frac{x_1}{f}$')
  hold on;
  plot(freqs, abs(squeeze(freqresp(G('x1', 'f'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Giff('x1', 'f'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');

  ax3 = subplot(2, 3, 3);
  title('$\displaystyle \frac{x_1}{F}$')
  hold on;
  plot(freqs, abs(squeeze(freqresp(G('x1', 'F'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Giff('x1', 'F'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');

  ax4 = subplot(2, 3, 4);
  title('$\displaystyle \frac{F_s}{w}$')
  hold on;
  plot(freqs, abs(squeeze(freqresp(G('Fs', 'w'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Giff('Fs', 'w'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]');xlabel('Frequency [Hz]');

  ax5 = subplot(2, 3, 5);
  title('$\displaystyle \frac{F_s}{f}$')
  hold on;
  plot(freqs, abs(squeeze(freqresp(G('Fs', 'f'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Giff('Fs', 'f'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');

  ax6 = subplot(2, 3, 6);
  title('$\displaystyle \frac{F_s}{F}$')
  hold on;
  plot(freqs, abs(squeeze(freqresp(G('Fs', 'F'), freqs, 'Hz'))));
  plot(freqs, abs(squeeze(freqresp(Giff('Fs', 'F'), freqs, 'Hz'))));
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');

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

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

#+name: fig:amplified_piezo_tf_ol_and_cl
#+caption: Matrix of transfer functions from input to output in open loop (blue) and closed loop (red)
#+RESULTS:
[[file:figs/amplified_piezo_tf_ol_and_cl.png]]

** Root Locus
#+begin_src matlab :exports none :post
  figure;

  gains = logspace(1, 6, 500);

  hold on;
  plot(real(pole(G('Fs', 'f'))),  imag(pole(G('Fs', 'f'))), 'kx');
  plot(real(tzero(G('Fs', 'f'))),  imag(tzero(G('Fs', 'f'))), 'ko');
  for k = 1:length(gains)
      cl_poles = pole(feedback(G('Fs', 'f'), -gains(k)/s));
      plot(real(cl_poles), imag(cl_poles), 'k.');
  end
  hold off;
  axis square;
  xlim([-2500, 100]); ylim([0, 2600]);

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

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

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

* Rotating X-Y platform
** 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
  simulinkproject('../');
#+END_SRC

#+begin_src matlab
  open 'amplified_piezo_xy_rotating_stage.slx'
#+end_src

** Parameters
#+begin_src matlab
  m = 1; % [kg]

  ke = 4.8e6; % [N/m]
  ce = 5; % [N/(m/s)]
  me = 0.001; % [kg]

  k1 = 0.96e6; % [N/m]
  c1 = 10; % [N/(m/s)]

  ka = 65e6; % [N/m]
  ca = 5; % [N/(m/s)]
  ma = 0.001; % [kg]

  h = 0.2; % [m]
#+end_src

#+begin_src matlab
  Kiff = tf(0);
#+end_src

** Identification
Rotating speed in rad/s:
#+begin_src matlab
  Ws = 2*pi*[0, 1, 10, 100];
#+end_src

#+begin_src matlab
  Gs = {zeros(length(Ws), 1)};
#+end_src

Identification in open-loop.
#+begin_src matlab
  %% Name of the Simulink File
  mdl = 'amplified_piezo_xy_rotating_stage';

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

  io(io_i) = linio([mdl, '/Fs'], 1, 'openoutput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Fs'], 2, 'openoutput'); io_i = io_i + 1;

  for i = 1:length(Ws)
      ws = Ws(i);
      G = linearize(mdl, io, 0);
      G.InputName  = {'fx', 'fy'};
      G.OutputName = {'Fsx', 'Fsy'};
      Gs(i) = {G};
  end
#+end_src

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

  figure;

  ax1 = subplot(2, 2, 1);
  title('$\displaystyle \frac{F_{s,x}}{f_x}$')
  hold on;
  for i = 1:length(Ws)
    plot(freqs, abs(squeeze(freqresp(Gs{i}('Fsx', 'fx'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]');xlabel('Frequency [Hz]');

  ax2 = subplot(2, 2, 2);
  title('$\displaystyle \frac{F_{s,y}}{f_x}$')
  hold on;
  for i = 1:length(Ws)
    plot(freqs, abs(squeeze(freqresp(Gs{i}('Fsy', 'fx'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');

  ax3 = subplot(2, 2, 3);
  title('$\displaystyle \frac{F_{s,x}}{f_y}$')
  hold on;
  for i = 1:length(Ws)
    plot(freqs, abs(squeeze(freqresp(Gs{i}('Fsx', 'fy'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');

  ax4 = subplot(2, 2, 4);
  title('$\displaystyle \frac{F_{s,y}}{f_y}$')
  hold on;
  for i = 1:length(Ws)
    plot(freqs, abs(squeeze(freqresp(Gs{i}('Fsy', 'fy'), freqs, 'Hz'))));
  end
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude [m/m]');xlabel('Frequency [Hz]');

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

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

#+name: fig:amplitifed_piezo_xy_rotation_plant_iff
#+caption: Transfer function matrix from forces to force sensors for multiple rotation speed
#+RESULTS:
[[file:figs/amplitifed_piezo_xy_rotation_plant_iff.png]]

** Root Locus
#+begin_src matlab :exports none :post
  figure;

  gains = logspace(1, 6, 500);

  hold on;
  for i = 1:length(Ws)
      set(gca,'ColorOrderIndex',i);
      plot(real(pole(Gs{i})),  imag(pole(Gs{i})), 'x');
      set(gca,'ColorOrderIndex',i);
      plot(real(tzero(Gs{i})),  imag(tzero(Gs{i})), 'o');
      for k = 1:length(gains)
          set(gca,'ColorOrderIndex',i);
          cl_poles = pole(feedback(Gs{i}, -gains(k)/s*eye(2)));
          plot(real(cl_poles), imag(cl_poles), '.');
      end
  end
  hold off;
  axis square;
  xlim([-2900, 100]); ylim([0, 3000]);

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

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

#+name: fig:amplified_piezo_xy_rotation_root_locus
#+caption: Root locus for 3 rotating speed
#+RESULTS:
[[file:figs/amplified_piezo_xy_rotation_root_locus.png]]

** Analysis
The negative stiffness induced by the rotation is equal to $m \omega_0^2$.
Thus, the maximum rotation speed where IFF can be applied is:
\[ \omega_\text{max} = \sqrt{\frac{k_1}{m}} \approx 156\,[Hz] \]

Let's verify that.
#+begin_src matlab
  Ws = 2*pi*[140, 160];
#+end_src

#+begin_src matlab :exports none
  Gs = {zeros(length(Ws), 1)};
#+end_src

Identification
#+begin_src matlab
  %% Name of the Simulink File
  mdl = 'amplified_piezo_xy_rotating_stage';

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

  io(io_i) = linio([mdl, '/Fs'], 1, 'openoutput'); io_i = io_i + 1;
  io(io_i) = linio([mdl, '/Fs'], 2, 'openoutput'); io_i = io_i + 1;

  for i = 1:length(Ws)
      ws = Ws(i);
      G = linearize(mdl, io, 0);
      G.InputName  = {'fx', 'fy'};
      G.OutputName = {'Fsx', 'Fsy'};
      Gs(i) = {G};
  end
#+end_src

#+begin_src matlab :exports none
  figure;

  gains = logspace(1, 6, 500);

  hold on;
  for i = 1:length(Ws)
      set(gca,'ColorOrderIndex',i);
      plot(real(pole(Gs{i})),  imag(pole(Gs{i})), 'x');
      set(gca,'ColorOrderIndex',i);
      plot(real(tzero(Gs{i})),  imag(tzero(Gs{i})), 'o');
      for k = 1:length(gains)
          set(gca,'ColorOrderIndex',i);
          cl_poles = pole(feedback(Gs{i}, -gains(k)/s*eye(2)));
          plot(real(cl_poles), imag(cl_poles), '.');
      end
  end
  hold off;
  axis square;
  xlim([-100, 50]); ylim([0, 150]);

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

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

#+name: fig:amplified_piezo_xy_rotating_unstable_root_locus
#+caption: Root Locus for the two considered rotation speed. For the red curve, the system is unstable.
#+RESULTS:
[[file:figs/amplified_piezo_xy_rotating_unstable_root_locus.png]]

* Stewart Platform with Amplified Actuators
** 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
  simulinkproject('../');
#+end_src

#+begin_src matlab
  open('nass_model.slx')
#+end_src