diff --git a/matlab/index.html b/matlab/index.html index 6ea41fa..c6922ad 100644 --- a/matlab/index.html +++ b/matlab/index.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Matlab Computation @@ -34,108 +34,1446 @@

Table of Contents

+ + +
+

1 System Description and Analysis

+

-Notations: + +

+
+ +
+

1.1 System description

+
+

+The system consists of one 2 degree of freedom translation stage on top of a spindle (figure 1). +

+ + +
+

rotating_xy_platform.png +

+

Figure 1: Figure caption

+
+ +

+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\). +

+ +

+The measurement is either the \(x-y\) displacement of the object located on top of the translation stage or the \(u-v\) displacement of the sample with respect to a fixed reference frame. +

+
+
+ +
+

1.2 Equations

+
+

+Based on the Figure 1, the equations of motions are: +

+
+\begin{equation} +\begin{bmatrix} d_u \\ d_v \end{bmatrix} = +\frac{\frac{1}{k}}{\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} +\begin{bmatrix} + \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} & 2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0} \\ + -2 \frac{\Omega}{\omega_0}\frac{s}{\omega_0} & \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} \\ +\end{bmatrix} +\begin{bmatrix} F_u \\ F_v \end{bmatrix} +\end{equation} + +
+ +

+Explain Coriolis and Centrifugal Forces (negative Stiffness) +

+
+
+ +
+

1.3 Numerical Values

+
+

+Let’s define initial values for the model. +

+
+
k = 1;    % Actuator Stiffness [N/m]
+c = 0.05; % Actuator Damping [N/(m/s)]
+m = 1;    % Payload mass [kg]
+
+
+ +
+
xi = c/(2*sqrt(k*m));
+w0 = sqrt(k/m); % [rad/s]
+
+
+
+
+ +
+

1.4 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 Figure 2, and one can see that the system becomes unstable for \(\Omega > \omega_0\) (the real part of one of the poles becomes positive). +

+ + +
+

campbell_diagram.png +

+

Figure 2: Campbell Diagram

+
+
+
+ +
+

1.5 Simscape Model

+
+

+Define the rotating speed for the Simscape Model. +

+
+
W = 0.1; % Rotation Speed [rad/s]
+
+
+ +

+The transfer function from \([F_u, F_v]\) to \([d_u, d_v]\) is identified from the Simscape model. +

+ +
+
%% 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'], 3, 'openoutput'); io_i = io_i + 1;
+
+
+ +
+
G = linearize(mdl, io, 0);
+
+%% Input/Output definition
+G.InputName  = {'Fu', 'Fv'};
+G.OutputName = {'du', 'dv'};
+
+
+
+
+ +
+

1.6 Comparison of the Analytical Model and the Simscape Model

+
+

+The same transfer function from \([F_u, F_v]\) to \([d_u, d_v]\) is written down from the analytical model. +

+
+
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 3 and are found to perfectly match. +

+ + +
+

plant_simscape_analytical.png +

+

Figure 3: 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

+
+
+
+ +
+

1.7 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. +

+
+
Ws = [0, 0.1, 0.5, 0.8, 1.1]*w0; % [rad/s]
+
+
+ +
+
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 Figure 4. +

+ + +
+

plant_compare_rotating_speed.png +

+

Figure 4: Comparison of the transfer functions from \([F_u, F_v]\) to \([d_u, d_v]\) for several rotating speed

+
+
+
+
+ +
+

2 Problem with pure Integral Force Feedback

+
+

+

+
+ +
+

2.1 Plant Parameters

+
+

+Let’s define initial values for the model. +

+
+
k = 1;    % Actuator Stiffness [N/m]
+c = 0.05; % Actuator Damping [N/(m/s)]
+m = 1;    % Payload mass [kg]
+
+
+ +
+
xi = c/(2*sqrt(k*m));
+w0 = sqrt(k/m); % [rad/s]
+
+
+
+
+ +
+

2.2 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: +

+
+\begin{equation} +\begin{bmatrix} f_{u} \\ f_{v} \end{bmatrix} = +\frac{1}{\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} +\begin{bmatrix} + (\frac{s^2}{{\omega_0}^2} - \frac{\Omega^2}{{\omega_0}^2}) (\frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2}) + (2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0})^2 & - (2 \xi \frac{s}{\omega_0} + 1) 2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0} \\ + (2 \xi \frac{s}{\omega_0} + 1) 2 \frac{\Omega}{\omega_0}\frac{s}{\omega_0} & (\frac{s^2}{{\omega_0}^2} - \frac{\Omega^2}{{\omega_0}^2}) (\frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2}) + (2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0})^2 \\ +\end{bmatrix} +\begin{bmatrix} F_u \\ F_v \end{bmatrix} +\end{equation} + +
+
+
+ +
+

2.3 Simscape Model

+
+

+The rotation speed is set to \(\Omega = 0.1 \omega_0\). +

+
+
W = 0.1*w0; % [rad/s]
+
+
+ +

+And the transfer function from \([F_u, F_v]\) to \([f_u, f_v]\) is identified using the Simscape model. +

+
+
%% 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;
+
+
+ +
+
Giff = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Giff.InputName  = {'Fu', 'Fv'};
+Giff.OutputName = {'fu', 'fv'};
+
+
+
+
+ +
+

2.4 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. +

+
+
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 5 and found to perfectly match. +

+ + +
+

plant_iff_comp_simscape_analytical.png +

+

Figure 5: Comparison of the transfer functions from \([F_u, F_v]\) to \([f_u, f_v]\) between the Simscape model and the analytical one

+
+
+
+ +
+

2.5 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. +

+
+
Ws = [0, 0.1, 0.5, 0.8, 1.1]*w0; % Rotating Speeds [rad/s]
+
+
+ +
+
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 6. +

+ +
+

plant_iff_compare_rotating_speed.png +

+

Figure 6: Comparison of the transfer functions from \([F_u, F_v]\) to \([f_u, f_v]\) for several rotating speed

+
+
+
+ +
+

2.6 Decentralized Integral Force Feedback

+
+

+Let’s take \(\Omega = \frac{\omega_0}{10}\). +

+
+
W = w0/10;
+
+
+ +
+
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 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} + +
+
g = 2;
+
+Kiff = g/s*tf(eye(2));
+
+
+ +

+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 7. +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. +

+ + +
+

root_locus_pure_iff.png +

+

Figure 7: Root Locus for the Decentralized Integral Force Feedback controller. Several rotating speed are shown.

+
+
+
+
+ +
+

3 Modified IFF (pseudo integrator)

+
+

+ +

+ +
+ +
+

3.1 Plant Parameters

+
+

+Let’s define initial values for the model. +

+
+
k = 1;    % Actuator Stiffness [N/m]
+c = 0.05; % Actuator Damping [N/(m/s)]
+m = 1;    % Payload mass [kg]
+
+
+ +
+
xi = c/(2*sqrt(k*m));
+w0 = sqrt(k/m); % [rad/s]
+
+
+
+
+ +
+

3.2 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: +

+
+
g = 2;
+wi = 0.1*w0;
+
+
+ +

+And the following rotating speed. +

+
+
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 8. +

+ +
+

loop_gain_modified_iff.png +

+

Figure 8: Loop Gain for the modified IFF controller

+
+
+
+ +
+

3.3 Root Locus

+
+

+As shown in the Root Locus plot (Figure 9), for some value of the gain, the system remains stable. +

+ + +
+

root_locus_modified_iff.png +

+

Figure 9: Root Locus for the modified IFF controller

