Add analysis for voice coil actuators (uniaxial model)

uniaxial/index.org

@@ -47,6 +47,8 @@ The idea is to use the same model as the full Simscape Model but to restrict the
 This is done in order to more easily study the system and evaluate control techniques.
 
 * Simscape Model
+<<sec:simscape_model>>
+
 A schematic of the uniaxial model used for simulations is represented in figure [[fig:uniaxial-model-nass-flexible]].
 
 The perturbations $w$ are:
@@ -623,6 +625,7 @@ Schematics of the active damping techniques are displayed in figure [[fig:uniaxi
 [[file:figs/uniaxial-model-nass-flexible-active-damping.png]]
 
 * Undamped System
+<<sec:undamped>>
 ** Introduction                                                     :ignore:
 Let's start by study the undamped system.
 
@@ -2184,6 +2187,7 @@ And initialize the controllers.
 Direct Velocity Feedback:
 #+end_important
 * With Cedrat Piezo-electric Actuators
+<<sec:cedrat_actuator>>
 ** 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>>
@@ -2749,3 +2753,400 @@ It is important to note that the effect of direct forces applied to the sample a
 | Sensitivity ($F_s$)       | - (at low freq) | +               | +        |
 | Sensitivity ($F_{ty,rz}$) | +               | -               | +        |
 | Overall RMS of $D$        | =               | =               | =        |
+* Voice Coil
+<<sec:voice_coil>>
+** 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 'simscape/sim_nano_station_uniaxial.slx'
+#+end_src
+
+** Init
+We initialize all the stages with the default parameters.
+The nano-hexapod is an hexapod with voice coils and the sample has a mass of 50kg.
+
+#+begin_src matlab :exports none
+  initializeGround();
+  initializeGranite();
+  initializeTy();
+  initializeRy();
+  initializeRz();
+  initializeMicroHexapod();
+  initializeAxisc();
+  initializeMirror();
+  initializeNanoHexapod(struct('actuator', 'lorentz'));
+  initializeSample(struct('mass', 50));
+#+end_src
+
+All the controllers are set to 0 (Open Loop).
+#+begin_src matlab :exports none
+  K = tf(0);
+  save('./mat/controllers.mat', 'K', '-append');
+  K_iff = tf(0);
+  save('./mat/controllers.mat', 'K_iff', '-append');
+  K_rmc = tf(0);
+  save('./mat/controllers.mat', 'K_rmc', '-append');
+  K_dvf = tf(0);
+  save('./mat/controllers.mat', 'K_dvf', '-append');
+#+end_src
+
+** Identification
+We identify the dynamics of the system.
+#+begin_src matlab
+  %% Options for Linearized
+  options = linearizeOptions;
+  options.SampleTime = 0;
+
+  %% Name of the Simulink File
+  mdl = 'sim_nano_station_uniaxial';
+#+end_src
+
+The inputs and outputs are defined below and corresponds to the name of simulink blocks.
+#+begin_src matlab
+  %% Input/Output definition
+  io(1)  = linio([mdl, '/Dw'],   1, 'input'); % Ground Motion
+  io(2)  = linio([mdl, '/Fs'],   1, 'input'); % Force applied on the sample
+  io(3)  = linio([mdl, '/Fnl'],  1, 'input'); % Force applied by the NASS
+  io(4)  = linio([mdl, '/Fdty'], 1, 'input'); % Parasitic force Ty
+  io(5)  = linio([mdl, '/Fdrz'], 1, 'input'); % Parasitic force Rz
+
+  io(6)  = linio([mdl, '/Dsm'],  1, 'output'); % Displacement of the sample
+  io(7)  = linio([mdl, '/Fnlm'], 1, 'output'); % Force sensor in NASS's legs
+  io(8)  = linio([mdl, '/Dnlm'], 1, 'output'); % Displacement of NASS's legs
+  io(9)  = linio([mdl, '/Dgm'],  1, 'output'); % Absolute displacement of the granite
+  io(10) = linio([mdl, '/Vlm'],  1, 'output'); % Measured absolute velocity of the top NASS platform
+#+end_src
+
+Finally, we use the =linearize= Matlab function to extract a state space model from the simscape model.
+#+begin_src matlab
+  %% Run the linearization
+  G_vc = linearize(mdl, io, options);
+  G_vc.InputName  = {'Dw',   ... % Ground Motion [m]
+                     'Fs',   ... % Force Applied on Sample [N]
+                     'Fn',   ... % Force applied by NASS [N]
+                     'Fty',  ... % Parasitic Force Ty [N]
+                     'Frz'};     % Parasitic Force Rz [N]
+  G_vc.OutputName = {'D',    ... % Measured sample displacement x.r.t. granite [m]
+                     'Fnm',  ... % Force Sensor in NASS [N]
+                     'Dnm',  ... % Displacement Sensor in NASS [m]
+                     'Dgm',  ... % Asbolute displacement of Granite [m]
+                     'Vlm'}; ... % Absolute Velocity of NASS [m/s]
+#+end_src
+
+Finally, we save the identified system dynamics for further analysis.
+#+begin_src matlab
+  save('./uniaxial/mat/plants.mat', 'G_vc', '-append');
+#+end_src
+
+** Sensitivity to Disturbances
+We load the dynamics when using a piezo-electric nano hexapod to compare the results.
+#+begin_src matlab
+  load('./uniaxial/mat/plants.mat', 'G');
+#+end_src
+
+We show several plots representing the sensitivity to disturbances:
+- in figure [[fig:uniaxial-sensitivity-vc-disturbances]] the transfer functions from ground motion $D_w$ to the sample position $D$ and the transfer function from direct force on the sample $F_s$ to the sample position $D$ are shown
+- in figure [[fig:uniaxial-sensitivity-vc-force-dist]], it is the effect of parasitic forces of the positioning stages ($F_{ty}$ and $F_{rz}$) on the position $D$ of the sample that are shown
+
+#+begin_src matlab :exports none
+  freqs = logspace(0, 3, 1000);
+
+  figure;
+
+  subplot(2, 1, 1);
+  title('$D_w$ to $D$');
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G_vc('D',   'Dw'), freqs, 'Hz'))), 'k-', 'DisplayName', 'G - VC');
+  plot(freqs, abs(squeeze(freqresp(G('D',   'Dw'), freqs, 'Hz'))), 'k--', 'DisplayName', 'G - PZ');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  legend('location', 'northeast');
+  ylabel('Amplitude [m/m]'); xlabel('Frequency [Hz]');
+
+  subplot(2, 1, 2);
+  title('$F_s$ to $D$');
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G_vc('D', 'Fs'), freqs, 'Hz'))), 'k-');
+  plot(freqs, abs(squeeze(freqresp(G('D', 'Fs'), freqs, 'Hz'))), 'k--');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/uniaxial-sensitivity-vc-disturbances.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:uniaxial-sensitivity-vc-disturbances
+#+CAPTION: Sensitivity to disturbances ([[./figs/uniaxial-sensitivity-vc-disturbances.png][png]], [[./figs/uniaxial-sensitivity-vc-disturbances.pdf][pdf]])
+[[file:figs/uniaxial-sensitivity-vc-disturbances.png]]
+
+#+begin_src matlab :exports none
+  freqs = logspace(0, 3, 1000);
+
+  figure;
+
+  subplot(2, 1, 1);
+  title('$F_{ty}$ to $D$');
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G_vc('D', 'Fty'), freqs, 'Hz'))), 'k-', 'DisplayName', 'G - VC');
+  plot(freqs, abs(squeeze(freqresp(G('D', 'Fty'), freqs, 'Hz'))), 'k--', 'DisplayName', 'G - PZ');
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  legend('location', 'northeast');
+  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
+
+  subplot(2, 1, 2);
+  title('$F_{rz}$ to $D$');
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G_vc('D', 'Frz'), freqs, 'Hz'))), 'k-');
+  plot(freqs, abs(squeeze(freqresp(G('D', 'Frz'), freqs, 'Hz'))), 'k--');
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/uniaxial-sensitivity-vc-force-dist.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:uniaxial-sensitivity-vc-force-dist
+#+CAPTION: Sensitivity to disturbances ([[./figs/uniaxial-sensitivity-vc-force-dist.png][png]], [[./figs/uniaxial-sensitivity-vc-force-dist.pdf][pdf]])
+[[file:figs/uniaxial-sensitivity-vc-force-dist.png]]
+
+** Noise Budget
+We first load the measured PSD of the disturbance.
+#+begin_src matlab
+  load('./disturbances/mat/dist_psd.mat', 'dist_f');
+#+end_src
+
+The effect of these disturbances on the distance $D$ is computed below.
+#+begin_src matlab :exports none
+  % Power Spectral Density of the relative Displacement [m^2/Hz]
+  psd_vc_gm_d = dist_f.psd_gm.*abs(squeeze(freqresp(G_vc('D', 'Dw'),  dist_f.f, 'Hz'))).^2;
+  psd_vc_ty_d = dist_f.psd_ty.*abs(squeeze(freqresp(G_vc('D', 'Fty'), dist_f.f, 'Hz'))).^2;
+  psd_vc_rz_d = dist_f.psd_rz.*abs(squeeze(freqresp(G_vc('D', 'Frz'), dist_f.f, 'Hz'))).^2;
+#+end_src
+
+#+begin_src matlab :exports none
+  % Power Spectral Density of the relative Displacement [m^2/Hz]
+  psd_gm_d = dist_f.psd_gm.*abs(squeeze(freqresp(G('D', 'Dw'),  dist_f.f, 'Hz'))).^2;
+  psd_ty_d = dist_f.psd_ty.*abs(squeeze(freqresp(G('D', 'Fty'), dist_f.f, 'Hz'))).^2;
+  psd_rz_d = dist_f.psd_rz.*abs(squeeze(freqresp(G('D', 'Frz'), dist_f.f, 'Hz'))).^2;
+#+end_src
+
+The PSD of the obtain distance $D$ due to each of the perturbation is shown in figure [[fig:uniaxial-vc-psd-dist]] and the Cumulative Amplitude Spectrum is shown in figure [[fig:uniaxial-vc-cas-dist]].
+
+The Root Mean Square value of the obtained displacement $D$ is computed below and can be determined from the figure [[fig:uniaxial-vc-cas-dist]].
+#+begin_src matlab :results value replace :exports results
+  cas_tot_d = sqrt(cumtrapz(dist_f.f, psd_vc_rz_d+psd_vc_ty_d+psd_vc_gm_d)); cas_tot_d(end)
+#+end_src
+
+#+RESULTS:
+: 4.8793e-06
+
+#+begin_src matlab :exports none
+  freqs = logspace(0, 3, 1000);
+
+  figure;
+  hold on;
+  plot(dist_f.f, psd_vc_gm_d, 'DisplayName', '$D_w$');
+  plot(dist_f.f, psd_vc_ty_d, 'DisplayName', '$T_y$');
+  plot(dist_f.f, psd_vc_rz_d, 'DisplayName', '$R_z$');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('CAS of the effect of disturbances on $D$ $\left[\frac{m^2}{Hz}\right]$');  xlabel('Frequency [Hz]');
+  legend('location', 'northeast')
+  xlim([0.5, 500]);
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/uniaxial-vc-psd-dist.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:uniaxial-vc-psd-dist
+#+CAPTION: PSD of the displacement $D$ due to disturbances ([[./figs/uniaxial-vc-psd-dist.png][png]], [[./figs/uniaxial-vc-psd-dist.pdf][pdf]])
+[[file:figs/uniaxial-vc-psd-dist.png]]
+
+#+begin_src matlab :exports none
+  freqs = logspace(0, 3, 1000);
+
+  figure;
+  hold on;
+  plot(dist_f.f, flip(sqrt(-cumtrapz(flip(dist_f.f), flip(psd_vc_gm_d)))), 'DisplayName', '$D_w$ - VC');
+  plot(dist_f.f, flip(sqrt(-cumtrapz(flip(dist_f.f), flip(psd_vc_ty_d)))), 'DisplayName', '$T_y$ - VC');
+  plot(dist_f.f, flip(sqrt(-cumtrapz(flip(dist_f.f), flip(psd_vc_rz_d)))), 'DisplayName', '$R_z$ - VC');
+  plot(dist_f.f, flip(sqrt(-cumtrapz(flip(dist_f.f), flip(psd_vc_gm_d+psd_vc_ty_d+psd_vc_rz_d)))), 'k-', 'DisplayName', 'tot - VC');
+  plot(dist_f.f, flip(sqrt(-cumtrapz(flip(dist_f.f), flip(psd_gm_d+psd_ty_d+psd_rz_d)))), 'k--', 'DisplayName', 'tot - PZ');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('CAS of the effect of disturbances on $D$ [m]');  xlabel('Frequency [Hz]');
+  legend('location', 'northeast')
+  xlim([0.5, 500]); ylim([1e-12, 5e-6]);
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/uniaxial-vc-cas-dist.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:uniaxial-vc-cas-dist
+#+CAPTION: CAS of the displacement $D$ due the disturbances ([[./figs/uniaxial-vc-cas-dist.png][png]], [[./figs/uniaxial-vc-cas-dist.pdf][pdf]])
+[[file:figs/uniaxial-vc-cas-dist.png]]
+
+#+begin_important
+  Even though the RMS value of the displacement $D$ is lower when using a piezo-electric actuator, the motion is mainly due to high frequency disturbances which are more difficult to control (an higher control bandwidth is required).
+
+  Thus, it may be desirable to use voice coil actuators.
+#+end_important
+** Integral Force Feedback
+#+begin_src matlab
+  K_iff = -20/s;
+#+end_src
+
+#+begin_src matlab :exports none
+  freqs = logspace(-1, 2, 1000);
+
+  figure;
+
+  ax1 = subplot(2, 1, 1);
+  plot(freqs, abs(squeeze(freqresp(K_iff*G_vc('Fnm', 'Fn'), freqs, 'Hz'))), 'k-');
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]);
+
+  ax2 = subplot(2, 1, 2);
+  plot(freqs, 180/pi*angle(squeeze(freqresp(K_iff*G_vc('Fnm', 'Fn'), freqs, 'Hz'))), 'k-');
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
+  ylim([-180, 180]);
+  yticks([-180, -90, 0, 90, 180]);
+
+  linkaxes([ax1,ax2],'x');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/uniaxial_iff_vc_open_loop.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:uniaxial_iff_vc_open_loop
+#+CAPTION: Open Loop Transfer Function for IFF control when using a voice coil actuator ([[./figs/uniaxial_iff_vc_open_loop.png][png]], [[./figs/uniaxial_iff_vc_open_loop.pdf][pdf]])
+[[file:figs/uniaxial_iff_vc_open_loop.png]]
+
+** Identification of the Damped Plant
+Let's initialize the system prior to identification.
+#+begin_src matlab
+  initializeGround();
+  initializeGranite();
+  initializeTy();
+  initializeRy();
+  initializeRz();
+  initializeMicroHexapod();
+  initializeAxisc();
+  initializeMirror();
+  initializeNanoHexapod(struct('actuator', 'lorentz'));
+  initializeSample(struct('mass', 50));
+#+end_src
+
+All the controllers are set to 0.
+#+begin_src matlab
+  K = tf(0);
+  save('./mat/controllers.mat', 'K', '-append');
+  K_iff = -K_iff;
+  save('./mat/controllers.mat', 'K_iff', '-append');
+  K_rmc = tf(0);
+  save('./mat/controllers.mat', 'K_rmc', '-append');
+  K_dvf = tf(0);
+  save('./mat/controllers.mat', 'K_dvf', '-append');
+#+end_src
+
+#+begin_src matlab
+  %% Options for Linearized
+  options = linearizeOptions;
+  options.SampleTime = 0;
+
+  %% Name of the Simulink File
+  mdl = 'sim_nano_station_uniaxial';
+#+end_src
+
+#+begin_src matlab
+  %% Input/Output definition
+  io(1)  = linio([mdl, '/Dw'],    1, 'input');  % Ground Motion
+  io(2)  = linio([mdl, '/Fs'],    1, 'input');  % Force applied on the sample
+  io(3)  = linio([mdl, '/Fnl'],   1, 'input');  % Force applied by the NASS
+  io(4)  = linio([mdl, '/Fdty'],  1, 'input');  % Parasitic force Ty
+  io(5)  = linio([mdl, '/Fdrz'],  1, 'input');  % Parasitic force Rz
+
+  io(6)  = linio([mdl, '/Dsm'],  1, 'output'); % Displacement of the sample
+  io(7)  = linio([mdl, '/Fnlm'], 1, 'output'); % Force sensor in NASS's legs
+  io(8)  = linio([mdl, '/Dnlm'], 1, 'output'); % Displacement of NASS's legs
+  io(9)  = linio([mdl, '/Dgm'],  1, 'output'); % Absolute displacement of the granite
+  io(10) = linio([mdl, '/Vlm'],  1, 'output'); % Measured absolute velocity of the top NASS platform
+#+end_src
+
+#+begin_src matlab
+  %% Run the linearization
+  G_vc_iff = linearize(mdl, io, options);
+  G_vc_iff.InputName  = {'Dw',   ... % Ground Motion [m]
+                         'Fs',   ... % Force Applied on Sample [N]
+                         'Fn',   ... % Force applied by NASS [N]
+                         'Fty',  ... % Parasitic Force Ty [N]
+                         'Frz'};     % Parasitic Force Rz [N]
+  G_vc_iff.OutputName = {'D',    ... % Measured sample displacement x.r.t. granite [m]
+                         'Fnm',  ... % Force Sensor in NASS [N]
+                         'Dnm',  ... % Displacement Sensor in NASS [m]
+                         'Dgm',  ... % Asbolute displacement of Granite [m]
+                         'Vlm'}; ... % Absolute Velocity of NASS [m/s]
+#+end_src
+
+** Noise Budget
+We compute the obtain PSD of the displacement $D$ when using IFF.
+#+begin_src matlab :exports none
+  % Power Spectral Density of the relative Displacement [m^2/Hz]
+  psd_vc_iff_gm_d = dist_f.psd_gm.*abs(squeeze(freqresp(G_vc_iff('D', 'Dw'),  dist_f.f, 'Hz'))).^2;
+  psd_vc_iff_ty_d = dist_f.psd_ty.*abs(squeeze(freqresp(G_vc_iff('D', 'Fty'), dist_f.f, 'Hz'))).^2;
+  psd_vc_iff_rz_d = dist_f.psd_rz.*abs(squeeze(freqresp(G_vc_iff('D', 'Frz'), dist_f.f, 'Hz'))).^2;
+#+end_src
+
+#+begin_src matlab :exports none
+  freqs = logspace(0, 3, 1000);
+
+  figure;
+  hold on;
+  plot(dist_f.f, flip(sqrt(-cumtrapz(flip(dist_f.f), flip(psd_gm_d+psd_ty_d+psd_rz_d)))), '-', 'DisplayName', 'OL - PZ');
+  plot(dist_f.f, flip(sqrt(-cumtrapz(flip(dist_f.f), flip(psd_vc_gm_d+psd_vc_ty_d+psd_vc_rz_d)))), 'k-', 'DisplayName', 'OL - VC');
+  plot(dist_f.f, flip(sqrt(-cumtrapz(flip(dist_f.f), flip(psd_vc_iff_gm_d+psd_vc_iff_ty_d+psd_vc_iff_rz_d)))), 'k--', 'DisplayName', 'IFF - VC');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('CAS of the effect of disturbances on $D$ [m]');  xlabel('Frequency [Hz]');
+  legend('location', 'northeast')
+  xlim([0.5, 500]); ylim([1e-12, 5e-6]);
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/uniaxial-cas-iff-vc.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:uniaxial-cas-iff-vc
+#+CAPTION: CAS of the displacement $D$ ([[./figs/uniaxial-cas-iff-vc.png][png]], [[./figs/uniaxial-cas-iff-vc.pdf][pdf]])
+[[file:figs/uniaxial-cas-iff-vc.png]]
+
+** Conclusion
+#+begin_important
+  The use of voice coil actuators would allow a better disturbance rejection for a fixed bandwidth compared with a piezo-electric hexapod.
+
+  Similarly, it would require much lower bandwidth to attain the same level of disturbance rejection for $D$.
+#+end_important