Publications: add paper pages (dehaeze18, brumund21, dehaeze20, dehaeze21 x2), drop Fastjack, look for PDFs in journal/

2026-09-27 19:45:28 +02:00
parent 8b90d69df2
commit 7195b9e67a
102 changed files with 5352 additions and 41 deletions
@@ -0,0 +1,959 @@
+++
+title = "Active Damping of Rotating Platforms using Integral Force Feedback - Matlab Computation"
+author = ["Dehaeze Thomas"]
+draft = false
+++
+
+<hr>
+<p>This report is also available as a <a href="./index.pdf">pdf</a>.</p>
+<hr>
+
+This document gathers the Matlab code used to for the conference paper (<a href="#citeproc_bib_item_1">Dehaeze and Collette 2020</a>) and the journal paper (<a href="#citeproc_bib_item_3">Dehaeze and Collette 2021</a>).
+
+It is structured in several sections:
+
+-   Section : presents a simple model of a rotating suspended platform that will be used throughout this study.
+-   Section : explains how the unconditional stability of IFF is lost due to Gyroscopic effects induced by the rotation.
+-   Section : suggests a simple modification of the control law such that damping can be added to the suspension modes in a robust way.
+-   Section : proposes to add springs in parallel with the force sensors to regain the unconditional stability of IFF.
+-   Section : compares both proposed modifications to the classical IFF in terms of damping authority and closed-loop system behavior.
+-   Section : contains the notations used for both the Matlab code and the paper
+
+The matlab code is accessible on [Zonodo](https://zenodo.org/record/3894343) and [Github](https://github.com/tdehaeze/dehaeze20_contr_stewa_platf) (<a href="#citeproc_bib_item_2">Dehaeze 2020</a>). It can also be download as a `.zip` file [here](https://git.tdehaeze.xyz/tdehaeze/dehaeze21_activ_dampin_rotat_platf_using/archive/master.zip).
+
+To run the Matlab code, go in the `matlab` directory and run the following Matlab files corresponding to each section.
+
+<div class="table-caption">
+  <span class="table-number">Table 1:</span>
+  Paper's sections and corresponding Matlab files
+</div>
+
+| Sections | Matlab File                |
+|----------|----------------------------|
+| Section  | `s1_system_description.m`  |
+| Section  | `s2_iff_pure_int.m`        |
+| Section  | `s3_iff_hpf.m`             |
+| Section  | `s4_iff_kp.m`              |
+| Section  | `s5_act_damp_comparison.m` |
+
+
+## System Description and Analysis {#system-description-and-analysis}
+
+<span class="org-target" id="org-target--sec-system-description"></span>
+
+
+### System description {#system-description}
+
+The system consists of one 2 degree of freedom translation stage on top of a spindle (figure [Figure 1](#figure--fig:system)).
+
+<a id="figure--fig:system"></a>
+
+{{< figure src="figs-paper/system.png" caption="<span class='figure-number'>Figure 1: </span>Schematic of the studied system" >}}
+
+The control inputs are the forces applied by the actuators of the translation stage (\\(F\_u\\) and \\(F\_v\\)).
+As the translation stage is rotating around the Z axis due to the spindle, the forces are applied along \\(\vec{i}\_u\\) and \\(\vec{i}\_v\\).
+
+
+### Equations {#equations}
+
+Based on the Figure [Figure 1](#figure--fig:system), the equations of motions are:
+
+<div class="important">
+
+\begin{equation}
+\begin{bmatrix} d\_u \\\ d\_v \end{bmatrix} =
+\bm{G}\_d
+\begin{bmatrix} F\_u \\\ F\_v \end{bmatrix}
+\end{equation}
+
+Where \\(\bm{G}\_d\\) is a \\(2 \times 2\\) transfer function matrix.
+
+\begin{equation}
+\bm{G}\_d = \frac{1}{k} \frac{1}{G\_{dp}}
+\begin{bmatrix}
+   G\_{dz} & G\_{dc} \\\\
+  -G\_{dc} & G\_{dz}
+\end{bmatrix}
+\end{equation}
+
+With:
+
+\begin{align}
+  G\_{dp} &= \left( \frac{s^2}{{\omega\_0}^2} + 2 \xi \frac{s}{\omega\_0} + 1 - \frac{{\Omega}^2}{{\omega\_0}^2} \right)^2 + \left( 2 \frac{\Omega}{\omega\_0} \frac{s}{\omega\_0} \right)^2 \\\\
+  G\_{dz} &= \frac{s^2}{{\omega\_0}^2} + 2 \xi \frac{s}{\omega\_0} + 1 - \frac{{\Omega}^2}{{\omega\_0}^2} \\\\
+  G\_{dc} &= 2 \frac{\Omega}{\omega\_0} \frac{s}{\omega\_0}
+\end{align}
+
+</div>
+
+
+### Numerical Values {#numerical-values}
+
+Let's define initial values for the model.
+
+```matlab
+  k = 1;    % Actuator Stiffness [N/m]
+  c = 0.05; % Actuator Damping [N/(m/s)]
+  m = 1;    % Payload mass [kg]
+```
+
+```matlab
+  xi = c/(2*sqrt(k*m));
+  w0 = sqrt(k/m); % [rad/s]
+```
+
+
+### Campbell Diagram {#campbell-diagram}
+
+The Campbell Diagram displays the evolution of the real and imaginary parts of the system as a function of the rotating speed.
+
+It is shown in Figures [Figure 2](#figure--fig:campbell-diagram-real) and [Figure 3](#figure--fig:campbell-diagram-imag), and one can see that the system becomes unstable for \\(\Omega > \omega\_0\\) (the real part of one of the poles becomes positive).
+
+<a id="figure--fig:campbell-diagram-real"></a>
+
+{{< figure src="figs/campbell_diagram_real.png" caption="<span class='figure-number'>Figure 2: </span>Campbell Diagram - Real Part" >}}
+
+<a id="figure--fig:campbell-diagram-imag"></a>
+
+{{< figure src="figs/campbell_diagram_imag.png" caption="<span class='figure-number'>Figure 3: </span>Campbell Diagram - Imaginary Part" >}}
+
+
+### Simscape Model {#simscape-model}
+
+In order to validate all the equations of motion, a Simscape model of the same system has been developed.
+The dynamics of the system can be identified from the Simscape model and compare with the analytical model.
+
+The rotating speed for the Simscape Model is defined.
+
+```matlab
+  W = 0.1; % Rotation Speed [rad/s]
+```
+
+```matlab
+  open('rotating_frame.slx');
+```
+
+The transfer function from \\([F\_u, F\_v]\\) to \\([d\_u, d\_v]\\) is identified from the Simscape model.
+
+```matlab
+  %% Name of the Simulink File
+  mdl = 'rotating_frame';
+
+  %% Input/Output definition
+  clear io; io_i = 1;
+  io(io_i) = linio([mdl, '/K'], 1, 'openinput');  io_i = io_i + 1;
+  io(io_i) = linio([mdl, '/G'], 2, 'openoutput'); io_i = io_i + 1;
+```
+
+```matlab
+  G = linearize(mdl, io, 0);
+
+  %% Input/Output definition
+  G.InputName  = {'Fu', 'Fv'};
+  G.OutputName = {'du', 'dv'};
+```
+
+The same transfer function from \\([F\_u, F\_v]\\) to \\([d\_u, d\_v]\\) is written down from the analytical model.
+
+```matlab
+  Gth = (1/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ...
+        [(s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2), 2*W*s/(w0^2) ; ...
+         -2*W*s/(w0^2), (s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)];
+```
+
+Both transfer functions are compared in Figure [Figure 4](#figure--fig:plant-simscape-analytical) and are found to perfectly match.
+
+<a id="figure--fig:plant-simscape-analytical"></a>
+
+{{< figure src="figs/plant_simscape_analytical.png" caption="<span class='figure-number'>Figure 4: </span>Bode plot of the transfer function from \\([F\_u, F\_v]\\) to \\([d\_u, d\_v]\\) as identified from the Simscape model and from an analytical model" >}}
+
+
+### Effect of the rotation speed {#effect-of-the-rotation-speed}
+
+The transfer functions from \\([F\_u, F\_v]\\) to \\([d\_u, d\_v]\\) are identified for the following rotating speeds.
+
+```matlab
+  Ws = [0, 0.2, 0.7, 1.1]*w0; % Rotating Speeds [rad/s]
+```
+
+```matlab
+  Gs = {zeros(2, 2, length(Ws))};
+
+  for W_i = 1:length(Ws)
+      W = Ws(W_i);
+
+      Gs(:, :, W_i) = {(1/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ...
+                       [(s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2), 2*W*s/(w0^2) ; ...
+                        -2*W*s/(w0^2), (s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)]};
+  end
+```
+
+They are compared in Figures [Figure 5](#figure--fig:plant-compare-rotating-speed-direct) and [Figure 6](#figure--fig:plant-compare-rotating-speed-coupling).
+
+<a id="figure--fig:plant-compare-rotating-speed-direct"></a>
+
+{{< figure src="figs/plant_compare_rotating_speed_direct.png" caption="<span class='figure-number'>Figure 5: </span>Comparison of the transfer functions from \\([F\_u, F\_v]\\) to \\([d\_u, d\_v]\\) for several rotating speed - Direct Terms" >}}
+
+<a id="figure--fig:plant-compare-rotating-speed-coupling"></a>
+
+{{< figure src="figs/plant_compare_rotating_speed_coupling.png" caption="<span class='figure-number'>Figure 6: </span>Comparison of the transfer functions from \\([F\_u, F\_v]\\) to \\([d\_u, d\_v]\\) for several rotating speed - Coupling Terms" >}}
+
+
+## Problem with pure Integral Force Feedback {#problem-with-pure-integral-force-feedback}
+
+<span class="org-target" id="org-target--sec-iff-pure-int"></span>
+
+Force sensors are added in series with the two actuators (Figure [Figure 7](#figure--fig:system-iff)).
+
+Two identical controllers \\(K\_F\\) are used to feedback each of the sensed force to its associated actuator.
+
+<a id="figure--fig:system-iff"></a>
+
+{{< figure src="figs-paper/system_iff.png" caption="<span class='figure-number'>Figure 7: </span>System with added Force Sensor in series with the actuators" >}}
+
+
+### Plant Parameters {#plant-parameters}
+
+Let's define initial values for the model.
+
+```matlab
+  k = 1;    % Actuator Stiffness [N/m]
+  c = 0.05; % Actuator Damping [N/(m/s)]
+  m = 1;    % Payload mass [kg]
+```
+
+```matlab
+  xi = c/(2*sqrt(k*m));
+  w0 = sqrt(k/m); % [rad/s]
+```
+
+
+### Equations {#equations}
+
+The sensed forces are equal to:
+
+\begin{equation}
+\begin{bmatrix} f\_{u} \\\ f\_{v} \end{bmatrix} =
+\begin{bmatrix}
+  1 & 0 \\\\
+  0 & 1
+\end{bmatrix}
+\begin{bmatrix} F\_u \\\ F\_v \end{bmatrix} - (c s + k)
+\begin{bmatrix} d\_u \\\ d\_v \end{bmatrix}
+\end{equation}
+
+Which then gives:
+
+<div class="important">
+
+\begin{equation}
+\begin{bmatrix} f\_{u} \\\ f\_{v} \end{bmatrix} =
+\bm{G}\_{f}
+\begin{bmatrix} F\_u \\\ F\_v \end{bmatrix}
+\end{equation}
+
+\begin{equation}
+\begin{bmatrix} f\_{u} \\\ f\_{v} \end{bmatrix} =
+\frac{1}{G\_{fp}}
+\begin{bmatrix}
+  G\_{fz} & -G\_{fc} \\\\
+  G\_{fc} &  G\_{fz}
+\end{bmatrix}
+\begin{bmatrix} F\_u \\\ F\_v \end{bmatrix}
+\end{equation}
+
+\begin{align}
+  G\_{fp} &= \left( \frac{s^2}{{\omega\_0}^2} + 2 \xi \frac{s}{\omega\_0} + 1 - \frac{{\Omega}^2}{{\omega\_0}^2} \right)^2 + \left( 2 \frac{\Omega}{\omega\_0} \frac{s}{\omega\_0} \right)^2 \\\\
+  G\_{fz} &= \left( \frac{s^2}{{\omega\_0}^2} - \frac{\Omega^2}{{\omega\_0}^2} \right) \left( \frac{s^2}{{\omega\_0}^2} + 2 \xi \frac{s}{\omega\_0} + 1 - \frac{{\Omega}^2}{{\omega\_0}^2} \right) + \left( 2 \frac{\Omega}{\omega\_0} \frac{s}{\omega\_0} \right)^2 \\\\
+  G\_{fc} &= \left( 2 \xi \frac{s}{\omega\_0} + 1 \right) \left( 2 \frac{\Omega}{\omega\_0} \frac{s}{\omega\_0} \right)
+\end{align}
+
+</div>
+
+
+### Comparison of the Analytical Model and the Simscape Model {#comparison-of-the-analytical-model-and-the-simscape-model}
+
+The rotation speed is set to \\(\Omega = 0.1 \omega\_0\\).
+
+```matlab
+  W = 0.1*w0; % [rad/s]
+```
+
+```matlab
+  open('rotating_frame.slx');
+```
+
+And the transfer function from \\([F\_u, F\_v]\\) to \\([f\_u, f\_v]\\) is identified using the Simscape model.
+
+```matlab
+  %% Name of the Simulink File
+  mdl = 'rotating_frame';
+
+  %% Input/Output definition
+  clear io; io_i = 1;
+  io(io_i) = linio([mdl, '/K'], 1, 'openinput');  io_i = io_i + 1;
+  io(io_i) = linio([mdl, '/G'], 1, 'openoutput'); io_i = io_i + 1;
+```
+
+```matlab
+  Giff = linearize(mdl, io, 0);
+
+  %% Input/Output definition
+  Giff.InputName  = {'Fu', 'Fv'};
+  Giff.OutputName = {'fu', 'fv'};
+```
+
+The same transfer function from \\([F\_u, F\_v]\\) to \\([f\_u, f\_v]\\) is written down from the analytical model.
+
+```matlab
+  Giff_th = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ...
+            [(s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)) + (2*W*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*W*s/(w0^2) ; ...
+             (2*xi*s/w0 + 1)*2*W*s/(w0^2), (s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))+ (2*W*s/(w0^2))^2];
+```
+
+The two are compared in Figure [Figure 8](#figure--fig:plant-iff-comp-simscape-analytical) and found to perfectly match.
+
+<a id="figure--fig:plant-iff-comp-simscape-analytical"></a>
+
+{{< figure src="figs/plant_iff_comp_simscape_analytical.png" caption="<span class='figure-number'>Figure 8: </span>Comparison of the transfer functions from \\([F\_u, F\_v]\\) to \\([f\_u, f\_v]\\) between the Simscape model and the analytical one" >}}
+
+
+### Effect of the rotation speed {#effect-of-the-rotation-speed}
+
+The transfer functions from \\([F\_u, F\_v]\\) to \\([f\_u, f\_v]\\) are identified for the following rotating speeds.
+
+```matlab
+  Ws = [0, 0.2, 0.7]*w0; % Rotating Speeds [rad/s]
+```
+
+```matlab
+  Gsiff = {zeros(2, 2, length(Ws))};
+
+  for W_i = 1:length(Ws)
+      W = Ws(W_i);
+
+      Gsiff(:, :, W_i) = {1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ...
+                        [(s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)) + (2*W*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*W*s/(w0^2) ; ...
+                         (2*xi*s/w0 + 1)*2*W*s/(w0^2), (s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))+ (2*W*s/(w0^2))^2]};
+  end
+```
+
+The obtained transfer functions are shown in Figure [Figure 9](#figure--fig:plant-iff-compare-rotating-speed).
+
+<a id="figure--fig:plant-iff-compare-rotating-speed"></a>
+
+{{< figure src="figs/plant_iff_compare_rotating_speed.png" caption="<span class='figure-number'>Figure 9: </span>Comparison of the transfer functions from \\([F\_u, F\_v]\\) to \\([f\_u, f\_v]\\) for several rotating speed" >}}
+
+
+### Decentralized Integral Force Feedback {#decentralized-integral-force-feedback}
+
+The decentralized IFF controller consists of pure integrators:
+
+\begin{equation}
+  \bm{K}\_{\text{IFF}}(s) = \frac{g}{s} \begin{bmatrix}
+    1 & 0 \\\\
+    0 & 1
+  \end{bmatrix}
+\end{equation}
+
+The Root Locus (evolution of the poles of the closed loop system in the complex plane as a function of \\(g\\)) is shown in Figure [Figure 10](#figure--fig:root-locus-pure-iff).
+It is shown that for non-null rotating speed, one pole is bound to the right-half plane, and thus the closed loop system is unstable.
+
+<a id="figure--fig:root-locus-pure-iff"></a>
+
+{{< figure src="figs/root_locus_pure_iff.png" caption="<span class='figure-number'>Figure 10: </span>Root Locus for the Decentralized Integral Force Feedback controller. Several rotating speed are shown." >}}
+
+
+## Integral Force Feedback with an High Pass Filter {#integral-force-feedback-with-an-high-pass-filter}
+
+<span class="org-target" id="org-target--sec-iff-pseudo-int"></span>
+
+
+### Plant Parameters {#plant-parameters}
+
+Let's define initial values for the model.
+
+```matlab
+  k = 1;    % Actuator Stiffness [N/m]
+  c = 0.05; % Actuator Damping [N/(m/s)]
+  m = 1;    % Payload mass [kg]
+```
+
+```matlab
+  xi = c/(2*sqrt(k*m));
+  w0 = sqrt(k/m); % [rad/s]
+```
+
+
+### Modified Integral Force Feedback Controller {#modified-integral-force-feedback-controller}
+
+Let's modify the initial Integral Force Feedback Controller ; instead of using pure integrators, pseudo integrators (i.e. low pass filters) are used:
+
+\begin{equation}
+  K\_{\text{IFF}}(s) = g\frac{1}{\omega\_i + s} \begin{bmatrix}
+  1 & 0 \\\\
+  0 & 1
+\end{bmatrix}
+\end{equation}
+
+where \\(\omega\_i\\) characterize down to which frequency the signal is integrated.
+
+Let's arbitrary choose the following control parameters:
+
+```matlab
+  g = 2;
+  wi = 0.1*w0;
+```
+
+And the following rotating speed.
+
+```matlab
+  Giff = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ...
+          [(s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)) + (2*W*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*W*s/(w0^2) ; ...
+           (2*xi*s/w0 + 1)*2*W*s/(w0^2), (s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))+ (2*W*s/(w0^2))^2];
+```
+
+The obtained Loop Gain is shown in Figure [Figure 11](#figure--fig:loop-gain-modified-iff).
+
+<a id="figure--fig:loop-gain-modified-iff"></a>
+
+{{< figure src="figs/loop_gain_modified_iff.png" caption="<span class='figure-number'>Figure 11: </span>Loop Gain for the modified IFF controller" >}}
+
+
+### Root Locus {#root-locus}
+
+As shown in the Root Locus plot (Figure [Figure 12](#figure--fig:root-locus-modified-iff)), for some value of the gain, the system remains stable.
+
+<a id="figure--fig:root-locus-modified-iff"></a>
+
+{{< figure src="figs/root_locus_modified_iff.png" caption="<span class='figure-number'>Figure 12: </span>Root Locus for the modified IFF controller" >}}
+
+<a id="figure--fig:root-locus-modified-iff-zoom"></a>
+
+{{< figure src="figs/root_locus_modified_iff_zoom.png" caption="<span class='figure-number'>Figure 13: </span>Root Locus for the modified IFF controller - Zoom" >}}
+
+
+### What is the optimal \\(\omega\_i\\) and \\(g\\)? {#what-is-the-optimal-omega-i-and-g}
+
+In order to visualize the effect of \\(\omega\_i\\) on the attainable damping, the Root Locus is displayed in Figure [Figure 14](#figure--fig:root-locus-wi-modified-iff) for the following \\(\omega\_i\\):
+
+```matlab
+  wis = [0.01, 0.1, 0.5, 1]*w0; % [rad/s]
+```
+
+<a id="figure--fig:root-locus-wi-modified-iff"></a>
+
+{{< figure src="figs/root_locus_wi_modified_iff.png" caption="<span class='figure-number'>Figure 14: </span>Root Locus for the modified IFF controller (zoomed plot on the left)" >}}
+
+<a id="figure--fig:root-locus-wi-modified-iff-zoom"></a>
+
+{{< figure src="figs/root_locus_wi_modified_iff_zoom.png" caption="<span class='figure-number'>Figure 15: </span>Root Locus for the modified IFF controller (zoomed plot on the left)" >}}
+
+For the controller
+
+\begin{equation}
+  K\_{\text{IFF}}(s) = g\frac{1}{\omega\_i + s} \begin{bmatrix}
+  1 & 0 \\\\
+  0 & 1
+\end{bmatrix}
+\end{equation}
+
+The gain at which the system becomes unstable is
+
+\begin{equation}
+  g\_\text{max} = \omega\_i \left( \frac{{\omega\_0}^2}{\Omega^2} - 1 \right) \label{eq:iff\_gmax}
+\end{equation}
+
+While it seems that small \\(\omega\_i\\) do allow more damping to be added to the system (Figure [Figure 14](#figure--fig:root-locus-wi-modified-iff)), the control gains may be limited to small values due to \ref{eq:iff\_gmax} thus reducing the attainable damping.
+
+There must be an optimum for \\(\omega\_i\\).
+To find the optimum, the gain that maximize the simultaneous damping of the mode is identified for a wide range of \\(\omega\_i\\) (Figure [Figure 16](#figure--fig:mod-iff-damping-wi)).
+
+```matlab
+  wis = logspace(-2, 1, 100)*w0; % [rad/s]
+
+  opt_xi = zeros(1, length(wis)); % Optimal simultaneous damping
+  opt_gain = zeros(1, length(wis)); % Corresponding optimal gain
+
+  for wi_i = 1:length(wis)
+      wi = wis(wi_i);
+      Kiff = 1/(s + wi)*eye(2);
+
+      fun = @(g)computeSimultaneousDamping(g, Giff, Kiff);
+
+      [g_opt, xi_opt] = fminsearch(fun, 0.5*wi*((w0/W)^2 - 1));
+      opt_xi(wi_i) = 1/xi_opt;
+      opt_gain(wi_i) = g_opt;
+  end
+```
+
+<a id="figure--fig:mod-iff-damping-wi"></a>
+
+{{< figure src="figs/mod_iff_damping_wi.png" caption="<span class='figure-number'>Figure 16: </span>Simultaneous attainable damping of the closed loop poles as a function of \\(\omega\_i\\)" >}}
+
+
+## IFF with a stiffness in parallel with the force sensor {#iff-with-a-stiffness-in-parallel-with-the-force-sensor}
+
+<span class="org-target" id="org-target--sec-iff-parallel-stiffness"></span>
+
+
+### Schematic {#schematic}
+
+In this section additional springs in parallel with the force sensors are added to counteract the negative stiffness induced by the rotation.
+
+<a id="figure--fig:system-parallel-springs"></a>
+
+{{< figure src="figs-paper/system_parallel_springs.png" caption="<span class='figure-number'>Figure 17: </span>Studied system with additional springs in parallel with the actuators and force sensors" >}}
+
+In order to keep the overall stiffness \\(k = k\_a + k\_p\\) constant, a scalar parameter \\(\alpha\\) (\\(0 \le \alpha < 1\\)) is defined to describe the fraction of the total stiffness in parallel with the actuator and force sensor
+
+\begin{equation}
+  k\_p = \alpha k, \quad k\_a = (1 - \alpha) k
+\end{equation}
+
+
+### Equations {#equations}
+
+<div class="important">
+
+\begin{equation}
+\begin{bmatrix} f\_u \\\ f\_v \end{bmatrix} =
+\bm{G}\_k
+\begin{bmatrix} F\_u \\\ F\_v \end{bmatrix}
+\end{equation}
+
+\begin{equation}
+\begin{bmatrix} f\_u \\\ f\_v \end{bmatrix} =
+\frac{1}{G\_{kp}}
+\begin{bmatrix}
+   G\_{kz} & -G\_{kc} \\\\
+   G\_{kc} &  G\_{kz}
+\end{bmatrix}
+\begin{bmatrix} F\_u \\\ F\_v \end{bmatrix}
+\end{equation}
+
+With:
+
+\begin{align}
+  G\_{kp} &= \left( \frac{s^2}{{\omega\_0}^2} + 2\xi \frac{s}{{\omega\_0}^2} + 1 - \frac{\Omega^2}{{\omega\_0}^2} \right)^2 + \left( 2 \frac{\Omega}{\omega\_0}\frac{s}{\omega\_0} \right)^2 \\\\
+  G\_{kz} &= \left( \frac{s^2}{{\omega\_0}^2} - \frac{\Omega^2}{{\omega\_0}^2} + \alpha \right) \left( \frac{s^2}{{\omega\_0}^2} + 2\xi \frac{s}{{\omega\_0}^2} + 1 - \frac{\Omega^2}{{\omega\_0}^2} \right) + \left( 2 \frac{\Omega}{\omega\_0}\frac{s}{\omega\_0} \right)^2 \\\\
+  G\_{kc} &= \left( 2 \xi \frac{s}{\omega\_0} + 1 - \alpha \right) \left( 2 \frac{\Omega}{\omega\_0}\frac{s}{\omega\_0} \right)
+\end{align}
+
+</div>
+
+If we compare \\(G\_{kz}\\) and \\(G\_{fz}\\), we see that the spring in parallel adds a term \\(\alpha\\).
+In order to have two complex conjugate zeros (instead of real zeros):
+
+\begin{equation}
+  \alpha > \frac{\Omega^2}{{\omega\_0}^2} \quad \Leftrightarrow \quad k\_p > m \Omega^2
+\end{equation}
+
+
+### Plant Parameters {#plant-parameters}
+
+Let's define initial values for the model.
+
+```matlab
+  k = 1;    % Actuator Stiffness [N/m]
+  c = 0.05; % Actuator Damping [N/(m/s)]
+  m = 1;    % Payload mass [kg]
+```
+
+```matlab
+  xi = c/(2*sqrt(k*m));
+  w0 = sqrt(k/m); % [rad/s]
+```
+
+
+### Comparison of the Analytical Model and the Simscape Model {#comparison-of-the-analytical-model-and-the-simscape-model}
+
+The same transfer function from \\([F\_u, F\_v]\\) to \\([f\_u, f\_v]\\) is written down from the analytical model.
+
+```matlab
+  W = 0.1*w0; % [rad/s]
+
+  kp = 1.5*m*W^2;
+  cp = 0;
+```
+
+```matlab
+  open('rotating_frame.slx');
+```
+
+```matlab
+  %% Name of the Simulink File
+  mdl = 'rotating_frame';
+
+  %% Input/Output definition
+  clear io; io_i = 1;
+  io(io_i) = linio([mdl, '/K'], 1, 'openinput');  io_i = io_i + 1;
+  io(io_i) = linio([mdl, '/G'], 1, 'openoutput'); io_i = io_i + 1;
+
+  Giff = linearize(mdl, io, 0);
+
+  %% Input/Output definition
+  Giff.InputName  = {'Fu', 'Fv'};
+  Giff.OutputName = {'fu', 'fv'};
+```
+
+```matlab
+  w0p = sqrt((k + kp)/m);
+  xip = c/(2*sqrt((k+kp)*m));
+
+  Giff_th = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ...
+                     (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p));
+                     (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2 ];
+  Giff_th.InputName  = {'Fu', 'Fv'};
+  Giff_th.OutputName = {'fu', 'fv'};
+```
+
+<a id="figure--fig:plant-iff-kp-comp-simscape-analytical"></a>
+
+{{< figure src="figs/plant_iff_kp_comp_simscape_analytical.png" caption="<span class='figure-number'>Figure 18: </span>Comparison of the transfer functions from \\([F\_u, F\_v]\\) to \\([f\_u, f\_v]\\) between the Simscape model and the analytical one" >}}
+
+
+### Effect of the parallel stiffness on the IFF plant {#effect-of-the-parallel-stiffness-on-the-iff-plant}
+
+The rotation speed is set to \\(\Omega = 0.1 \omega\_0\\).
+
+```matlab
+  W = 0.1*w0; % [rad/s]
+```
+
+And the IFF plant (transfer function from \\([F\_u, F\_v]\\) to \\([f\_u, f\_v]\\)) is identified in three different cases:
+
+-   without parallel stiffness
+-   with a small parallel stiffness \\(k\_p < m \Omega^2\\)
+-   with a large parallel stiffness \\(k\_p > m \Omega^2\\)
+
+The results are shown in Figure [Figure 19](#figure--fig:plant-iff-kp).
+
+One can see that for \\(k\_p > m \Omega^2\\), the systems shows alternating complex conjugate poles and zeros.
+
+```matlab
+  kp = 0;
+
+  w0p = sqrt((k + kp)/m);
+  xip = c/(2*sqrt((k+kp)*m));
+
+  Giff = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ...
+      (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p));
+      (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2];
+```
+
+```matlab
+  kp = 0.5*m*W^2;
+  k = 1 - kp;
+
+  w0p = sqrt((k + kp)/m);
+  xip = c/(2*sqrt((k+kp)*m));
+
+  Giff_s = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ...
+      (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p));
+      (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2];
+```
+
+```matlab
+  kp = 1.5*m*W^2;
+  k = 1 - kp;
+
+  w0p = sqrt((k + kp)/m);
+  xip = c/(2*sqrt((k+kp)*m));
+
+  Giff_l = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ...
+      (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p));
+      (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2];
+```
+
+<a id="figure--fig:plant-iff-kp"></a>
+
+{{< figure src="figs/plant_iff_kp.png" caption="<span class='figure-number'>Figure 19: </span>Transfer function from \\([F\_u, F\_v]\\) to \\([f\_u, f\_v]\\) for \\(k\_p = 0\\), \\(k\_p < m \Omega^2\\) and \\(k\_p > m \Omega^2\\)" >}}
+
+
+### IFF when adding a spring in parallel {#iff-when-adding-a-spring-in-parallel}
+
+In Figure [Figure 20](#figure--fig:root-locus-iff-kp) is displayed the Root Locus in the three considered cases with
+
+\begin{equation}
+  K\_{\text{IFF}} = \frac{g}{s} \begin{bmatrix}
+  1 & 0 \\\\
+  0 & 1
+\end{bmatrix}
+\end{equation}
+
+One can see that for \\(k\_p > m \Omega^2\\), the root locus stays in the left half of the complex plane and thus the control system is unconditionally stable.
+
+Thus, decentralized IFF controller with pure integrators can be used if:
+
+\begin{equation}
+  k\_{p} > m \Omega^2
+\end{equation}
+
+<a id="figure--fig:root-locus-iff-kp"></a>
+
+{{< figure src="figs/root_locus_iff_kp.png" caption="<span class='figure-number'>Figure 20: </span>Root Locus" >}}
+
+<a id="figure--fig:root-locus-iff-kp-zoom"></a>
+
+{{< figure src="figs/root_locus_iff_kp_zoom.png" caption="<span class='figure-number'>Figure 21: </span>Root Locus" >}}
+
+
+### Effect of \\(k\_p\\) on the attainable damping {#effect-of-k-p-on-the-attainable-damping}
+
+However, having large values of \\(k\_p\\) may decrease the attainable damping.
+
+To study the second point, Root Locus plots for the following values of \\(k\_p\\) are shown in Figure [Figure 22](#figure--fig:root-locus-iff-kps).
+
+```matlab
+  kps = [2, 20, 40]*m*W^2;
+```
+
+It is shown that large values of \\(k\_p\\) decreases the attainable damping.
+
+<a id="figure--fig:root-locus-iff-kps"></a>
+
+{{< figure src="figs/root_locus_iff_kps.png" caption="<span class='figure-number'>Figure 22: </span>Root Locus plot" >}}
+
+```matlab
+  alphas = logspace(-2, 0, 100);
+
+  opt_xi = zeros(1, length(alphas)); % Optimal simultaneous damping
+  opt_gain = zeros(1, length(alphas)); % Corresponding optimal gain
+
+  Kiff = 1/s*eye(2);
+
+  for alpha_i = 1:length(alphas)
+      kp = alphas(alpha_i);
+      k = 1 - alphas(alpha_i);
+
+      w0p = sqrt((k + kp)/m);
+      xip = c/(2*sqrt((k+kp)*m));
+
+      Giff = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ...
+          (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p));
+          (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2];
+
+      fun = @(g)computeSimultaneousDamping(g, Giff, Kiff);
+
+      [g_opt, xi_opt] = fminsearch(fun, 2);
+      opt_xi(alpha_i) = 1/xi_opt;
+      opt_gain(alpha_i) = g_opt;
+  end
+```
+
+<a id="figure--fig:opt-damp-alpha"></a>
+
+{{< figure src="figs/opt_damp_alpha.png" caption="<span class='figure-number'>Figure 23: </span>Attainable damping ratio and corresponding controller gain for different parameter \\(\alpha\\)" >}}
+
+
+## Comparison {#comparison}
+
+<span class="org-target" id="org-target--sec-comparison"></span>
+
+Two modifications to adapt the IFF control strategy to rotating platforms have been proposed.
+These two methods are now compared in terms of added damping, closed-loop compliance and transmissibility.
+
+
+### Plant Parameters {#plant-parameters}
+
+Let's define initial values for the model.
+
+```matlab
+  k = 1;    % Actuator Stiffness [N/m]
+  c = 0.05; % Actuator Damping [N/(m/s)]
+  m = 1;    % Payload mass [kg]
+```
+
+```matlab
+  xi = c/(2*sqrt(k*m));
+  w0 = sqrt(k/m); % [rad/s]
+```
+
+The rotating speed is set to \\(\Omega = 0.1 \omega\_0\\).
+
+```matlab
+  W = 0.1*w0;
+```
+
+
+### Root Locus {#root-locus}
+
+IFF with High Pass Filter
+
+```matlab
+  wi = 0.1*w0; % [rad/s]
+
+  Giff = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ...
+          [(s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)) + (2*W*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*W*s/(w0^2) ; ...
+           (2*xi*s/w0 + 1)*2*W*s/(w0^2), (s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))+ (2*W*s/(w0^2))^2];
+```
+
+IFF With parallel Stiffness
+
+```matlab
+  kp = 5*m*W^2;
+  k = k - kp;
+
+  w0p = sqrt((k + kp)/m);
+  xip = c/(2*sqrt((k+kp)*m));
+
+  Giff_kp = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ...
+                     (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p));
+                     (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2 ];
+
+  k = k + kp;
+```
+
+<a id="figure--fig:comp-root-locus"></a>
+
+{{< figure src="figs/comp_root_locus.png" caption="<span class='figure-number'>Figure 24: </span>Root Locus plot - Comparison of IFF with additional high pass filter, IFF with additional parallel stiffness" >}}
+
+
+### Controllers - Optimal Gains {#controllers-optimal-gains}
+
+In order to compare to three considered Active Damping techniques, gains that yield maximum damping of all the modes are computed for each case.
+
+The obtained damping ratio and control are shown below.
+
+|                     | Obtained \\(\xi\\) | Control Gain |
+|---------------------|--------------------|--------------|
+| Modified IFF        | 0.83               | 1.99         |
+| IFF with \\(k\_p\\) | 0.83               | 2.02         |
+
+
+### Passive Damping - Critical Damping {#passive-damping-critical-damping}
+
+\begin{equation}
+  \xi = \frac{c}{2 \sqrt{km}}
+\end{equation}
+
+Critical Damping corresponds to to \\(\xi = 1\\), and thus:
+
+\begin{equation}
+  c\_{\text{crit}} = 2 \sqrt{km}
+\end{equation}
+
+```matlab
+  c_opt = 2*sqrt(k*m);
+```
+
+
+### Transmissibility And Compliance {#transmissibility-and-compliance}
+
+<span class="org-target" id="org-target--sec-comp-transmissibilty"></span>
+
+```matlab
+  open('rotating_frame.slx');
+```
+
+```matlab
+  %% Name of the Simulink File
+  mdl = 'rotating_frame';
+
+  %% Input/Output definition
+  clear io; io_i = 1;
+  io(io_i) = linio([mdl, '/dw'], 1, 'input');  io_i = io_i + 1;
+  io(io_i) = linio([mdl, '/fd'], 1, 'input');  io_i = io_i + 1;
+  io(io_i) = linio([mdl, '/Meas'], 1, 'output');  io_i = io_i + 1;
+```
+
+```matlab
+  G_ol = linearize(mdl, io, 0);
+
+  %% Input/Output definition
+  G_ol.InputName  = {'Dwx', 'Dwy', 'Fdx', 'Fdy'};
+  G_ol.OutputName = {'Dx', 'Dy'};
+```
+
+
+#### Passive Damping {#passive-damping}
+
+```matlab
+  kp = 0;
+  cp = 0;
+```
+
+```matlab
+  c_old = c;
+  c = c_opt;
+```
+
+```matlab
+  G_pas = linearize(mdl, io, 0);
+
+  %% Input/Output definition
+  G_pas.InputName  = {'Dwx', 'Dwy', 'Fdx', 'Fdy'};
+  G_pas.OutputName = {'Dx', 'Dy'};
+```
+
+```matlab
+  c = c_old;
+```
+
+```matlab
+  Kiff = opt_gain_iff/(wi + s)*tf(eye(2));
+```
+
+```matlab
+  G_iff = linearize(mdl, io, 0);
+
+  %% Input/Output definition
+  G_iff.InputName  = {'Dwx', 'Dwy', 'Fdx', 'Fdy'};
+  G_iff.OutputName = {'Dx', 'Dy'};
+```
+
+```matlab
+  kp = 5*m*W^2;
+  cp = 0.01;
+```
+
+```matlab
+  Kiff = opt_gain_kp/s*tf(eye(2));
+```
+
+```matlab
+  G_kp = linearize(mdl, io, 0);
+
+  %% Input/Output definition
+  G_kp.InputName  = {'Dwx', 'Dwy', 'Fdx', 'Fdy'};
+  G_kp.OutputName = {'Dx', 'Dy'};
+```
+
+<a id="figure--fig:comp-transmissibility"></a>
+
+{{< figure src="figs/comp_transmissibility.png" caption="<span class='figure-number'>Figure 25: </span>Comparison of the transmissibility" >}}
+
+<a id="figure--fig:comp-compliance"></a>
+
+{{< figure src="figs/comp_compliance.png" caption="<span class='figure-number'>Figure 26: </span>Comparison of the obtained Compliance" >}}
+
+
+## Notations {#notations}
+
+<span class="org-target" id="org-target--sec-notations"></span>
+
+|                                       | Mathematical Notation            | Matlab        | Unit    |
+|---------------------------------------|----------------------------------|---------------|---------|
+| Actuator Stiffness                    | \\(k\\)                          | `k`           | N/m     |
+| Actuator Damping                      | \\(c\\)                          | `c`           | N/(m/s) |
+| Payload Mass                          | \\(m\\)                          | `m`           | kg      |
+| Damping Ratio                         | \\(\xi = \frac{c}{2\sqrt{km}}\\) | `xi`          |         |
+| Actuator Force                        | \\(\bm{F}, F\_u, F\_v\\)         | `F` `Fu` `Fv` | N       |
+| Force Sensor signal                   | \\(\bm{f}, f\_u, f\_v\\)         | `f` `fu` `fv` | N       |
+| Relative Displacement                 | \\(\bm{d}, d\_u, d\_v\\)         | `d` `du` `dv` | m       |
+| Resonance freq. when \\(\Omega = 0\\) | \\(\omega\_0\\)                  | `w0`          | rad/s   |
+| Rotation Speed                        | \\(\Omega = \dot{\theta}\\)      | `W`           | rad/s   |
+| Low Pass Filter corner frequency      | \\(\omega\_i\\)                  | `wi`          | rad/s   |
+
+|                  | Mathematical Notation | Matlab | Unit    |
+|------------------|-----------------------|--------|---------|
+| Laplace variable | \\(s\\)               | `s`    |         |
+| Complex number   | \\(j\\)               | `j`    |         |
+| Frequency        | \\(\omega\\)          | `w`    | [rad/s] |
+
+<style>.csl-entry{text-indent: -1.5em; margin-left: 1.5em;}</style><div class="csl-bib-body">
+  <div class="csl-entry"><a id="citeproc_bib_item_1"></a>Dehaeze, T., and C. Collette. 2020. “Active Damping of Rotating Platforms Using Integral Force Feedback.” In <i>Proceedings of the International Conference on Modal Analysis Noise and Vibration Engineering (ISMA)</i>.</div>
+  <div class="csl-entry"><a id="citeproc_bib_item_2"></a>Dehaeze, Thomas. 2020. “Active Damping of Rotating Positioning Platforms.” Source Code on Zonodo. doi:<a href="https://doi.org/10.5281/zenodo.3894342">10.5281/zenodo.3894342</a>.</div>
+  <div class="csl-entry"><a id="citeproc_bib_item_3"></a>Dehaeze, Thomas, and Christophe Collette. 2021. “Active Damping of Rotating Platforms Using Integral Force Feedback.” <i>Engineering Research Express</i>. <a href="http://iopscience.iop.org/article/10.1088/2631-8695/abe803">http://iopscience.iop.org/article/10.1088/2631-8695/abe803</a>.</div>
+</div>