+
+
+
+ +
+

3.4 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 10 for the following \(\omega_i\): +

+
+
wis = [0.01, 0.1, 0.5, 1]*w0; % [rad/s]
+
+
+ + +
+

root_locus_wi_modified_iff.png +

+

Figure 10: 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} +\label{org69ac90f} + g_\text{max} = \omega_i \left( \frac{{\omega_0}^2}{\Omega_2} - 1 \right) +\end{equation} + +

+While it seems that small \(\omega_i\) do allow more damping to be added to the system (Figure 10), the control gains may be limited to small values due to \eqref{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 11). +

+
+
wis = logspace(-2, 1, 31)*w0; % [rad/s]
+
+opt_zeta = 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);
+    gains = linspace(0, (w0^2/W^2 - 1)*wi, 100);
+
+    for g = gains
+        Kiff = (g/(wi+s))*eye(2);
+
+        [w, zeta] = damp(minreal(feedback(Giff, Kiff)));
+
+        if min(zeta) > opt_zeta(wi_i) && all(zeta > 0)
+            opt_zeta(wi_i) = min(zeta);
+            opt_gain(wi_i) = g;
+        end
+    end
+end
+
+
+ + +
+

mod_iff_damping_wi.png +

+

Figure 11: Simultaneous attainable damping of the closed loop poles as a function of \(\omega_i\)

+
+
+
+
+ + +
+

4 IFF with a stiffness in parallel with the force sensor

+
+

+ +

+
+ +
+

4.1 Plant Parameters

+
+

+Let’s define initial values for the model. +

+
+
k = 1;    % Actuator Stiffness [N/m]
+c = 0.05; % Actuator Damping [N/(m/s)]
+m = 1;    % Payload mass [kg]
+
+
+ +
+
xi = c/(2*sqrt(k*m));
+w0 = sqrt(k/m); % [rad/s]
+
+
+
+
+ +
+

4.2 Schematic

+
+ +
+

rotating_xy_platform_springs.png +

+

Figure 12: Figure caption

+
+
+
+ +
+

4.3 Physical Explanation

+
+
    +
  • Negative stiffness induced by gyroscopic effects
    • +
  • Zeros of the open-loop <=> Poles of the subsystem with the force sensors removes
    • +
  • As the zeros are the poles of the closed loop system for high gains, we want them to be in the left-half plane
    • +
  • Thus we want the zeros to be in the left half plant and thus the system with the force sensors stable
    • +
  • This can be done by adding springs in parallel with the force sensors with a stiffness larger than the virtual negative stiffness added by the gyroscopic effects
    • +
+ +

+The negative stiffness induced by the rotation is: +

+\begin{equation} + k_{n} = - m \Omega^2 +\end{equation} + +

+And thus, the stiffness in parallel should be such that: +

+\begin{equation} + k_{p} > m \Omega^2 +\end{equation} +
+
+ +
+

4.4 Equations

+
+

+The equations should be the same as before by taking into account the additional stiffness. +It then may be better to write it in terms of \(k\), \(c\), \(m\) instead of \(\omega_0\) and \(\xi\). +

+ +

+I just have to determine the measured force by the sensor +

+
+
+ +
+

4.5 Effect of the parallel stiffness on the IFF plant

+
+

+The rotation speed is set to \(\Omega = 0.1 \omega_0\). +

+
+
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 13. +

+ +

+One can see that for \(k_p > m \Omega^2\), the systems shows alternating complex conjugate poles and zeros. +

+ +
+
kp = 0;
+cp = 0;
+
+Giff = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Giff.InputName  = {'Fu', 'Fv'};
+Giff.OutputName = {'fu', 'fv'};
+
+
+ +
+
kp = 0.5*m*W^2;
+cp = 0.001;
+
+Giff_s = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Giff_s.InputName  = {'Fu', 'Fv'};
+Giff_s.OutputName = {'fu', 'fv'};
+
+
+ +
+
kp = 1.5*m*W^2;
+cp = 0.001;
+
+Giff_l = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Giff_l.InputName  = {'Fu', 'Fv'};
+Giff_l.OutputName = {'fu', 'fv'};
+
+
+ + +
+

plant_iff_kp.png +

+

Figure 13: 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\)

+
+
+
+
+

4.6 IFF when adding a spring in parallel

+
+

+In Figure 14 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} + + +
+

root_locus_iff_kp.png +

+

Figure 14: Root Locus

+
+
+
+ +
+

4.7 Effect of \(k_p\) on the attainable damping

+
+

+However, having large values of \(k_p\) may: +

+
    +
  • decrease the actuator force authority
    • +
  • decrease the attainable damping
    • +
+ +

+To study the second point, Root Locus plots for the following values of \(k_p\) are shown in Figure 15. +

+
+
kps = [1, 5, 10, 50]*m*W^2;
+cp = 0.01;
+
+
+ +

+It is shown that large values of \(k_p\) decreases the attainable damping. +

+ + +
+

root_locus_iff_kps.png +

+
+
+
+ +
+

4.8 Optimal Gain

+
+

+Let’s take \(k_p = 5 m \Omega^2\) and find the optimal IFF control gain \(g\) such that maximum damping are added to the poles of the closed loop system. +

+ +
+
kp = 5*m*W^2;
+cp = 0.01;
+
+Giff = linearize(mdl, io, 0);
+
+
+ +
+
opt_zeta = 0;
+opt_gain = 0;
+
+gains = logspace(-2, 4, 100);
+
+for g = gains
+    Kiff = (g/s)*eye(2);
+
+    [w, zeta] = damp(minreal(feedback(Giff, Kiff)));
+
+    if min(zeta) > opt_zeta && all(zeta > 0)
+        opt_zeta = min(zeta);
+        opt_gain = min(g);
+    end
+end
+
+
+ + +
+

root_locus_opt_gain_iff_kp.png +

+
+
+
+
+ +
+

5 Direct Velocity Feedback

+
+

+ +

+
+ +
+

5.1 Equations

+
+

+The sensed relative velocity are equal to: +

+
+\begin{equation} +\begin{bmatrix} \dot{d}_u \\ \dot{d}_v \end{bmatrix} = +\frac{s \frac{1}{k}}{\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} +\begin{bmatrix} + \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} & 2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0} \\ + -2 \frac{\Omega}{\omega_0}\frac{s}{\omega_0} & \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} \\ +\end{bmatrix} +\begin{bmatrix} F_u \\ F_v \end{bmatrix} +\end{equation} + +
+
+
+ +
+

5.2 Plant Parameters

+
+

+Let’s define initial values for the model. +

+
+
k = 1;    % Actuator Stiffness [N/m]
+c = 0.05; % Actuator Damping [N/(m/s)]
+m = 1;    % Payload mass [kg]
+
+
+ +
+
xi = c/(2*sqrt(k*m));
+w0 = sqrt(k/m); % [rad/s]
+
+
+
+
+ +
+

5.3 Plant - Bode Plot

+
+

+The rotating speed is set to \(\Omega = 0.1 \omega_0\). +

+
+
W = 0.1*w0;
+
+
+ +

+And the transfer function from \([F_u, F_v]\) to \([v_u, v_v]\) is identified using the Simscape model. +

+
+
%% 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;
+
+
+ +
+
Gdvf = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Gdvf.InputName  = {'Fu', 'Fv'};
+Gdvf.OutputName = {'Vu', 'Vv'};
+
+
+
+
+ +
+

5.4 Comparison of the Analytical Model and the Simscape Model

+
+

+The same transfer function from \([F_u, F_v]\) to \([v_u, v_v]\) is written down from the analytical model. +

+
+
Gdvf_th = (s/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)];
+
+Gdvf_th.InputName  = {'Fu', 'Fv'};
+Gdvf_th.OutputName = {'vu', 'vv'};
+
+
+ +

+The two are compared in Figure 5 and found to perfectly match. +

+ + +
+

plant_dvf_comp_simscape_analytical.png +

+

Figure 17: Comparison of the transfer functions from \([F_u, F_v]\) to \([v_u, v_v]\) between the Simscape model and the analytical one

+
+
+
+ +
+

5.5 Root Locus

+
+

+The Decentralized Direct Velocity Feedback controller consist of a pure gain on the diagonal: +

+\begin{equation} + K_{\text{DVF}}(s) = g \begin{bmatrix} + 1 & 0 \\ + 0 & 1 +\end{bmatrix} +\end{equation} + +

+The corresponding Root Locus plots for the following rotating speeds are shown in Figure 18. +

+
+
Ws = [0, 0.1, 0.5, 0.8, 1.1]*w0; % Rotating Speeds [rad/s]
+
+
+ +

+It is shown that for rotating speed \(\Omega < \omega_0\), the closed loop system is unconditionally stable and arbitrary damping can be added to the poles. +

+ + +
+

root_locus_dvf.png +

+

Figure 18: Root Locus for the Decentralized Direct Velocity Feedback controller. Several rotating speed are shown.

+
+
+
+
+ +
+

6 Comparison

+
+

+ +

+
+ +
+

6.1 Plant Parameters

+
+

+Let’s define initial values for the model. +

+
+
k = 1;    % Actuator Stiffness [N/m]
+c = 0.05; % Actuator Damping [N/(m/s)]
+m = 1;    % Payload mass [kg]
+
+
+ +
+
xi = c/(2*sqrt(k*m));
+w0 = sqrt(k/m); % [rad/s]
+
+
+ +

+The rotating speed is set to \(\Omega = 0.1 \omega_0\). +

+
+
W = 0.1*w0;
+
+
+
+
+ +
+

6.2 Root Locus

+
+
+
wi = 0.1*w0;
+
+
+ +
+
%% 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;
+
+
+ +
+
Giff = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Giff.InputName  = {'Fu', 'Fv'};
+Giff.OutputName = {'Fmu', 'Fmv'};
+
+
+
+
kp = 5*m*W^2;
+cp = 0.01;
+
+
+ +
+
%% 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;
+
+
+ +
+
Giff_kp = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Giff_kp.InputName  = {'Fu', 'Fv'};
+Giff_kp.OutputName = {'Fmu', 'Fmv'};
+
+
+
+
%% 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;
+
+
+ +
+
Gdvf = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Gdvf.InputName  = {'Fu', 'Fv'};
+Gdvf.OutputName = {'Vu', 'Vv'};
+
+
+ +
+

comp_root_locus.png +

+
+
+
+ +
+

6.3 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 IFF0.832.0
IFF with \(k_p\)0.842.01
DVF0.851.67
+
+
+ +
+

6.4 Transmissibility

+
+

+ +

+
+
%% 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, '/Meas'], 1, 'output');  io_i = io_i + 1;
+
+
+ +
+
Tol = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Tol.InputName  = {'Dwx', 'Dwy'};
+Tol.OutputName = {'Dx', 'Dy'};
+
+
+
+
Kiff = opt_gain_iff/(wi + s)*tf(eye(2));
+
+
+ +
+
%% 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, '/Meas'], 1, 'output');  io_i = io_i + 1;
+
+
+ +
+
Tiff = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Tiff.InputName  = {'Dwx', 'Dwy'};
+Tiff.OutputName = {'Dx', 'Dy'};
+
+
+
+
kp = 5*m*W^2;
+cp = 0.01;
+
+
+ +
+
Kiff = opt_gain_kp/s*tf(eye(2));
+
+
+ +
+
%% 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, '/Meas'], 1, 'output');  io_i = io_i + 1;
+
+
+ +
+
Tiff_kp = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Tiff_kp.InputName  = {'Dwx', 'Dwy'};
+Tiff_kp.OutputName = {'Dx', 'Dy'};
+
+
+
+
Kdvf = opt_gain_kp*tf(eye(2));
+
+
+ +
+
%% 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, '/Meas'], 1, 'output');  io_i = io_i + 1;
+
+
+ +
+
Tdvf = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Tdvf.InputName  = {'Dwx', 'Dwy'};
+Tdvf.OutputName = {'Dx', 'Dy'};
+
+
+ +
+

comp_transmissibility.png +

+
+
+
+ +
+

6.5 Compliance

+
+

+ +

+
+
%% Name of the Simulink File
+mdl = 'rotating_frame';
+
+%% Input/Output definition
+clear io; 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;
+
+
+ +
+
Col = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Col.InputName  = {'Fdx', 'Fdy'};
+Col.OutputName = {'Dx', 'Dy'};
+
+
+
+
Kiff = opt_gain_iff/(wi + s)*tf(eye(2));
+
+
+ +
+
Ciff = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Ciff.InputName  = {'Fdx', 'Fdy'};
+Ciff.OutputName = {'Dx', 'Dy'};
+
+
+
+
kp = 5*m*W^2;
+cp = 0.01;
+
+
+ +
+
Kiff = opt_gain_kp/s*tf(eye(2));
+
+
+ +
+
Ciff_kp = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Ciff_kp.InputName  = {'Fdx', 'Fdy'};
+Ciff_kp.OutputName = {'Dx', 'Dy'};
+
+
+
+
Kdvf = opt_gain_kp*tf(eye(2));
+
+
+ +
+
Cdvf = linearize(mdl, io, 0);
+
+%% Input/Output definition
+Cdvf.InputName  = {'Fdx', 'Fdy'};
+Cdvf.OutputName = {'Dx', 'Dy'};
+
+
+ +
+

comp_compliance.png +

+

Figure 21: Comparison of the obtained Compliance

+
+
+
+
+ +
+

7 Notations

+
+

+ +

@@ -171,1270 +1509,172 @@ Notations: + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
c N/(m/s)
Payload Mass\(m\)mkg
Damping Ratio\(\xi = \frac{c}{2\sqrt{km}}\)xi 
Actuator Force\(\bm{F}, F_u, F_v\)F Fu FvN
Force Sensor signal\(\bm{f}, f_u, f_v\)f fu fvN
Relative Displacement\(\bm{d}, d_u, d_v\)d du dvm
Relative Velocity\(\bm{v}, v_u, v_v\)v vu vvm/s
Resonance freq. when \(\Omega = 0\)\(\omega_0\)w0rad/s
Rotation Speed\(\Omega = \dot{\theta}\)Wrad/s
Low Pass Filter corner frequency\(\omega_i\)wirad/s
-
-

1 System Description and Analysis

-
-

- -

-
- -
-

1.1 System description

-
-

-The system consists of one 2 degree of freedom translation stage on top of a spindle (figure 1). -

- - -
-

rotating_xy_platform.png -

-

Figure 1: Figure caption

-
- -

-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\). -

- -

-The measurement is either the \(x-y\) displacement of the object located on top of the translation stage or the \(u-v\) displacement of the sample with respect to a fixed reference frame. -

-
-
- -
-

1.2 Equations

-
-

-Based on the Figure 1. -

- -\begin{equation} -\begin{bmatrix} d_u \\ d_v \end{bmatrix} = -\frac{\frac{1}{k}}{\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} -\begin{bmatrix} - \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} & 2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0} \\ - -2 \frac{\Omega}{\omega_0}\frac{s}{\omega_0} & \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} \\ -\end{bmatrix} -\begin{bmatrix} F_u \\ F_v \end{bmatrix} -\end{equation} -
-
- -
-

1.3 Numerical Values

-
-
-
k = 2;
-m = 1;
-c = 0.05;
-xi = c/(2*sqrt(k*m));
-
-w0 = sqrt(k/m);
-
-
-
-
- -
-

1.4 Campbell Diagram

-
-

-Compute the poles -

-
-
wrs = linspace(0, 2, 51); % [rad/s]
-
-polesvc = zeros(4, length(wrs));
-
-for i = 1:length(wrs)
-    wr = wrs(i);
-    polei = pole(1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wrs(i)^2)/(w0^2))^2 + (2*wrs(i)*s/(w0^2))^2));
-    [~, i_sort] = sort(imag(polei));
-    polesvc(:, i) = polei(i_sort);
-end
-
-
-
-
- -
-

1.5 Simscape Model

-
-
-
open('rotating_frame.slx');
-
-
- -
-
k = 2;
-m = 1;
-c = 0.05;
-xi = c/(2*sqrt(k*m));
-
-w0 = sqrt(k/m);
-wr = 0.1;
-
-
- -
-
%% Name of the Simulink File
-mdl = 'rotating_frame';
-
-%% Input/Output definition
-clear io; io_i = 1;
-io(io_i) = linio([mdl, '/Fu'], 1, 'openinput');  io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Fv'], 1, 'openinput');  io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Tuv'], 1, 'openoutput'); io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Tuv'], 2, 'openoutput'); io_i = io_i + 1;
-
-
- -
-
G = linearize(mdl, io, 0);
-
-%% Input/Output definition
-G.InputName  = {'Fu', 'Fv'};
-G.OutputName = {'du', 'dv'};
-
-
-
-
- -
-

1.6 Comparison with the model

-
-
-
Ga = (1/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2))^2 + (2*wr*s/(w0^2))^2) * ...
-     [(s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2), 2*wr*s/(w0^2) ; ...
-      -2*wr*s/(w0^2), (s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2)];
-
-
- -
-
figure; bode(G, 'k-', Ga, 'r--')
-
-
-
-
-
- -
-

2 Integral Force Feedback

-
-
-
-

2.1 IFF with pure integrator

-
-
-
-

2.1.1 Numerical Values

-
-
-
k = 1;
-m = 1;
-xi = 0.01;
-c = 2*xi*sqrt(k*m);
-w0 = sqrt(k/m);
-wr = 0.1*w0;
-
-
-
-
- -
-

2.1.2 Equations

-
-

-The sensed forces are equal to: -

-\begin{equation} -\begin{bmatrix} F_{um} \\ F_{vm} \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} - -
-\begin{equation} -\begin{bmatrix} F_{um} \\ F_{vm} \end{bmatrix} = -\frac{1}{\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} -\begin{bmatrix} - (\frac{s^2}{{\omega_0}^2} - \frac{\Omega^2}{{\omega_0}^2}) (\frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2}) + (2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0})^2 & - (2 \xi \frac{s}{\omega_0} + 1) 2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0} \\ - (2 \xi \frac{s}{\omega_0} + 1) 2 \frac{\Omega}{\omega_0}\frac{s}{\omega_0} & (\frac{s^2}{{\omega_0}^2} - \frac{\Omega^2}{{\omega_0}^2}) (\frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2}) + (2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0})^2 \\ -\end{bmatrix} -\begin{bmatrix} F_u \\ F_v \end{bmatrix} -\end{equation} - -
- -
-
Giffa = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2))^2 + (2*wr*s/(w0^2))^2) * ...
-       [(s^2/w0^2 - wr^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2)) + (2*wr*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*wr*s/(w0^2) ; ...
-        (2*xi*s/w0 + 1)*2*wr*s/(w0^2), (s^2/w0^2 - wr^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2))+ (2*wr*s/(w0^2))^2];
-
-
-
-
- -
-

2.1.3 Poles and Zeros

-
-
-
syms wr w0 xi positive
-assumealso(w0 > wr)
-syms x
-
-
- -
-
z = (x^2/w0^2 - wr^2/w0^2)*((x^2)/(w0^2) + 1 - (wr^2)/(w0^2)) + (2*wr*x/(w0^2))^2 == 0
-p = ((x^2)/(w0^2) + 1 - (wr^2)/(w0^2))^2 + (2*wr*x/(w0^2))^2 == 0
-
-
- -
-
solve(p, x)
-
-
- -
-
solve(z, x)
-
-
- -

-The zeros are the roots of: -

-\begin{equation} - \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 = 0 -\end{equation} - -

-Poles (without damping) -

-\begin{equation} - \left(\begin{array}{c} -w_{0}\,1{}\mathrm{i}-\mathrm{wr}\,1{}\mathrm{i}\\ -w_{0}\,1{}\mathrm{i}+\mathrm{wr}\,1{}\mathrm{i}\\ w_{0}\,1{}\mathrm{i}-\mathrm{wr}\,1{}\mathrm{i}\\ w_{0}\,1{}\mathrm{i}+\mathrm{wr}\,1{}\mathrm{i} \end{array}\right) -\end{equation} - -

-Zeros (without damping) -

-\begin{equation} - \left(\begin{array}{c} -\sqrt{-\frac{w_{0}\,\sqrt{{w_{0}}^2+8\,{\mathrm{wr}}^2}}{2}-\frac{{w_{0}}^2}{2}-{\mathrm{wr}}^2}\\ -\sqrt{\frac{w_{0}\,\sqrt{{w_{0}}^2+8\,{\mathrm{wr}}^2}}{2}-\frac{{w_{0}}^2}{2}-{\mathrm{wr}}^2}\\ \sqrt{-\frac{w_{0}\,\sqrt{{w_{0}}^2+8\,{\mathrm{wr}}^2}}{2}-\frac{{w_{0}}^2}{2}-{\mathrm{wr}}^2}\\ \sqrt{\frac{w_{0}\,\sqrt{{w_{0}}^2+8\,{\mathrm{wr}}^2}}{2}-\frac{{w_{0}}^2}{2}-{\mathrm{wr}}^2} \end{array}\right) -\end{equation} -
-
- -
-

2.1.4 Simscape Model

-
-
-
%% Name of the Simulink File
-mdl = 'rotating_frame';
-
-%% Input/Output definition
-clear io; io_i = 1;
-io(io_i) = linio([mdl, '/Fu'], 1, 'openinput');  io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Fv'], 1, 'openinput');  io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Tuv'], 5, 'openoutput'); io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Tuv'], 6, 'openoutput'); io_i = io_i + 1;
-
-
- -
-
Giff = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Giff.InputName  = {'Fu', 'Fv'};
-Giff.OutputName = {'Fmu', 'Fmv'};
-
-
- -
-
figure; bode(Giff, Giffa)
-
-
-
-
- -
-

2.1.5 IFF Plant

-
-
-

2.1.6 Loop Gain

-
-

-Let’s take \(\Omega = \frac{\omega_0}{10}\). -

-
-
ws = 0.1*w0;
-G_iff = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (ws^2)/(w0^2))^2 + (2*ws*s/(w0^2))^2) * ...
-        [(s^2/w0^2 - ws^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (ws^2)/(w0^2)) + (2*ws*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*ws*s/(w0^2) ; ...
-         (2*xi*s/w0 + 1)*2*ws*s/(w0^2), (s^2/w0^2 - ws^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (ws^2)/(w0^2))+ (2*ws*s/(w0^2))^2];
-
-
-
-
-
- -
-

2.1.7 Root Locus

-
-
-
-

2.2 Modified IFF (pseudo integrator)

-
-
-
-

2.2.1 Control Law

-
-

-Let’s take the integral feedback controller as a low pass filter (pseudo integrator): -

-\begin{equation} - K_{\text{IFF}} = \begin{bmatrix} - g\frac{\omega_i}{\omega_i + s} & 0 \\ - 0 & g\frac{\omega_i}{\omega_i + s} -\end{bmatrix} -\end{equation} - -
-
xi = 0.005;
-w0 = 1;
-ws = 0.1*w0;
-
-G_iff = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (ws^2)/(w0^2))^2 + (2*ws*s/(w0^2))^2) * ...
-        [(s^2/w0^2 - ws^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (ws^2)/(w0^2)) + (2*ws*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*ws*s/(w0^2) ; ...
-         (2*xi*s/w0 + 1)*2*ws*s/(w0^2), (s^2/w0^2 - ws^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (ws^2)/(w0^2))+ (2*ws*s/(w0^2))^2];
-
-
-
- -
-
g = 100;
-wi = ws;
-
-K_iff = (g/(1+s/wi))*eye(2);
-
-
-
-
- -
-

2.2.2 Loop Gain

-
-
-

2.2.3 Root Locus

-
-
-

2.2.4 Optimal Gain

-
-

-The DC gain for Giff is (for \(\Omega < \omega_0\)): -

-\begin{equation} - G_{\text{IFF}}(\omega = 0) = \frac{1}{1 - \frac{{\omega_0}^2}{\Omega^2}} \begin{bmatrix} - 1 & 0 \\ - 0 & 1 - \end{bmatrix} -\end{equation} - -

-The maximum gain where is system is still stable is -

-\begin{equation} - g_\text{max} = \frac{{\omega_0}^2}{\Omega^2} - 1 -\end{equation} - - -

-Let’s find the gain that maximize the simultaneous damping of the two modes. -

-
-
K_opt = (opt_gain/(1+s/wi))*eye(2);
-
-G_cl = feedback(G_iff, K_opt);
-
-
-
-
-
- -
-

2.3 Stiffness in parallel with the force sensor

-
-
-
-

2.3.1 Schematic

-
- -
-

2.3.2 Equations

-
-

-The equations should be the same as before by taking \(k = k^\prime + k_a\). -I just have to determine the measured force by the sensor -

-
-
- -
-

2.3.3 Parameters

-
-
-
k = 1;
-m = 1;
-c = 0.05;
-xi = c/(2*sqrt(k*m));
-
-w0 = sqrt(k/m);
-wr = 0.1;
-
-
-
-
- -
-

2.3.4 IFF Plant

-
-
-
open('rotating_frame.slx');
-
-
- -
-
%% Name of the Simulink File
-mdl = 'rotating_frame';
-
-%% Input/Output definition
-clear io; io_i = 1;
-io(io_i) = linio([mdl, '/Fu'], 1, 'openinput');  io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Fv'], 1, 'openinput');  io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Tuv'], 5, 'openoutput'); io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Tuv'], 6, 'openoutput'); io_i = io_i + 1;
-
-
- -
-
kp = 0;
-cp = 0;
-
-
- -
-
Giff = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Giff.InputName  = {'Fu', 'Fv'};
-Giff.OutputName = {'Fmu', 'Fmv'};
-
-
- -
-
kp = 0.5*m*wr^2;
-cp = 0.01;
-
-
- -
-
Giffa = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Giffa.InputName  = {'Fu', 'Fv'};
-Giffa.OutputName = {'Fmu', 'Fmv'};
-
-
- -
-
kp = 1.5*m*wr^2;
-cp = 0.01;
-
-
- -
-
Giffb = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Giffb.InputName  = {'Fu', 'Fv'};
-Giffb.OutputName = {'Fmu', 'Fmv'};
-
-
- -

-Comparison with the model -

-
-
figure; bode(Giff, 'k-', Giffa, 'b--', Giffb, 'r--')
-
-
-
-
- -
-

2.3.5 Parallel Stiffness effect

-
-

-Pure IFF controller can be used if: -

-\begin{equation} - k_{p} > m \Omega^2 -\end{equation} - -

-However, having large values of \(k_p\) may: -

-
    -
  • decrease the actuator stroke
    • -
  • decrease the attainable damping (section about optimal value)
    • -
-
-
- -
-

2.3.6 Root locus

-
-
-

2.3.7 Optimal value of \(k_p\)

-
-
-
kps = [0, 0.5, 1, 2, 10]*m*wr^2;
-cp = 0.01;
-
-
-
-
- -
-

2.3.8 Physical Explanation

-
-

-Negative stiffness -Zeros are for high loop gain: remove the force sensor -Thus the stiffness in parallel should be higher than the virtual negative stiffness added by the gyroscopic effects -

-
-
-
- -
-

2.4 Comparison

-
-

-Fix the parameters, and see how the root locus change and the final damped system -

-
-
k = 1;
-m = 1;
-c = 0.05;
-
-xi = c/(2*sqrt(k*m));
-
-w0 = sqrt(k/m);
-wr = 0.1;
-
-
- -
-
kp = 0;
-cp = 0;
-
-
- -
-
%% Name of the Simulink File
-mdl = 'rotating_frame';
-
-%% Input/Output definition
-clear io; io_i = 1;
-io(io_i) = linio([mdl, '/Fu'], 1, 'openinput');  io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Fv'], 1, 'openinput');  io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Tuv'], 5, 'openoutput'); io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Tuv'], 6, 'openoutput'); io_i = io_i + 1;
-
-
- -
-
Giff = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Giff.InputName  = {'Fu', 'Fv'};
-Giff.OutputName = {'Fmu', 'Fmv'};
-
-
- -

-Modified IFF -

-
-
wi = w0;
-
-
- -

-With parallel Stiffness -

-
-
kp = 2*m*wr^2;
-cp = 0.01;
-
-
- -
-
Giff_kp = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Giff_kp.InputName  = {'Fu', 'Fv'};
-Giff_kp.OutputName = {'Fmu', 'Fmv'};
-
-
- -

-With parallel stiffness => unconditional stability and more damping but have to add mechanical parts -

-
-
-
- -
-

3 Direct Velocity Feedback

-
-
-
-

3.1 Equations

-
-

-The sensed relative velocity are equal to: -

- -
-\begin{equation} -\begin{bmatrix} \dot{d}_u \\ \dot{d}_v \end{bmatrix} = -\frac{s \frac{1}{k}}{\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} -\begin{bmatrix} - \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} & 2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0} \\ - -2 \frac{\Omega}{\omega_0}\frac{s}{\omega_0} & \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} \\ -\end{bmatrix} -\begin{bmatrix} F_u \\ F_v \end{bmatrix} -\end{equation} - -
-
-
- -
-

3.2 Numerical Values

-
-
-
k = 1;
-m = 1;
-c = 0.05;
-xi = c/(2*sqrt(k*m));
-
-w0 = sqrt(k/m);
-
-
- -
-
Ga = (s/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2))^2 + (2*wr*s/(w0^2))^2) * ...
-     [(s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2), 2*wr*s/(w0^2) ; ...
-      -2*wr*s/(w0^2), (s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2)];
-
-
-
-
- -
-

3.3 Simscape Model

-
-
-
%% Name of the Simulink File
-mdl = 'rotating_frame';
-
-%% Input/Output definition
-clear io; io_i = 1;
-io(io_i) = linio([mdl, '/Fu'], 1, 'openinput');  io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Fv'], 1, 'openinput');  io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Tuv'], 3, 'openoutput'); io_i = io_i + 1;
-io(io_i) = linio([mdl, '/Tuv'], 4, 'openoutput'); io_i = io_i + 1;
-
-
- -
-
G = linearize(mdl, io, 0);
-
-%% Input/Output definition
-G.InputName  = {'Fu', 'Fv'};
-G.OutputName = {'Vu', 'Vv'};
-
-
-
-
- -
-

3.4 DVF Plant

-
-
-
figure; bode(G, Ga)
-
-
-
-
- -
-

3.5 Loop Gain

-
-

-Let’s take \(\Omega = \frac{\omega_0}{10}\). -

-
-
ws = 0.1*w0;
-G_dvf = (s/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2))^2 + (2*wr*s/(w0^2))^2) * ...
-     [(s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2), 2*wr*s/(w0^2) ; ...
-      -2*wr*s/(w0^2), (s^2)/(w0^2) + 2*xi*s/w0 + 1 - (wr^2)/(w0^2)];
-
-
-
-
-
- -
-

3.6 Root Locus

-
-
-
-

4 Comparison

-
-
-
-

4.1 Parameters

-
-
-
k = 1;
-m = 1;
-c = 0.05;
-
-xi = c/(2*sqrt(k*m));
-
-w0 = sqrt(k/m);
-wr = 0.1;
-
-
-
-
- -
-

4.2 Root Locus

-
-
-
-

4.2.1 Pseudo Integrator IFF

-
-
-
kp = 0;
-cp = 0;
-
-
- -
-
%% 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;
-
-
- -
-
Kiff = tf(zeros(2));
-Kdvf = tf(zeros(2));
-
-
- -
-
Giff = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Giff.InputName  = {'Fu', 'Fv'};
-Giff.OutputName = {'Fmu', 'Fmv'};
-
-
-
-
- -
-

4.2.2 IFF With parallel Stiffness

-
-
-
kp = 2*m*wr^2;
-cp = 0.01;
-
-
- -
-
%% 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;
-
-
- -
-
Giff_kp = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Giff_kp.InputName  = {'Fu', 'Fv'};
-Giff_kp.OutputName = {'Fmu', 'Fmv'};
-
-
-
-
- -
-

4.2.3 DVF

-
-
-
%% 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;
-
-
- -
-
Gdvf = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Gdvf.InputName  = {'Fu', 'Fv'};
-Gdvf.OutputName = {'Vu', 'Vv'};
-
-
-
-
- -
-

4.2.4 Root Locus

-
-
-
wi = w0;
-
-
-
-
-
- -
-

4.3 Controllers - Optimal Gains

-
-

-Estimate the controller gain that yields good damping in all cases. -

- -

-Pseudo integrator IFF: -Parallel Stiffness -DVF: -

-
-
opt_zeta_iff, opt_zeta_kp, opt_zeta_dvf
-
-
-
-
- -
-

4.4 Transmissibility

-
-
-
-

4.4.1 Open Loop

-
-
-
wi = w0;
-
-
- -
-
Kdvf = tf(zeros(2));
-Kiff = tf(zeros(2));
-
-
- -
-
kp = 0;
-cp = 0;
-
-
- -
-
%% 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, '/Meas'], 1, 'output');  io_i = io_i + 1;
-
-
- -
-
Tol = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Tol.InputName  = {'Dwx', 'Dwy'};
-Tol.OutputName = {'Dx', 'Dy'};
-
-
-
-
- -
-

4.4.2 Pseudo Integrator IFF

-
-
-
wi = w0;
-
-
- -
-
Kdvf = tf(zeros(2));
-Kiff = opt_gain_iff/(1 + s/wi)*tf(eye(2));
-
-
- -
-
kp = 0;
-cp = 0;
-
-
- -
-
%% 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, '/Meas'], 1, 'output');  io_i = io_i + 1;
-
-
- -
-
Tiff = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Tiff.InputName  = {'Dwx', 'Dwy'};
-Tiff.OutputName = {'Dx', 'Dy'};
-
-
-
-
- -
-

4.4.3 IFF With parallel Stiffness

-
-
-
kp = 2*m*wr^2;
-cp = 0.01;
-
-
- -
-
Kdvf = tf(zeros(2));
-Kiff = opt_gain_kp/s*tf(eye(2));
-
-
- -
-
%% 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, '/Meas'], 1, 'output');  io_i = io_i + 1;
-
-
- -
-
Tiff_kp = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Tiff_kp.InputName  = {'Dwx', 'Dwy'};
-Tiff_kp.OutputName = {'Dx', 'Dy'};
-
-
-
-
- -
-

4.4.4 DVF

-
-
-
kp = 0;
-cp = 0;
-
-
- -
-
Kdvf = opt_gain_kp*tf(eye(2));
-Kiff = tf(zeros(2));
-
-
- -
-
%% 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, '/Meas'], 1, 'output');  io_i = io_i + 1;
-
-
- -
-
Tdvf = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Tdvf.InputName  = {'Dwx', 'Dwy'};
-Tdvf.OutputName = {'Dx', 'Dy'};
-
-
-
-
- -
-

4.4.5 Transmissibility

-
-
-
freqs = logspace(-2, 1, 1000);
-
-figure;
-hold on;
-plot(freqs, abs(squeeze(freqresp(Tiff(1,1), freqs))), ...
-     'DisplayName', 'IFF Pseudo int')
-plot(freqs, abs(squeeze(freqresp(Tiff_kp(1,1), freqs))), ...
-     'DisplayName', 'IFF Paral. stiff')
-plot(freqs, abs(squeeze(freqresp(Tdvf(1,1), freqs))), ...
-     'DisplayName', 'DVF')
-plot(freqs, abs(squeeze(freqresp(Tol(1,1), freqs))), 'k-', ...
-     'DisplayName', 'IFF Pseudo int')
-hold off;
-set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-ylabel('Frequency [rad/s]'); ylabel('Amplitude [m/m]');
-legend('location', 'northwest');
-
-
-
-
-
- -
-

4.5 Compliance

-
-
-
-

4.5.1 Open Loop

-
-
-
Kdvf = tf(zeros(2));
-Kiff = tf(zeros(2));
-
-
- -
-
kp = 0;
-cp = 0;
-
-
- -
-
%% Name of the Simulink File
-mdl = 'rotating_frame';
-
-%% Input/Output definition
-clear io; 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;
-
-
- -
-
Col = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Col.InputName  = {'Fdx', 'Fdy'};
-Col.OutputName = {'Dx', 'Dy'};
-
-
-
-
- -
-

4.5.2 Pseudo Integrator IFF

-
-
-
wi = w0;
-
-
- -
-
Kdvf = tf(zeros(2));
-Kiff = opt_gain_iff/(1 + s/wi)*tf(eye(2));
-
-
- -
-
kp = 0;
-cp = 0;
-
-
- -
-
%% Name of the Simulink File
-mdl = 'rotating_frame';
-
-%% Input/Output definition
-clear io; 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;
-
-
- -
-
Ciff = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Ciff.InputName  = {'Fdx', 'Fdy'};
-Ciff.OutputName = {'Dx', 'Dy'};
-
-
-
-
- -
-

4.5.3 IFF With parallel Stiffness

-
-
-
kp = 2*m*wr^2;
-cp = 0.01;
-
-
- -
-
Kdvf = tf(zeros(2));
-Kiff = opt_gain_kp/s*tf(eye(2));
-
-
- -
-
%% Name of the Simulink File
-mdl = 'rotating_frame';
-
-%% Input/Output definition
-clear io; 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;
-
-
- -
-
Ciff_kp = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Ciff_kp.InputName  = {'Fdx', 'Fdy'};
-Ciff_kp.OutputName = {'Dx', 'Dy'};
-
-
-
-
- -
-

4.5.4 DVF

-
-
-
kp = 0;
-cp = 0;
-
-
- -
-
Kdvf = opt_gain_kp*tf(eye(2));
-Kiff = tf(zeros(2));
-
-
- -
-
%% Name of the Simulink File
-mdl = 'rotating_frame';
-
-%% Input/Output definition
-clear io; 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;
-
-
- -
-
Cdvf = linearize(mdl, io, 0);
-
-%% Input/Output definition
-Cdvf.InputName  = {'Fdx', 'Fdy'};
-Cdvf.OutputName = {'Dx', 'Dy'};
-
-
-
-
- -
-

4.5.5 Compliance

-
-
-
freqs = logspace(-2, 1, 1000);
-
-figure;
-hold on;
-plot(freqs, abs(squeeze(freqresp(Ciff(1,1), freqs))), ...
-     'DisplayName', 'IFF Pseudo int')
-plot(freqs, abs(squeeze(freqresp(Ciff_kp(1,1), freqs))), ...
-     'DisplayName', 'IFF Paral. stiff')
-plot(freqs, abs(squeeze(freqresp(Cdvf(1,1), freqs))), ...
-     'DisplayName', 'DVF')
-plot(freqs, abs(squeeze(freqresp(Col(1,1), freqs))), 'k-', ...
-     'DisplayName', 'IFF Pseudo int')
-hold off;
-set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-ylabel('Frequency [rad/s]'); ylabel('Compliance [m/N]');
-legend('location', 'northwest');
-% ylim([0, 1e3]);
-% xlim([freqs(1), freqs(end)]);
-
-
-
-
+ + + +++ ++ ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
 Mathematical NotationMatlabUnit
Laplace variable\(s\)s 
Complex number\(j\)j 
Frequency\(\omega\)w[rad/s]
+ + + + +++ ++ ++ ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
 Mathematical NotationMatlabUnit
IFF Plant\(\bm{G}_\text{IFF}(s) = \frac{\bm{f}}{\bm{F}}\)GiffN/N
DVF Plant\(\bm{G}_\text{DVF}(s) = \frac{\bm{v}}{\bm{F}}\)Gdvf(m/s)/N
IFF Controller\(\bm{K}_\text{IFF}(s)\)Kiff 
DVF Controller\(\bm{K}_\text{DVF}(s)\)Kdvf 

Author: Thomas Dehaeze

-

Created: 2020-06-10 mer. 17:21

+

Created: 2020-06-12 ven. 19:29

diff --git a/matlab/index.org b/matlab/index.org index a76316b..ba137be 100644 --- a/matlab/index.org +++ b/matlab/index.org @@ -101,7 +101,7 @@ The Campbell Diagram displays the evolution of the real and imaginary parts of t It is shown in Figure [[fig:campbell_diagram]], and one can see that the system becomes unstable for $\Omega > \omega_0$ (the real part of one of the poles becomes positive). -#+begin_src matlab :exports code +#+begin_src matlab :exports none Ws = linspace(0, 2, 51); % Vector of considered rotation speed [rad/s] p_ws = zeros(4, length(Ws)); @@ -117,7 +117,7 @@ It is shown in Figure [[fig:campbell_diagram]], and one can see that the system clear pole_G; #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none figure; ax1 = subplot(1,2,1); @@ -185,7 +185,7 @@ Define the rotating speed for the Simscape Model. W = 0.1; % Rotation Speed [rad/s] #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none Kiff = tf(zeros(2)); Kdvf = tf(zeros(2)); @@ -223,7 +223,7 @@ The same transfer function from $[F_u, F_v]$ to $[d_u, d_v]$ is written down fro Both transfer functions are compared in Figure [[fig:plant_simscape_analytical]] and are found to perfectly match. -#+begin_src matlab :exports code +#+begin_src matlab :exports none freqs = logspace(-1, 1, 1000); figure; @@ -302,7 +302,7 @@ The transfer functions from $[F_u, F_v]$ to $[d_u, d_v]$ are identified for the They are compared in Figure [[fig:plant_compare_rotating_speed]]. -#+begin_src matlab :exports code +#+begin_src matlab :exports none freqs = logspace(-2, 1, 1000); figure; @@ -402,7 +402,7 @@ Let's define initial values for the model. w0 = sqrt(k/m); % [rad/s] #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none kp = 0; % [N/m] cp = 0; % [N/(m/s)] #+end_src @@ -438,7 +438,7 @@ The rotation speed is set to $\Omega = 0.1 \omega_0$. W = 0.1*w0; % [rad/s] #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none Kiff = tf(zeros(2)); Kdvf = tf(zeros(2)); #+end_src @@ -472,7 +472,7 @@ The same transfer function from $[F_u, F_v]$ to $[f_u, f_v]$ is written down fro The two are compared in Figure [[fig:plant_iff_comp_simscape_analytical]] and found to perfectly match. -#+begin_src matlab :exports code +#+begin_src matlab :exports none freqs = logspace(-1, 1, 1000); figure; @@ -550,7 +550,7 @@ The transfer functions from $[F_u, F_v]$ to $[f_u, f_v]$ are identified for the #+end_src The obtained transfer functions are shown in Figure [[fig:plant_iff_compare_rotating_speed]]. -#+begin_src matlab :exports code +#+begin_src matlab :exports none freqs = logspace(-2, 1, 1000); figure; @@ -649,7 +649,7 @@ It is shown that for non-null rotating speed, one pole is bound to the right-hal xlim([freqs(1), freqs(end)]); #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none figure; gains = logspace(-2, 4, 100); @@ -722,7 +722,7 @@ Let's define initial values for the model. w0 = sqrt(k/m); % [rad/s] #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none kp = 0; % [N/m] cp = 0; % [N/(m/s)] #+end_src @@ -743,12 +743,12 @@ Let's arbitrary choose the following control parameters: wi = 0.1*w0; #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none Kiff = (g/(wi+s))*eye(2); #+end_src And the following rotating speed. -#+begin_src matlab :exports code +#+begin_src matlab :exports none W = 0.1*w0; #+end_src @@ -759,7 +759,7 @@ And the following rotating speed. #+end_src The obtained Loop Gain is shown in Figure [[fig:loop_gain_modified_iff]]. -#+begin_src matlab :exports code +#+begin_src matlab :exports none freqs = logspace(-2, 1, 1000); figure; @@ -801,7 +801,7 @@ The obtained Loop Gain is shown in Figure [[fig:loop_gain_modified_iff]]. ** Root Locus As shown in the Root Locus plot (Figure [[fig:root_locus_modified_iff]]), for some value of the gain, the system remains stable. -#+begin_src matlab :exports code +#+begin_src matlab :exports none figure; gains = logspace(-2, 4, 100); @@ -877,7 +877,7 @@ In order to visualize the effect of $\omega_i$ on the attainable damping, the Ro wis = [0.01, 0.1, 0.5, 1]*w0; % [rad/s] #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none figure; gains = logspace(-2, 4, 100); @@ -975,7 +975,7 @@ To find the optimum, the gain that maximize the simultaneous damping of the mode end #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none figure; yyaxis left plot(wis, opt_zeta, '-o', 'DisplayName', '$\xi_{cl}$'); @@ -1040,7 +1040,7 @@ Let's define initial values for the model. w0 = sqrt(k/m); % [rad/s] #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none kp = 0; % [N/m] cp = 0; % [N/(m/s)] #+end_src @@ -1080,7 +1080,7 @@ The rotation speed is set to $\Omega = 0.1 \omega_0$. W = 0.1*w0; % [rad/s] #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none Kiff = tf(zeros(2)); Kdvf = tf(zeros(2)); #+end_src @@ -1094,7 +1094,7 @@ The results are shown in Figure [[fig:plant_iff_kp]]. One can see that for $k_p > m \Omega^2$, the systems shows alternating complex conjugate poles and zeros. -#+begin_src matlab :exports code +#+begin_src matlab :exports none %% Name of the Simulink File mdl = 'rotating_frame'; @@ -1137,7 +1137,7 @@ One can see that for $k_p > m \Omega^2$, the systems shows alternating complex c Giff_l.OutputName = {'fu', 'fv'}; #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none freqs = logspace(-2, 1, 1000); figure; @@ -1194,7 +1194,7 @@ Thus, decentralized IFF controller with pure integrators can be used if: k_{p} > m \Omega^2 \end{equation} -#+begin_src matlab :exports code +#+begin_src matlab :exports none figure; gains = logspace(-2, 2, 100); @@ -1306,7 +1306,7 @@ To study the second point, Root Locus plots for the following values of $k_p$ ar It is shown that large values of $k_p$ decreases the attainable damping. -#+begin_src matlab :exports code +#+begin_src matlab :exports none figure; gains = logspace(-2, 4, 100); @@ -1375,7 +1375,7 @@ Let's take $k_p = 5 m \Omega^2$ and find the optimal IFF control gain $g$ such t end #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none figure; gains = logspace(-2, 4, 100); @@ -1458,7 +1458,7 @@ Let's define initial values for the model. w0 = sqrt(k/m); % [rad/s] #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none kp = 0; % [N/m] cp = 0; % [N/(m/s)] #+end_src @@ -1469,7 +1469,7 @@ The rotating speed is set to $\Omega = 0.1 \omega_0$. W = 0.1*w0; #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none Kiff = tf(zeros(2)); Kdvf = tf(zeros(2)); #+end_src @@ -1506,7 +1506,7 @@ The same transfer function from $[F_u, F_v]$ to $[v_u, v_v]$ is written down fro The two are compared in Figure [[fig:plant_iff_comp_simscape_analytical]] and found to perfectly match. -#+begin_src matlab :exports code +#+begin_src matlab :exports none freqs = logspace(-1, 1, 1000); figure; @@ -1579,7 +1579,7 @@ The corresponding Root Locus plots for the following rotating speeds are shown i It is shown that for rotating speed $\Omega < \omega_0$, the closed loop system is unconditionally stable and arbitrary damping can be added to the poles. -#+begin_src matlab :exports code +#+begin_src matlab :exports none gains = logspace(-2, 1, 100); figure; @@ -1659,12 +1659,12 @@ Let's define initial values for the model. w0 = sqrt(k/m); % [rad/s] #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none kp = 0; % [N/m] cp = 0; % [N/(m/s)] #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none Kiff = tf(zeros(2)); Kdvf = tf(zeros(2)); #+end_src @@ -1676,7 +1676,7 @@ The rotating speed is set to $\Omega = 0.1 \omega_0$. ** Root Locus *** Pseudo Integrator IFF :ignore: -#+begin_src matlab :exports code +#+begin_src matlab :exports none kp = 0; cp = 0; #+end_src @@ -1728,7 +1728,7 @@ The rotating speed is set to $\Omega = 0.1 \omega_0$. #+end_src *** DVF :ignore: -#+begin_src matlab :exports code +#+begin_src matlab :exports none kp = 0; cp = 0; #+end_src @@ -1752,7 +1752,7 @@ The rotating speed is set to $\Omega = 0.1 \omega_0$. #+end_src *** Root Locus :ignore: -#+begin_src matlab :exports code +#+begin_src matlab :exports none figure; gains = logspace(-2, 2, 100); @@ -1819,7 +1819,7 @@ The rotating speed is set to $\Omega = 0.1 \omega_0$. ** 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. -#+begin_src matlab :exports code +#+begin_src matlab :exports none %% IFF with pseudo integrators gains = linspace(0, (w0^2/W^2 - 1)*wi, 100); opt_zeta_iff = 0; @@ -1837,7 +1837,7 @@ In order to compare to three considered Active Damping techniques, gains that yi end #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none %% IFF with Parallel Stiffness gains = logspace(-2, 4, 100); opt_zeta_kp = 0; @@ -1855,7 +1855,7 @@ In order to compare to three considered Active Damping techniques, gains that yi end #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none %% Direct Velocity Feedback gains = logspace(0, 2, 100); opt_zeta_dvf = 0; @@ -1889,7 +1889,7 @@ The obtained damping ratio and control are shown below. ** Transmissibility <> *** Open Loop :ignore: -#+begin_src matlab :exports code +#+begin_src matlab :exports none Kdvf = tf(zeros(2)); Kiff = tf(zeros(2)); @@ -1916,7 +1916,7 @@ The obtained damping ratio and control are shown below. #+end_src *** Pseudo Integrator IFF :ignore: -#+begin_src matlab :exports code +#+begin_src matlab :exports none kp = 0; cp = 0; @@ -1955,7 +1955,7 @@ The obtained damping ratio and control are shown below. Kiff = opt_gain_kp/s*tf(eye(2)); #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none Kdvf = tf(zeros(2)); #+end_src @@ -1978,7 +1978,7 @@ The obtained damping ratio and control are shown below. #+end_src *** DVF :ignore: -#+begin_src matlab :exports code +#+begin_src matlab :exports none kp = 0; cp = 0; @@ -2008,7 +2008,7 @@ The obtained damping ratio and control are shown below. #+end_src *** Transmissibility :ignore: -#+begin_src matlab :exports code +#+begin_src matlab :exports none freqs = logspace(-2, 1, 1000); figure; @@ -2039,7 +2039,7 @@ The obtained damping ratio and control are shown below. ** Compliance <> *** Open Loop :ignore: -#+begin_src matlab :exports code +#+begin_src matlab :exports none Kdvf = tf(zeros(2)); Kiff = tf(zeros(2)); @@ -2066,7 +2066,7 @@ The obtained damping ratio and control are shown below. #+end_src *** Pseudo Integrator IFF :ignore: -#+begin_src matlab :exports code +#+begin_src matlab :exports none kp = 0; cp = 0; @@ -2095,7 +2095,7 @@ The obtained damping ratio and control are shown below. Kiff = opt_gain_kp/s*tf(eye(2)); #+end_src -#+begin_src matlab :exports code +#+begin_src matlab :exports none Kdvf = tf(zeros(2)); #+end_src @@ -2108,7 +2108,7 @@ The obtained damping ratio and control are shown below. #+end_src *** DVF :ignore: -#+begin_src matlab :exports code +#+begin_src matlab :exports none kp = 0; cp = 0; @@ -2128,7 +2128,7 @@ The obtained damping ratio and control are shown below. #+end_src *** Compliance :ignore: -#+begin_src matlab :exports code +#+begin_src matlab :exports none freqs = logspace(-2, 1, 1000); figure;