1 System Description and Analysis
+
+
1 System Description and Analysis
-
-
1.1 System description
+
+
-
1.1 System description
-The system consists of one 2 degree of freedom translation stage on top of a spindle (figure 1). +The system consists of one 2 degree of freedom translation stage on top of a spindle (figure 1).
-
+
Figure 1: Figure caption
+Figure 1: Schematic of the studied system
@@ -141,33 +163,48 @@ The measurement is either the \(x-y\) displacement of the object located on top
-
1.2 Equations
+
+
-1.2 Equations
-Based on the Figure 1, the equations of motions are: +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}
+\bm{G}_d
+\begin{bmatrix} F_u \\ F_v \end{bmatrix}
+\end{equation}
+
+\begin{equation}
+\begin{bmatrix} d_u \\ d_v \end{bmatrix} =
+\frac{1}{k} \frac{1}{G_{dp}}
\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} \\
+ G_{dz} & G_{dc} \\
+ -G_{dc} & G_{dz}
\end{bmatrix}
\begin{bmatrix} F_u \\ F_v \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}Explain Coriolis and Centrifugal Forces (negative Stiffness) +=> First write the equations in terms of \(k\), \(m\) and \(c\) and explain the terms.
-
1.3 Numerical Values
+
+
-1.3 Numerical Values
Let’s define initial values for the model. @@ -187,19 +224,19 @@ w0 = sqrt(k/m); % [rad/s]
-
1.4 Campbell Diagram
+
+
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). +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).
-
+
-
Figure 2: Campbell Diagram
@@ -207,8 +244,8 @@ It is shown in Figure 2, and one can see that the syst
-
1.5 Simscape Model
+
+
-1.5 Simscape Model
Define the rotating speed for the Simscape Model. @@ -244,8 +281,8 @@ G.OutputName = {'du', 'dv'};
-
1.6 Comparison of the Analytical Model and the Simscape Model
+
+
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. @@ -258,11 +295,11 @@ The same transfer function from \([F_u, F_v]\) to \([d_u, d_v]\) is written down
-Both transfer functions are compared in Figure 3 and are found to perfectly match. +Both transfer functions are compared in Figure 3 and are found to perfectly match.
-
+
-
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
@@ -270,8 +307,8 @@ Both transfer functions are compared in Figure 3 and a
-
1.7 Effect of the rotation speed
+
+
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. @@ -295,11 +332,11 @@ end
-They are compared in Figure 4. +They are compared in Figure 4.
-
+
-
Figure 4: Comparison of the transfer functions from \([F_u, F_v]\) to \([d_u, d_v]\) for several rotating speed
@@ -308,11 +345,11 @@ They are compared in Figure 4.
-
2 Problem with pure Integral Force Feedback
+
+
-
2 Problem with pure Integral Force Feedback
- Diagram with the controller @@ -320,8 +357,8 @@ They are compared in Figure 4.
-
2.1 Plant Parameters
+
+
-2.1 Plant Parameters
Let’s define initial values for the model. @@ -341,8 +378,8 @@ w0 = sqrt(k/m); % [rad/s]
-
2.2 Equations
+
+
-2.2 Equations
The sensed forces are equal to: @@ -363,20 +400,32 @@ 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}
+\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}
- (\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 \\
+ 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}
+
-
2.3 Simscape Model
+
+
-2.3 Simscape Model
The rotation speed is set to \(\Omega = 0.1 \omega_0\). @@ -411,8 +460,8 @@ Giff.OutputName = {'fu', 'fv'};
-
2.4 Comparison of the Analytical Model and the Simscape Model
+
+
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. @@ -425,11 +474,11 @@ The same transfer function from \([F_u, F_v]\) to \([f_u, f_v]\) is written down
-The two are compared in Figure 5 and found to perfectly match. +The two are compared in Figure 5 and found to perfectly match.
-
+
-
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
@@ -437,8 +486,8 @@ The two are compared in Figure 5 and found to perfectl
-
2.5 Effect of the rotation speed
+
+
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. @@ -462,10 +511,10 @@ end
-The obtained transfer functions are shown in Figure 6. +The obtained transfer functions are shown in Figure 6.
-
+
-
Figure 6: Comparison of the transfer functions from \([F_u, F_v]\) to \([f_u, f_v]\) for several rotating speed
@@ -473,8 +522,8 @@ The obtained transfer functions are shown in Figure 6.
-
+
+2.6 Decentralized Integral Force Feedback
+
+
@@ -510,33 +558,191 @@ Kiff = g/s*tf(eye(2));
2.6 Decentralized Integral Force Feedback
Let’s take \(\Omega = \frac{\omega_0}{10}\). @@ -488,7 +537,6 @@ Let’s take \(\Omega = \frac{\omega_0}{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 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. +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.
-
+
Figure 7: Root Locus for the Decentralized Integral Force Feedback controller. Several rotating speed are shown.
+
+
+2.7 Sort poles
+
+
+
+
-W = 0.1; + +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]; +g = 1; + +gains = logspace(-2, 4, 100); + +poles = zeros(length(pole(minreal(feedback(Giff, g/s*eye(2))))), length(gains)); +poles(:, 1) = pole(minreal(feedback(Giff, gains(1)/s*eye(2)))); +poles(:, 2) = pole(minreal(feedback(Giff, gains(2)/s*eye(2)))); +% poles(:, 2) = poles(:, 1); + +for g_i = 3:length(gains) + poles_est = poles(:, g_i-1) + (poles(:, g_i-1) - poles(:, g_i-2))*(gains(g_i) - gains(g_i-1))/(gains(g_i-1) - gains(g_i - 2)); + poles_est = poles(:, g_i-1); + + poles_gi = pole(minreal(feedback(Giff, gains(g_i)/s*eye(2)))); + + % Array of distances between all the poles + poles_dist = (poles_est-poles_gi.').*conj(poles_est-poles_gi.'); + [~, c] = sort(min(poles_dist)); + poles_dist = poles_dist(:, c); + + for p_i = 1:size(poles_dist, 1) + [~, a_i] = min(poles_dist(:, p_i)); + + poles(c(p_i), g_i) = poles_gi(a_i); + poles_dist(a_i, :) = []; + poles_gi(a_i, :) = []; + end +end +
-
+3 Modified IFF (pseudo integrator)
++Working +
+
+
+
+W = 0.1; + +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]; + +gains = logspace(-2, 4, 500); + +poles = zeros(length(pole(feedback(Giff, 1/s*eye(2)))), length(gains)); +poles(:, 1) = pole(feedback(Giff, gains(1)/s*eye(2))); +poles(:, 2) = pole(feedback(Giff, gains(2)/s*eye(2))); +% poles(:, 2) = poles(:, 1); + +for g_i = 3:length(gains) + % poles_est = poles(:, g_i-1) + (poles(:, g_i-1) - poles(:, g_i-2))*(gains(g_i) - gains(g_i-1))/(gains(g_i-1) - gains(g_i - 2)); + poles_est = poles(:, g_i-1); + poles_gi = pole(feedback(Giff, gains(g_i)/s*eye(2))); + + % Array of distances between all the poles + poles_dist = sqrt((poles_est-poles_gi.').*conj(poles_est-poles_gi.')); + + % Columns of distances are sorted from lowest to highest + [~, c] = sort(min(poles_dist)); + poles_dist = poles_dist(:, c); + + % for each column of poles_dist corresponding to the i'th pole + % with closest previous poles + for p_i = 1:size(poles_dist, 2) + % Get the indice a_i of the previous pole that is the closest + % to pole c(p_i) + [~, a_i] = min(poles_dist(:, p_i)); + + poles(a_i, g_i) = poles_gi(c(p_i)); + + % poles_dist(a_i, :) = []; + % poles_gi(a_i, :) = []; + end +end ++
+
+
+W = 0.1; + +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]; + +gains = logspace(-2, 4, 500); + +poles = zeros(length(pole(feedback(Giff, 1/s*eye(2)))), length(gains)); +poles(:, 1) = pole(feedback(Giff, gains(1)/s*eye(2))); +poles(:, 2) = pole(feedback(Giff, gains(2)/s*eye(2))); + +for g_i = 3:length(gains) + % Estimated value of the poles + poles_est = poles(:, g_i-1) + (poles(:, g_i-1) - poles(:, g_i-2))*(gains(g_i) - gains(g_i-1))/(gains(g_i-1) - gains(g_i - 2)); + + % New values for the poles + poles_gi = pole(feedback(Giff, gains(g_i)/s*eye(2))); + + % Array of distances between all the poles + poles_dist = sqrt((poles_est-poles_gi.').*conj(poles_est-poles_gi.')); + + % Get indices corresponding to distances from lowest to highest + [~, c] = sort(min(poles_dist)); + + as = 1:length(poles_gi); + + % for each column of poles_dist corresponding to the i'th pole + % with closest previous poles + for p_i = c + % Get the indice a_i of the previous pole that is the closest + % to pole c(p_i) + [~, a_i] = min(poles_dist(:, p_i)); + + poles(as(a_i), g_i) = poles_gi(p_i); + + % Remove old poles that are already matched + % poles_gi(as(a_i), :) = []; + poles_dist(a_i, :) = []; + as(a_i) = []; + end +end ++
+How to remove poles that are not moving? +
+
+
+
+poles_rl = poles(max(abs(poles(:, 2:end) - poles(:, 1:end-1))') > 1e-8, :); +poles_rl = poles_rl(1:end/2, :); ++
-
+
- works best without minreal +
[ ]create a function
+
+
+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]; +Kiff = 1/s*eye(2); +gains = logspace(-2, 4, 500); + +[poles] = rootLocusPolesSorted(Giff, Kiff, gains, 'd_max', 1e-4); ++
+
3 Modified IFF (pseudo integrator)
-
-
3.1 Plant Parameters
+
+
-3.1 Plant Parameters
Let’s define initial values for the model. @@ -556,8 +762,8 @@ w0 = sqrt(k/m); % [rad/s]
-
3.2 Modified Integral Force Feedback Controller
+
+
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: @@ -592,10 +798,10 @@ And the following rotating speed.
-The obtained Loop Gain is shown in Figure 8. +The obtained Loop Gain is shown in Figure 8.
-
+
-
Figure 8: Loop Gain for the modified IFF controller
@@ -603,15 +809,15 @@ The obtained Loop Gain is shown in Figure 8.
-
3.3 Root Locus
+
+
3.3 Root Locus
-As shown in the Root Locus plot (Figure 9), for some value of the gain, the system remains stable. +As shown in the Root Locus plot (Figure 9), for some value of the gain, the system remains stable.
-
+
-
Figure 9: Root Locus for the modified IFF controller
@@ -619,11 +825,11 @@ As shown in the Root Locus plot (Figure 9), for some v
-
3.4 What is the optimal \(\omega_i\) and \(g\)?
+
+
-
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\): +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] @@ -631,7 +837,7 @@ In order to visualize the effect of \(\omega_i\) on the attainable damping, the
+
-
Figure 10: Root Locus for the modified IFF controller (zoomed plot on the left)
@@ -650,18 +856,17 @@ For the controller 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) + 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 10), the control gains may be limited to small values due to \eqref{eq:iff_gmax} thus reducing the attainable damping. +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). +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] @@ -688,7 +893,7 @@ end
+
-
Figure 11: Simultaneous attainable damping of the closed loop poles as a function of \(\omega_i\)
@@ -698,17 +903,80 @@ end
-
4 IFF with a stiffness in parallel with the force sensor
+
+
4 IFF with a stiffness in parallel with the force sensor
-
-
4.1 Plant Parameters
+
+
+
+4.1 Equations
+
+
+\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^\prime}^2} + 2\xi^\prime \frac{s}{{\omega_0^\prime}^2} + 1 - \frac{\Omega}{\omega_0^\prime} \right)^2 + \left( 2 \frac{\Omega}{\omega_0^\prime}\frac{s}{\omega_0^\prime} \right)^2 \\ + G_{kz} &= \left( \frac{s^2}{{\omega_0^\prime}^2} + \frac{k_p}{k + k_p} - \frac{\Omega^2}{{\omega_0^\prime}^2} \right) \left( \frac{s^2}{{\omega_0^\prime}^2} + 2\xi^\prime \frac{s}{{\omega_0^\prime}^2} + 1 - \frac{\Omega}{\omega_0^\prime} \right) + \left( 2 \frac{\Omega}{\omega_0^\prime}\frac{s}{\omega_0^\prime} \right)^2 \\ + G_{kc} &= \left( 2 \xi^\prime \frac{s}{\omega_0^\prime} + \frac{k}{k + k_p} \right) \left( 2 \frac{\Omega}{\omega_0^\prime}\frac{s}{\omega_0^\prime} \right) +\end{align} ++where: +
+-
+
- \(\omega_0^\prime = \frac{k + k_p}{m}\) +
- \(\xi^\prime = \frac{c}{2 \sqrt{(k + k_p) m}}\) +
+
+
+Giff_th = 1/( (m*s^2 + c*s + k + kp - m*W^2)^2 + (2*m*s*W)^2 )*[... + (m*s^2 + c*s + k + kp - m*W^2)^2 + (2*m*s*W)^2 - (c*s + k)*(m*s^2 + c*s + k + kp - m*W^2), -(c*s + k)*(2*m*s*W); + (c*s + k)*(2*m*s*W), (m*s^2 + c*s + k + kp - m*W^2)^2 + (2*m*s*W)^2 - (c*s + k)*(m*s^2 + c*s + k + kp - m*W^2) + ]; ++
+
+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_b = 1/( (m*s^2 + c*s + k + kp - m*W^2)^2 + (2*m*s*W)^2 )*[... + (m*s^2 + c*s + k + kp - m*W^2)^2 + (2*m*s*W)^2 - (c*s + k)*(m*s^2 + c*s + k + kp - m*W^2), -(c*s + k)*(2*m*s*W); + (c*s + k)*(2*m*s*W), (m*s^2 + c*s + k + kp - m*W^2)^2 + (2*m*s*W)^2 - (c*s + k)*(m*s^2 + c*s + k + kp - m*W^2) + ]; ++
+
-4.2 Plant Parameters
+Let’s define initial values for the model.
@@ -727,11 +995,11 @@ w0 = sqrt(k/m); % [rad/s]
-
4.2 Schematic
-
+
+
4.3 Schematic
+
-
+
-
Figure 12: Figure caption
@@ -739,9 +1007,9 @@ w0 = sqrt(k/m); % [rad/s]
-
4.3 Physical Explanation
-
+
+
-4.4 Physical Explanation
+- Negative stiffness induced by gyroscopic effects
- Zeros of the open-loop <=> Poles of the subsystem with the force sensors removes @@ -766,22 +1034,8 @@ And thus, the stiffness in parallel should be such that:
-
-
-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
+
+
-
4.5 Effect of the parallel stiffness on the IFF plant
The rotation speed is set to \(\Omega = 0.1 \omega_0\). @@ -801,7 +1055,7 @@ And the IFF plant (transfer function from \([F_u, F_v]\) to \([f_u, f_v]\)) is i
-The results are shown in Figure 13. +The results are shown in Figure 13.
@@ -822,7 +1076,7 @@ Giff.OutputName = {'fu', 'fv'};
kp = 0.5*m*W^2;
-cp = 0.001;
+cp = 0;
Giff_s = linearize(mdl, io, 0);
@@ -834,7 +1088,7 @@ Giff_s.OutputName = {'fu', 'fv'};
kp = 1.5*m*W^2;
-cp = 0.001;
+cp = 0;
Giff_l = linearize(mdl, io, 0);
@@ -845,18 +1099,19 @@ Giff_l.OutputName = {'fu', 'fv'};
-
+
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
+ +
+
4.6 IFF when adding a spring in parallel
-In Figure 14 is displayed the Root Locus in the three considered cases with +In Figure 14 is displayed the Root Locus in the three considered cases with
\begin{equation} K_{\text{IFF}} = \frac{g}{s} \begin{bmatrix} @@ -877,7 +1132,7 @@ Thus, decentralized IFF controller with pure integrators can be used if: \end{equation} -
+
-
Figure 14: Root Locus
@@ -885,8 +1140,8 @@ Thus, decentralized IFF controller with pure integrators can be used if:
-
4.7 Effect of \(k_p\) on the attainable damping
+
+
-
4.7 Effect of \(k_p\) on the attainable damping
However, having large values of \(k_p\) may: @@ -897,7 +1152,7 @@ However, having large values of \(k_p\) may:
-To study the second point, Root Locus plots for the following values of \(k_p\) are shown in Figure 15. +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; @@ -910,15 +1165,16 @@ It is shown that large values of \(k_p\) decreases the attainable damping. -++
Figure 15: Root Locus plot
-
-4.8 Optimal Gain
+
+
-
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. @@ -952,24 +1208,25 @@ end
+
Figure 16: Root Locus for \(k_p = 5 m \Omega^2\) and the poles corresponding to the identified optimal gain
-
5 Direct Velocity Feedback
+
+
5 Direct Velocity Feedback
-
-
5.1 Equations
+
+
-5.1 Equations
The sensed relative velocity are equal to: @@ -977,20 +1234,34 @@ 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}
+\bm{G}_v
\begin{bmatrix} F_u \\ F_v \end{bmatrix}
\end{equation}
+\begin{equation}
+\begin{bmatrix} \dot{d}_u \\ \dot{d}_v \end{bmatrix} =
+\frac{s}{k} \frac{1}{G_{vp}}
+\begin{bmatrix}
+ G_{vz} & G_{vc} \\
+ -G_{vc} & G_{vz}
+\end{bmatrix}
+\begin{bmatrix} F_u \\ F_v \end{bmatrix}
+\end{equation}
+
+With: +
+\begin{align} + G_{vp} &= \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_{vz} &= \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} \\ + G_{vc} &= 2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0} +\end{align} +
-
5.2 Plant Parameters
+
+
-5.2 Plant Parameters
Let’s define initial values for the model. @@ -1010,8 +1281,8 @@ w0 = sqrt(k/m); % [rad/s]
-
5.3 Plant - Bode Plot
+
+
-5.3 Plant - Bode Plot
The rotating speed is set to \(\Omega = 0.1 \omega_0\). @@ -1046,8 +1317,8 @@ Gdvf.OutputName = {'Vu', 'Vv'};
-
5.4 Comparison of the Analytical Model and the Simscape Model
+
+
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. @@ -1063,11 +1334,11 @@ Gdvf_th.OutputName = {'vu', 'vv'};
-The two are compared in Figure 5 and found to perfectly match. +The two are compared in Figure 5 and found to perfectly match.
-
+
-
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
@@ -1075,8 +1346,8 @@ The two are compared in Figure 5 and found to perfectl
-
5.5 Root Locus
+
+
+
5.5 Root Locus
The Decentralized Direct Velocity Feedback controller consist of a pure gain on the diagonal: @@ -1089,7 +1360,7 @@ The Decentralized Direct Velocity Feedback controller consist of a pure gain on \end{equation}
-The corresponding Root Locus plots for the following rotating speeds are shown in Figure 18. +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] @@ -1099,9 +1370,17 @@ The corresponding Root Locus plots for the following rotating speeds are shown iIt 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.
+++gains = logspace(-2, 2, 100); ++
+
-5.5.0.1 Root Locus Plots
+
-
+
+
Figure 18: Root Locus for the Decentralized Direct Velocity Feedback controller. Several rotating speed are shown.
@@ -1109,17 +1388,18 @@ It is shown that for rotating speed \(\Omega < \omega_0\), the closed loop syste
-
6 Comparison
+
+
6 Comparison
-
-
6.1 Plant Parameters
+
+
-6.1 Plant Parameters
Let’s define initial values for the model. @@ -1147,8 +1427,8 @@ The rotating speed is set to \(\Omega = 0.1 \omega_0\).
-
6.2 Root Locus
+
+
-
-
6.2 Root Locus
wi = 0.1*w0;
@@ -1219,15 +1499,16 @@ Gdvf.OutputName = {'Vu', 'Vv'};
+
Figure 19: Root Locus plot - Comparison of IFF with additional high pass filter, IFF with additional parallel stiffness and DVF
-
6.3 Controllers - Optimal Gains
+
+
-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. @@ -1277,11 +1558,11 @@ The obtained damping ratio and control are shown below.
-
6.4 Transmissibility
+
+
-
-
6.4 Transmissibility
%% Name of the Simulink File
@@ -1381,18 +1662,19 @@ Tdvf.OutputName = {'Dx', 'Dy'};
+
Figure 20: Comparison of the transmissibility
-
6.5 Compliance
+
+
->
@@ -383,10 +407,7 @@ They are compared in Figure [[fig:plant_compare_rotating_speed]].
#+begin_src matlab
addpath('./matlab/');
-#+end_src
-
-#+begin_src matlab
- open('rotating_frame.slx');
+ addpath('./src/');
#+end_src
** Plant Parameters
@@ -423,13 +444,25 @@ Which then gives:
#+begin_important
\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}
+\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}
- (\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 \\
+ 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}
#+end_important
** Simscape Model
@@ -443,6 +476,10 @@ The rotation speed is set to $\Omega = 0.1 \omega_0$.
Kdvf = tf(zeros(2));
#+end_src
+#+begin_src matlab
+ open('rotating_frame.slx');
+#+end_src
+
And the transfer function from $[F_u, F_v]$ to $[f_u, f_v]$ is identified using the Simscape model.
#+begin_src matlab
%% Name of the Simulink File
@@ -534,7 +571,7 @@ The two are compared in Figure [[fig:plant_iff_comp_simscape_analytical]] and fo
** 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.
#+begin_src matlab
- Ws = [0, 0.1, 0.5, 0.8, 1.1]*w0; % Rotating Speeds [rad/s]
+ Ws = [0, 0.2, 0.7, 1.1]*w0; % Rotating Speeds [rad/s]
#+end_src
#+begin_src matlab
@@ -559,7 +596,7 @@ The obtained transfer functions are shown in Figure [[fig:plant_iff_compare_rota
hold on;
for W_i = 1:length(Ws)
plot(freqs, abs(squeeze(freqresp(Gsiff{W_i}(1,1), freqs))), ...
- 'DisplayName', sprintf('$\\omega = %.1f \\omega_0 $', Ws(W_i)/w0))
+ 'DisplayName', sprintf('$\\Omega = %.1f \\omega_0 $', Ws(W_i)/w0))
end
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
@@ -590,19 +627,11 @@ The obtained transfer functions are shown in Figure [[fig:plant_iff_compare_rota
#+RESULTS:
[[file:figs/plant_iff_compare_rotating_speed.png]]
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/plant_iff_compare_rotating_speed.pdf', 'width', 'full', 'height', 'full', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
** Decentralized Integral Force Feedback
-Let's take $\Omega = \frac{\omega_0}{10}$.
-#+begin_src matlab
- W = w0/10;
-#+end_src
-
-#+begin_src 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];
-
-#+end_src
-
The decentralized IFF controller consists of pure integrators:
\begin{equation}
\bm{K}_{\text{IFF}}(s) = \frac{g}{s} \begin{bmatrix}
@@ -611,44 +640,9 @@ The decentralized IFF controller consists of pure integrators:
\end{bmatrix}
\end{equation}
-#+begin_src matlab
- g = 2;
-
- Kiff = g/s*tf(eye(2));
-#+end_src
-
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 [[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.
-#+begin_src matlab :exports none
- freqs = logspace(-2, 1, 1000);
-
- figure;
-
- ax1 = subplot(2, 1, 1);
- hold on;
- for W_i = 1:length(Ws)
- plot(freqs, abs(squeeze(freqresp(Gsiff{W_i}(1,1)*Kiff(1,1), freqs))))
- end
- hold off;
- set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
- set(gca, 'XTickLabel',[]); ylabel('Loop Gain');
-
- ax2 = subplot(2, 1, 2);
- hold on;
- for W_i = 1:length(Ws)
- plot(freqs, 180/pi*angle(squeeze(freqresp(Gsiff{W_i}(1,1)*Kiff(1,1), freqs))))
- end
- set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
- xlabel('Frequency [rad/s]'); ylabel('Phase [deg]');
- yticks(-180:90:180);
- ylim([-180 180]);
- hold off;
-
- linkaxes([ax1,ax2],'x');
- xlim([freqs(1), freqs(end)]);
-#+end_src
-
#+begin_src matlab :exports none
figure;
@@ -658,7 +652,7 @@ It is shown that for non-null rotating speed, one pole is bound to the right-hal
for W_i = 1:length(Ws)
set(gca,'ColorOrderIndex',W_i);
plot(real(pole(Gsiff{W_i})), imag(pole(Gsiff{W_i})), 'x', ...
- 'DisplayName', sprintf('$\\omega = %.1f \\omega_0 $', Ws(W_i)/w0));
+ 'DisplayName', sprintf('$\\Omega = %.1f \\omega_0 $', Ws(W_i)/w0));
set(gca,'ColorOrderIndex',W_i);
plot(real(tzero(Gsiff{W_i})), imag(tzero(Gsiff{W_i})), 'o', ...
'HandleVisibility', 'off');
@@ -686,7 +680,39 @@ It is shown that for non-null rotating speed, one pole is bound to the right-hal
#+RESULTS:
[[file:figs/root_locus_pure_iff.png]]
-* Modified IFF (pseudo integrator)
+#+begin_src matlab :exports none :tangle no
+ gains = logspace(-2, 4, 1000);
+
+ figure;
+ hold on;
+ for W_i = 1:length(Ws)
+
+ set(gca,'ColorOrderIndex',W_i);
+ plot(real(pole(Gsiff{W_i})), imag(pole(Gsiff{W_i})), 'x', ...
+ 'DisplayName', sprintf('$\\Omega = %.1f \\omega_0 $', Ws(W_i)/w0));
+ set(gca,'ColorOrderIndex',W_i);
+ plot(real(tzero(Gsiff{W_i})), imag(tzero(Gsiff{W_i})), 'o', ...
+ 'HandleVisibility', 'off');
+ poles = rootLocusPolesSorted(Gsiff{W_i}, 1/s*eye(2), gains, 'd_max', 1e-4);
+ for p_i = 1:size(poles, 2)
+ set(gca,'ColorOrderIndex',W_i);
+ plot(real(poles(:, p_i)), imag(poles(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+ end
+ hold off;
+ axis square;
+ xlim([-2, 0.5]); ylim([0, 2.5]);
+
+ xlabel('Real Part'); ylabel('Imaginary Part');
+ legend('location', 'northwest');
+#+end_src
+
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/root_locus_pure_iff.pdf', 'width', 'wide', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
+* Integral Force Feedback with an High Pass Filter
<>
** Introduction :ignore:
@@ -703,10 +729,7 @@ It is shown that for non-null rotating speed, one pole is bound to the right-hal
#+begin_src matlab
addpath('./matlab/');
-#+end_src
-
-#+begin_src matlab
- open('rotating_frame.slx');
+ addpath('./src/');
#+end_src
** Plant Parameters
@@ -775,9 +798,9 @@ The obtained Loop Gain is shown in Figure [[fig:loop_gain_modified_iff]].
ax2 = subplot(2, 1, 2);
hold on;
plot(freqs, 180/pi*angle(squeeze(freqresp(Giff(1,1)*(g/s), freqs))), ...
- 'DisplayName', 'Pure Integrator')
+ 'DisplayName', 'IFF')
plot(freqs, 180/pi*angle(squeeze(freqresp(Giff(1,1)*Kiff(1,1), freqs))), ...
- 'DisplayName', 'Pseudo Integrator')
+ 'DisplayName', 'IFF + HPF')
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [rad/s]'); ylabel('Phase [deg]');
yticks(-180:90:180);
@@ -798,6 +821,10 @@ The obtained Loop Gain is shown in Figure [[fig:loop_gain_modified_iff]].
#+RESULTS:
[[file:figs/loop_gain_modified_iff.png]]
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/loop_gain_modified_iff.pdf', 'width', 'full', 'height', 'full', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
** 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.
@@ -810,7 +837,7 @@ As shown in the Root Locus plot (Figure [[fig:root_locus_modified_iff]]), for so
hold on;
% Pure Integrator
set(gca,'ColorOrderIndex',1);
- plot(real(pole(Giff)), imag(pole(Giff)), 'x', 'DisplayName', 'Pure Int');
+ plot(real(pole(Giff)), imag(pole(Giff)), 'x', 'DisplayName', 'IFF');
set(gca,'ColorOrderIndex',1);
plot(real(tzero(Giff)), imag(tzero(Giff)), 'o', 'HandleVisibility', 'off');
for g = gains
@@ -820,7 +847,7 @@ As shown in the Root Locus plot (Figure [[fig:root_locus_modified_iff]]), for so
end
% Modified IFF
set(gca,'ColorOrderIndex',2);
- plot(real(pole(Giff)), imag(pole(Giff)), 'x', 'DisplayName', 'Pseudo Int');
+ plot(real(pole(Giff)), imag(pole(Giff)), 'x', 'DisplayName', 'IFF + HPF');
set(gca,'ColorOrderIndex',2);
plot(real(tzero(Giff)), imag(tzero(Giff)), 'o', 'HandleVisibility', 'off');
for g = gains
@@ -871,6 +898,74 @@ As shown in the Root Locus plot (Figure [[fig:root_locus_modified_iff]]), for so
#+RESULTS:
[[file:figs/root_locus_modified_iff.png]]
+#+begin_src matlab :exports none :tangle no
+ gains = logspace(-2, 3, 200);
+
+ poles_iff = rootLocusPolesSorted(Giff, 1/s*eye(2), gains, 'd_max', 1e-4);
+ poles_iff_hpf = rootLocusPolesSorted(Giff, 1/(s + wi)*eye(2), gains, 'd_max', 1e-4);
+
+ figure;
+
+ ax1 = subplot(1, 2, 1);
+ hold on;
+ % Pure Integrator
+ set(gca,'ColorOrderIndex',1);
+ plot(real(pole(Giff)), imag(pole(Giff)), 'x', 'DisplayName', 'IFF');
+ set(gca,'ColorOrderIndex',1);
+ plot(real(tzero(Giff)), imag(tzero(Giff)), 'o', 'HandleVisibility', 'off');
+ for p_i = 1:size(poles_iff, 2)
+ set(gca,'ColorOrderIndex',1);
+ plot(real(poles_iff(:, p_i)), imag(poles_iff(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+ % Modified IFF
+ set(gca,'ColorOrderIndex',2);
+ plot(real(pole(Giff)), imag(pole(Giff)), 'x', 'DisplayName', 'IFF + HPF');
+ set(gca,'ColorOrderIndex',2);
+ plot(real(tzero(Giff)), imag(tzero(Giff)), 'o', 'HandleVisibility', 'off');
+ for p_i = 1:size(poles_iff_hpf, 2)
+ set(gca,'ColorOrderIndex',2);
+ plot(real(poles_iff_hpf(:, p_i)), imag(poles_iff_hpf(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+ hold off;
+ axis square;
+ xlim([-2, 0.5]); ylim([-1.25, 1.25]);
+ legend('location', 'northwest');
+ xlabel('Real Part'); ylabel('Imaginary Part');
+
+ ax2 = subplot(1, 2, 2);
+ hold on;
+ % Pure Integrator
+ set(gca,'ColorOrderIndex',1);
+ plot(real(pole(Giff)), imag(pole(Giff)), 'x');
+ set(gca,'ColorOrderIndex',1);
+ plot(real(tzero(Giff)), imag(tzero(Giff)), 'o');
+ for p_i = 1:size(poles_iff, 2)
+ set(gca,'ColorOrderIndex',1);
+ plot(real(poles_iff(:, p_i)), imag(poles_iff(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+ % Modified IFF
+ set(gca,'ColorOrderIndex',2);
+ plot(real(pole(Giff)), imag(pole(Giff)), 'x');
+ set(gca,'ColorOrderIndex',2);
+ plot(real(tzero(Giff)), imag(tzero(Giff)), 'o');
+ for p_i = 1:size(poles_iff_hpf, 2)
+ set(gca,'ColorOrderIndex',2);
+ plot(real(poles_iff_hpf(:, p_i)), imag(poles_iff_hpf(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+ hold off;
+ axis square;
+ xlim([-0.2, 0.1]); ylim([-0.15, 0.15]);
+ xlabel('Real Part'); ylabel('Imaginary Part');
+#+end_src
+
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/root_locus_modified_iff.pdf', 'width', 'full', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
** 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 [[fig:root_locus_wi_modified_iff]] for the following $\omega_i$:
#+begin_src matlab
@@ -887,7 +982,7 @@ In order to visualize the effect of $\omega_i$ on the attainable damping, the Ro
for wi_i = 1:length(wis)
set(gca,'ColorOrderIndex',wi_i);
wi = wis(wi_i);
- L(wi_i) = plot(nan, nan, '.', 'DisplayName', sprintf('$\\omega_i = %.2f \\omega_0$', wi./w0));
+ L(wi_i) = plot(nan, nan, '.', 'DisplayName', sprintf('$\\Omega_i = %.2f \\omega_0$', wi./w0));
for g = gains
clpoles = pole(feedback(Giff, (g/(wi+s))*eye(2)));
set(gca,'ColorOrderIndex',wi_i);
@@ -934,6 +1029,64 @@ In order to visualize the effect of $\omega_i$ on the attainable damping, the Ro
#+RESULTS:
[[file:figs/root_locus_wi_modified_iff.png]]
+#+begin_src matlab :exports none
+ gains = logspace(-2, 4, 100);
+
+ poles_iff_hpf = rootLocusPolesSorted(Giff, 1/(s + wi)*eye(2), gains, 'd_max', 1e-4);
+
+ figure;
+
+ ax1 = subplot(1, 2, 1);
+ hold on;
+ for wi_i = 1:length(wis)
+ wi = wis(wi_i);
+
+ set(gca,'ColorOrderIndex',wi_i);
+ L(wi_i) = plot(nan, nan, '.', 'DisplayName', sprintf('$\\Omega_i = %.2f \\omega_0$', wi./w0));
+
+ poles = rootLocusPolesSorted(Giff, 1/(s + wi)*eye(2), gains, 'd_max', 1e-4);
+ for p_i = 1:size(poles, 2)
+ set(gca,'ColorOrderIndex',wi_i);
+ plot(real(poles(:, p_i)), imag(poles(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+ end
+ plot(real(pole(Giff)), imag(pole(Giff)), 'kx');
+ plot(real(tzero(Giff)), imag(tzero(Giff)), 'ko');
+ hold off;
+ axis square;
+ xlim([-2.3, 0.1]); ylim([-1.2, 1.2]);
+ xticks([-2:1:2]); yticks([-2:1:2]);
+ legend(L, 'location', 'northwest');
+ xlabel('Real Part'); ylabel('Imaginary Part');
+
+ clear L
+
+ ax2 = subplot(1, 2, 2);
+ hold on;
+ for wi_i = 1:length(wis)
+ wi = wis(wi_i);
+
+ poles = rootLocusPolesSorted(Giff, 1/(s + wi)*eye(2), gains, 'd_max', 1e-4);
+ for p_i = 1:size(poles, 2)
+ set(gca,'ColorOrderIndex',wi_i);
+ plot(real(poles(:, p_i)), imag(poles(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+ end
+ plot(real(pole(Giff)), imag(pole(Giff)), 'kx');
+ plot(real(tzero(Giff)), imag(tzero(Giff)), 'ko');
+ hold off;
+ axis square;
+ xlim([-0.2, 0.1]); ylim([-0.15, 0.15]);
+ xticks([-0.2:0.1:0.1]); yticks([-0.2:0.1:0.2]);
+ xlabel('Real Part'); ylabel('Imaginary Part');
+#+end_src
+
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/root_locus_wi_modified_iff.pdf', 'width', 'full', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
For the controller
\begin{equation}
K_{\text{IFF}}(s) = g\frac{1}{\omega_i + s} \begin{bmatrix}
@@ -942,9 +1095,8 @@ For the controller
\end{bmatrix}
\end{equation}
The gain at which the system becomes unstable is
-#+name: eq:iff_gmax
\begin{equation}
- g_\text{max} = \omega_i \left( \frac{{\omega_0}^2}{\Omega_2} - 1 \right)
+ 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 [[fig:root_locus_wi_modified_iff]]), the control gains may be limited to small values due to eqref:eq:iff_gmax thus reducing the attainable damping.
@@ -981,7 +1133,7 @@ To find the optimum, the gain that maximize the simultaneous damping of the mode
plot(wis, opt_zeta, '-o', 'DisplayName', '$\xi_{cl}$');
set(gca, 'YScale', 'lin');
ylim([0,1]);
- ylabel('Attainable Damping Ration $\xi$');
+ ylabel('Attainable Damping Ratio $\xi$');
yyaxis right
hold on;
@@ -1005,6 +1157,9 @@ To find the optimum, the gain that maximize the simultaneous damping of the mode
#+RESULTS:
[[file:figs/mod_iff_damping_wi.png]]
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/root_locus_wi_modified_iff.pdf', 'width', 'wide', 'height', 'normal', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
* IFF with a stiffness in parallel with the force sensor
<>
@@ -1021,11 +1176,69 @@ To find the optimum, the gain that maximize the simultaneous damping of the mode
#+begin_src matlab
addpath('./matlab/');
+ addpath('./src/');
#+end_src
-#+begin_src matlab
- open('rotating_frame.slx');
-#+end_src
+** Schematic
+
+#+name: fig:rotating_xy_platform_springs
+#+caption: Figure caption
+[[file:figs-tikz/rotating_xy_platform_springs.png]]
+
+** Equations
+#+begin_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^\prime}^2} + 2\xi^\prime \frac{s}{{\omega_0^\prime}^2} + 1 - \frac{\Omega^2}{{\omega_0^\prime}^2} \right)^2 + \left( 2 \frac{\Omega}{\omega_0^\prime}\frac{s}{\omega_0^\prime} \right)^2 \\
+ G_{kz} &= \left( \frac{s^2}{{\omega_0^\prime}^2} + \frac{k_p}{k + k_p} - \frac{\Omega^2}{{\omega_0^\prime}^2} \right) \left( \frac{s^2}{{\omega_0^\prime}^2} + 2\xi^\prime \frac{s}{{\omega_0^\prime}^2} + 1 - \frac{\Omega^2}{{\omega_0^\prime}^2} \right) + \left( 2 \frac{\Omega}{\omega_0^\prime}\frac{s}{\omega_0^\prime} \right)^2 \\
+ G_{kc} &= \left( 2 \xi^\prime \frac{s}{\omega_0^\prime} + \frac{k}{k + k_p} \right) \left( 2 \frac{\Omega}{\omega_0^\prime}\frac{s}{\omega_0^\prime} \right)
+\end{align}
+where:
+- $\omega_0^\prime = \frac{k + k_p}{m}$
+- $\xi^\prime = \frac{c}{2 \sqrt{(k + k_p) m}}$
+#+end_important
+
+If we compare $G_{kz}$ and $G_{fz}$, we see that the spring in parallel adds a term $\frac{k_p}{k + k_p}$.
+In order to have two complex conjugate zeros (instead of real zeros):
+\begin{equation}
+ \frac{k_p}{k + k_p} - \frac{\Omega^2}{{\omega_0^\prime}^2} > 0
+\end{equation}
+Which is equivalent to
+\begin{equation}
+ k_p > m \Omega^2
+\end{equation}
+
+** 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}
** Plant Parameters
Let's define initial values for the model.
@@ -1045,39 +1258,13 @@ Let's define initial values for the model.
cp = 0; % [N/(m/s)]
#+end_src
-** Schematic
-
-#+name: fig:rotating_xy_platform_springs
-#+caption: Figure caption
-[[file:figs-tikz/rotating_xy_platform_springs.png]]
-
-** 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}
-
-** 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
-
-** Effect of the parallel stiffness on the IFF plant
-The rotation speed is set to $\Omega = 0.1 \omega_0$.
+** 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.
#+begin_src matlab
W = 0.1*w0; % [rad/s]
+
+ kp = 1.5*m*W^2;
+ cp = 0;
#+end_src
#+begin_src matlab :exports none
@@ -1085,6 +1272,102 @@ The rotation speed is set to $\Omega = 0.1 \omega_0$.
Kdvf = tf(zeros(2));
#+end_src
+#+begin_src matlab
+ open('rotating_frame.slx');
+#+end_src
+
+#+begin_src 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;
+
+ Giff = linearize(mdl, io, 0);
+
+ %% Input/Output definition
+ Giff.InputName = {'Fu', 'Fv'};
+ Giff.OutputName = {'fu', 'fv'};
+#+end_src
+
+#+begin_src 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'};
+#+end_src
+
+#+begin_src matlab :exports none
+ freqs = logspace(-1, 1, 1000);
+
+ figure;
+ ax1 = subplot(2, 2, 1);
+ hold on;
+ plot(freqs, abs(squeeze(freqresp(Giff(1,1), freqs))), '-')
+ plot(freqs, abs(squeeze(freqresp(Giff_th(1,1), freqs))), '--')
+ hold off;
+ set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ set(gca, 'XTickLabel',[]); ylabel('Magnitude [N/N]');
+ title('$f_u/F_u$, $f_v/F_v$');
+
+ ax3 = subplot(2, 2, 3);
+ hold on;
+ plot(freqs, 180/pi*angle(squeeze(freqresp(Giff(1,1), freqs))), '-')
+ plot(freqs, 180/pi*angle(squeeze(freqresp(Giff_th(1,1), freqs))), '--')
+ set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+ xlabel('Frequency [rad/s]'); ylabel('Phase [deg]');
+ yticks(-180:90:180);
+ ylim([-180 180]);
+ hold off;
+
+ ax2 = subplot(2, 2, 2);
+ hold on;
+ plot(freqs, abs(squeeze(freqresp(Giff(1,2), freqs))), '-')
+ plot(freqs, abs(squeeze(freqresp(Giff_th(1,2), freqs))), '--')
+ hold off;
+ set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ set(gca, 'XTickLabel',[]); ylabel('Magnitude [N/N]');
+ title('$f_u/F_v$, $f_v/F_u$');
+
+ ax4 = subplot(2, 2, 4);
+ hold on;
+ plot(freqs, 180/pi*angle(squeeze(freqresp(Giff(1,2), freqs))), '-', ...
+ 'DisplayName', 'Simscape')
+ plot(freqs, 180/pi*angle(squeeze(freqresp(Giff_th(1,2), freqs))), '--', ...
+ 'DisplayName', 'Analytical')
+ set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+ xlabel('Frequency [rad/s]'); ylabel('Phase [deg]');
+ yticks(-180:90:180);
+ ylim([-180 180]);
+ hold off;
+ legend('location', 'northeast');
+
+ linkaxes([ax1,ax2,ax3,ax4],'x');
+ xlim([freqs(1), freqs(end)]);
+ linkaxes([ax1,ax2],'y');
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+ exportFig('figs/plant_iff_kp_comp_simscape_analytical.pdf', 'width', 'full', 'height', 'full');
+#+end_src
+
+#+name: fig:plant_iff_kp_comp_simscape_analytical
+#+caption: Comparison of the transfer functions from $[F_u, F_v]$ to $[f_u, f_v]$ between the Simscape model and the analytical one
+#+RESULTS:
+[[file:figs/plant_iff_kp_comp_simscape_analytical.png]]
+
+** Effect of the parallel stiffness on the IFF plant
+The rotation speed is set to $\Omega = 0.1 \omega_0$.
+#+begin_src matlab
+ W = 0.1*w0; % [rad/s]
+#+end_src
+
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$
@@ -1094,47 +1377,40 @@ 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 none
- %% 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;
-#+end_src
-
#+begin_src matlab
kp = 0;
cp = 0;
- Giff = linearize(mdl, io, 0);
+ w0p = sqrt((k + kp)/m);
+ xip = c/(2*sqrt((k+kp)*m));
- %% Input/Output definition
- Giff.InputName = {'Fu', 'Fv'};
- Giff.OutputName = {'fu', 'fv'};
+ 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];
#+end_src
#+begin_src matlab
kp = 0.5*m*W^2;
- cp = 0.001;
+ cp = 0;
- Giff_s = linearize(mdl, io, 0);
+ w0p = sqrt((k + kp)/m);
+ xip = c/(2*sqrt((k+kp)*m));
- %% Input/Output definition
- Giff_s.InputName = {'Fu', 'Fv'};
- Giff_s.OutputName = {'fu', 'fv'};
+ 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];
#+end_src
#+begin_src matlab
kp = 1.5*m*W^2;
- cp = 0.001;
+ cp = 0;
- Giff_l = linearize(mdl, io, 0);
+ w0p = sqrt((k + kp)/m);
+ xip = c/(2*sqrt((k+kp)*m));
- %% Input/Output definition
- Giff_l.InputName = {'Fu', 'Fv'};
- Giff_l.OutputName = {'fu', 'fv'};
+ 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];
#+end_src
#+begin_src matlab :exports none
@@ -1178,6 +1454,11 @@ One can see that for $k_p > m \Omega^2$, the systems shows alternating complex c
#+caption: 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$
#+RESULTS:
[[file:figs/plant_iff_kp.png]]
+
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/plant_iff_kp.pdf', 'width', 'full', 'height', 'full', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
** IFF when adding a spring in parallel
In Figure [[fig:root_locus_iff_kp]] is displayed the Root Locus in the three considered cases with
\begin{equation}
@@ -1293,6 +1574,100 @@ Thus, decentralized IFF controller with pure integrators can be used if:
#+RESULTS:
[[file:figs/root_locus_iff_kp.png]]
+#+begin_src matlab :exports none :tangle no
+ gains = logspace(-2, 2, 200);
+
+ poles_iff = rootLocusPolesSorted(Giff, 1/s*eye(2), gains, 'd_max', 1e-4);
+ poles_iff_s = rootLocusPolesSorted(Giff_s, 1/s*eye(2), gains, 'd_max', 1e-4);
+ poles_iff_l = rootLocusPolesSorted(Giff_l, 1/s*eye(2), gains, 'd_max', 1e-4);
+
+ figure;
+
+ subplot(1,2,1);
+ hold on;
+ set(gca,'ColorOrderIndex',1);
+ plot(real(pole(Giff)), imag(pole(Giff)), 'x', ...
+ 'DisplayName', '$k_p = 0$');
+ set(gca,'ColorOrderIndex',1);
+ plot(real(tzero(Giff)), imag(tzero(Giff)), 'o', ...
+ 'HandleVisibility', 'off');
+ for p_i = 1:size(poles_iff, 2)
+ set(gca,'ColorOrderIndex',1);
+ plot(real(poles_iff(:, p_i)), imag(poles_iff(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+
+ set(gca,'ColorOrderIndex',2);
+ plot(real(pole(Giff_s)), imag(pole(Giff_s)), 'x', ...
+ 'DisplayName', '$k_p < m\Omega^2$');
+ set(gca,'ColorOrderIndex',2);
+ plot(real(tzero(Giff_s)), imag(tzero(Giff_s)), 'o', ...
+ 'HandleVisibility', 'off');
+ for p_i = 1:size(poles_iff_s, 2)
+ set(gca,'ColorOrderIndex',2);
+ plot(real(poles_iff_s(:, p_i)), imag(poles_iff_s(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+
+ set(gca,'ColorOrderIndex',3);
+ plot(real(pole(Giff_l)), imag(pole(Giff_l)), 'x', ...
+ 'DisplayName', '$k_p > m\Omega^2$');
+ set(gca,'ColorOrderIndex',3);
+ plot(real(tzero(Giff_l)), imag(tzero(Giff_l)), 'o', ...
+ 'HandleVisibility', 'off');
+ for p_i = 1:size(poles_iff_l, 2)
+ set(gca,'ColorOrderIndex',3);
+ plot(real(poles_iff_l(:, p_i)), imag(poles_iff_l(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+ hold off;
+ axis square;
+ xlim([-1, 0.2]); ylim([0, 1.2]);
+
+ xlabel('Real Part'); ylabel('Imaginary Part');
+ legend('location', 'northwest');
+
+ subplot(1,2,2);
+ hold on;
+ set(gca,'ColorOrderIndex',1);
+ plot(real(pole(Giff)), imag(pole(Giff)), 'x');
+ set(gca,'ColorOrderIndex',1);
+ plot(real(tzero(Giff)), imag(tzero(Giff)), 'o');
+ for p_i = 1:size(poles_iff, 2)
+ set(gca,'ColorOrderIndex',1);
+ plot(real(poles_iff(:, p_i)), imag(poles_iff(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+
+ set(gca,'ColorOrderIndex',2);
+ plot(real(pole(Giff_s)), imag(pole(Giff_s)), 'x');
+ set(gca,'ColorOrderIndex',2);
+ plot(real(tzero(Giff_s)), imag(tzero(Giff_s)), 'o');
+ for p_i = 1:size(poles_iff_s, 2)
+ set(gca,'ColorOrderIndex',2);
+ plot(real(poles_iff_s(:, p_i)), imag(poles_iff_s(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+
+ set(gca,'ColorOrderIndex',3);
+ plot(real(pole(Giff_l)), imag(pole(Giff_l)), 'x');
+ set(gca,'ColorOrderIndex',3);
+ plot(real(tzero(Giff_l)), imag(tzero(Giff_l)), 'o');
+ for p_i = 1:size(poles_iff_l, 2)
+ set(gca,'ColorOrderIndex',3);
+ plot(real(poles_iff_l(:, p_i)), imag(poles_iff_l(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+ hold off;
+ axis square;
+ xlim([-0.04, 0.06]); ylim([0, 0.1]);
+ xlabel('Real Part'); ylabel('Imaginary Part');
+#+end_src
+
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/root_locus_iff_kp.pdf', 'width', 'full', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
** Effect of $k_p$ on the attainable damping
However, having large values of $k_p$ may:
- decrease the actuator force authority
@@ -1300,8 +1675,7 @@ However, having large values of $k_p$ may:
To study the second point, Root Locus plots for the following values of $k_p$ are shown in Figure [[fig:root_locus_iff_kps]].
#+begin_src matlab
- kps = [1, 5, 10, 50]*m*W^2;
- cp = 0.01;
+ kps = [2, 20, 40]*m*W^2;
#+end_src
It is shown that large values of $k_p$ decreases the attainable damping.
@@ -1309,12 +1683,18 @@ It is shown that large values of $k_p$ decreases the attainable damping.
#+begin_src matlab :exports none
figure;
- gains = logspace(-2, 4, 100);
+ gains = logspace(-2, 4, 500);
hold on;
for kp_i = 1:length(kps)
kp = kps(kp_i);
- Giff = linearize(mdl, io, 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 ];
set(gca,'ColorOrderIndex',kp_i);
plot(real(pole(Giff)), imag(pole(Giff)), 'x', ...
@@ -1343,25 +1723,73 @@ It is shown that large values of $k_p$ decreases the attainable damping.
#+end_src
#+name: fig:root_locus_iff_kps
-#+caption:
+#+caption: Root Locus plot
#+RESULTS:
[[file:figs/root_locus_iff_kps.png]]
+#+begin_src matlab :exports none :tangle no
+ gains = logspace(-2, 4, 500);
+
+ figure;
+
+ hold on;
+ for kp_i = 1:length(kps)
+ kp = kps(kp_i);
+ k = 1 - kp;
+
+ 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 ];
+
+ poles_iff = rootLocusPolesSorted(Giff, 1/s*eye(2), gains, 'd_max', 1e-4);
+
+ set(gca,'ColorOrderIndex',kp_i);
+ plot(real(pole(Giff)), imag(pole(Giff)), 'x', ...
+ 'DisplayName', sprintf('$k_p = %.1f m \\Omega^2$', kp/(m*W^2)));
+ set(gca,'ColorOrderIndex',kp_i);
+ plot(real(tzero(Giff)), imag(tzero(Giff)), 'o', ...
+ 'HandleVisibility', 'off');
+ for p_i = 1:size(poles_iff, 2)
+ set(gca,'ColorOrderIndex', kp_i);
+ plot(real(poles_iff(:, p_i)), imag(poles_iff(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+ end
+ hold off;
+ axis square;
+ xlim([-1.2, 0.2]); ylim([0, 1.4]);
+
+ xlabel('Real Part'); ylabel('Imaginary Part');
+ legend('location', 'northwest');
+#+end_src
+
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/root_locus_iff_kps.pdf', 'width', 'wide', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
** 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.
#+begin_src matlab
kp = 5*m*W^2;
cp = 0.01;
-
- Giff = linearize(mdl, io, 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 ];
#+end_src
#+begin_src matlab
opt_zeta = 0;
opt_gain = 0;
- gains = logspace(-2, 4, 100);
+ gains = logspace(-2, 4, 1000);
for g = gains
Kiff = (g/s)*eye(2);
@@ -1378,17 +1806,17 @@ Let's take $k_p = 5 m \Omega^2$ and find the optimal IFF control gain $g$ such t
#+begin_src matlab :exports none
figure;
- gains = logspace(-2, 4, 100);
+ gains = logspace(-2, 4, 1000);
hold on;
plot(real(pole(Giff)), imag(pole(Giff)), 'kx');
plot(real(tzero(Giff)), imag(tzero(Giff)), 'ko');
for g = gains
- clpoles = pole(feedback(Giff, (g/s)*eye(2)));
+ clpoles = pole(minreal(feedback(Giff, (g/s)*eye(2))));
plot(real(clpoles), imag(clpoles), 'k.');
end
% Optimal Gain
- clpoles = pole(feedback(Giff, (opt_gain/s)*eye(2)));
+ clpoles = pole(minreal(feedback(Giff, (opt_gain/s)*eye(2))));
set(gca,'ColorOrderIndex',1);
plot(real(clpoles), imag(clpoles), 'x');
for clpole = clpoles'
@@ -1406,10 +1834,43 @@ Let's take $k_p = 5 m \Omega^2$ and find the optimal IFF control gain $g$ such t
#+end_src
#+name: fig:root_locus_opt_gain_iff_kp
-#+caption:
+#+caption: Root Locus for $k_p = 5 m \Omega^2$ and the poles corresponding to the identified optimal gain
#+RESULTS:
[[file:figs/root_locus_opt_gain_iff_kp.png]]
+#+begin_src matlab :exports none :tangle no
+ gains = logspace(-2, 4, 1000);
+
+ poles_iff = rootLocusPolesSorted(Giff, 1/s*eye(2), gains, 'd_max', 1e-4);
+
+ figure;
+
+ hold on;
+ plot(real(pole(Giff)), imag(pole(Giff)), 'kx');
+ plot(real(tzero(Giff)), imag(tzero(Giff)), 'ko');
+ for p_i = 1:size(poles_iff, 2)
+ plot(real(poles_iff(:, p_i)), imag(poles_iff(:, p_i)), 'k-', ...
+ 'HandleVisibility', 'off');
+ end
+ % Optimal Gain
+ clpoles = pole(minreal(feedback(Giff, (opt_gain/s)*eye(2))));
+ set(gca,'ColorOrderIndex',1);
+ plot(real(clpoles), imag(clpoles), 'x');
+ for clpole = clpoles'
+ set(gca,'ColorOrderIndex',1);
+ plot([0, real(clpole)], [0, imag(clpole)], '-', 'LineWidth', 1);
+ end
+ hold off;
+ axis square;
+ xlim([-1.2, 0.2]); ylim([0, 1.4]);
+
+ xlabel('Real Part'); ylabel('Imaginary Part');
+#+end_src
+
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/root_locus_opt_gain_iff_kp.pdf', 'width', 'wide', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
* Direct Velocity Feedback
<>
@@ -1425,10 +1886,7 @@ Let's take $k_p = 5 m \Omega^2$ and find the optimal IFF control gain $g$ such t
#+begin_src matlab
addpath('./matlab/');
-#+end_src
-
-#+begin_src matlab
- open('rotating_frame.slx');
+ addpath('./src/');
#+end_src
** Equations
@@ -1436,13 +1894,25 @@ The sensed relative velocity are equal to:
#+begin_important
\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}
+\bm{G}_v
+\begin{bmatrix} F_u \\ F_v \end{bmatrix}
+\end{equation}
+
+\begin{equation}
+\begin{bmatrix} \dot{d}_u \\ \dot{d}_v \end{bmatrix} =
+\frac{s}{k} \frac{1}{G_{vp}}
\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} \\
+ G_{vz} & G_{vc} \\
+ -G_{vc} & G_{vz}
\end{bmatrix}
\begin{bmatrix} F_u \\ F_v \end{bmatrix}
\end{equation}
+With:
+\begin{align}
+ G_{vp} &= \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_{vz} &= \frac{s^2}{{\omega_0}^2} + 2 \xi \frac{s}{\omega_0} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} \\
+ G_{vc} &= 2 \frac{\Omega}{\omega_0} \frac{s}{\omega_0}
+\end{align}
#+end_important
** Plant Parameters
@@ -1463,7 +1933,7 @@ Let's define initial values for the model.
cp = 0; % [N/(m/s)]
#+end_src
-** Plant - Bode Plot
+** Comparison of the Analytical Model and the Simscape Model
The rotating speed is set to $\Omega = 0.1 \omega_0$.
#+begin_src matlab
W = 0.1*w0;
@@ -1474,6 +1944,10 @@ The rotating speed is set to $\Omega = 0.1 \omega_0$.
Kdvf = tf(zeros(2));
#+end_src
+#+begin_src matlab
+ open('rotating_frame.slx');
+#+end_src
+
And the transfer function from $[F_u, F_v]$ to $[v_u, v_v]$ is identified using the Simscape model.
#+begin_src matlab
%% Name of the Simulink File
@@ -1493,7 +1967,6 @@ And the transfer function from $[F_u, F_v]$ to $[v_u, v_v]$ is identified using
Gdvf.OutputName = {'Vu', 'Vv'};
#+end_src
-** 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.
#+begin_src matlab
Gdvf_th = (s/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ...
@@ -1574,11 +2047,10 @@ The Decentralized Direct Velocity Feedback controller consist of a pure gain on
The corresponding Root Locus plots for the following rotating speeds are shown in Figure [[fig:root_locus_dvf]].
#+begin_src matlab
- Ws = [0, 0.1, 0.5, 0.8, 1.1]*w0; % Rotating Speeds [rad/s]
+ Ws = [0, 0.2, 0.7, 1.1]*w0; % Rotating Speeds [rad/s]
#+end_src
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 none
gains = logspace(-2, 1, 100);
@@ -1593,7 +2065,7 @@ It is shown that for rotating speed $\Omega < \omega_0$, the closed loop system
set(gca,'ColorOrderIndex',W_i);
plot(real(pole(Gdvf)), imag(pole(Gdvf)), 'x', ...
- 'DisplayName', sprintf('$\\omega = %.2f \\omega_0 $', W/w0));
+ 'DisplayName', sprintf('$\\Omega = %.2f \\omega_0 $', W/w0));
set(gca,'ColorOrderIndex',W_i);
plot(real(tzero(Gdvf)), imag(tzero(Gdvf)), 'o', ...
@@ -1624,6 +2096,45 @@ It is shown that for rotating speed $\Omega < \omega_0$, the closed loop system
#+RESULTS:
[[file:figs/root_locus_dvf.png]]
+#+begin_src matlab :exports none
+ gains = logspace(-2, 1, 1000);
+
+ figure;
+ hold on;
+ for W_i = 1:length(Ws)
+ W = Ws(W_i);
+
+ Gdvf = (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)];
+
+ set(gca,'ColorOrderIndex',W_i);
+ plot(real(pole(Gdvf)), imag(pole(Gdvf)), 'x', ...
+ 'DisplayName', sprintf('$\\Omega = %.2f \\omega_0 $', W/w0));
+
+ set(gca,'ColorOrderIndex',W_i);
+ plot(real(tzero(Gdvf)), imag(tzero(Gdvf)), 'o', ...
+ 'HandleVisibility', 'off');
+
+ poles_dvf = rootLocusPolesSorted(Gdvf, eye(2), gains, 'd_max', 1e-4);
+ for p_i = 1:size(poles_dvf, 2)
+ set(gca,'ColorOrderIndex', W_i);
+ plot(real(poles_dvf(:, p_i)), imag(poles_dvf(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+ end
+ hold off;
+ axis square;
+ xlim([-2, 0.5]); ylim([0, 2.5]);
+
+ xlabel('Real Part'); ylabel('Imaginary Part');
+ legend('location', 'northwest');
+#+end_src
+
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/root_locus_dvf.pdf', 'width', 'wide', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
* Comparison
<>
@@ -1640,10 +2151,7 @@ It is shown that for rotating speed $\Omega < \omega_0$, the closed loop system
#+begin_src matlab
addpath('./matlab/');
-#+end_src
-
-#+begin_src matlab
- open('rotating_frame.slx');
+ addpath('./src/');
#+end_src
** Plant Parameters
@@ -1664,94 +2172,43 @@ Let's define initial values for the model.
cp = 0; % [N/(m/s)]
#+end_src
-#+begin_src matlab :exports none
- Kiff = tf(zeros(2));
- Kdvf = tf(zeros(2));
-#+end_src
-
The rotating speed is set to $\Omega = 0.1 \omega_0$.
#+begin_src matlab
W = 0.1*w0;
#+end_src
** Root Locus
-*** Pseudo Integrator IFF :ignore:
-#+begin_src matlab :exports none
- kp = 0;
- cp = 0;
-#+end_src
-
+IFF with High Pass Filter
#+begin_src matlab
- wi = 0.1*w0;
+ 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];
#+end_src
-#+begin_src 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;
-#+end_src
-
-#+begin_src matlab
- Giff = linearize(mdl, io, 0);
-
- %% Input/Output definition
- Giff.InputName = {'Fu', 'Fv'};
- Giff.OutputName = {'Fmu', 'Fmv'};
-#+end_src
-
-*** IFF With parallel Stiffness :ignore:
+IFF With parallel Stiffness
#+begin_src matlab
kp = 5*m*W^2;
- cp = 0.01;
+ 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;
#+end_src
+DVF
#+begin_src 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;
+ Gdvf = (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)];
#+end_src
-#+begin_src matlab
- Giff_kp = linearize(mdl, io, 0);
-
- %% Input/Output definition
- Giff_kp.InputName = {'Fu', 'Fv'};
- Giff_kp.OutputName = {'Fmu', 'Fmv'};
-#+end_src
-
-*** DVF :ignore:
-#+begin_src matlab :exports none
- kp = 0;
- cp = 0;
-#+end_src
-
-#+begin_src 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;
-#+end_src
-
-#+begin_src matlab
- Gdvf = linearize(mdl, io, 0);
-
- %% Input/Output definition
- Gdvf.InputName = {'Fu', 'Fv'};
- Gdvf.OutputName = {'Vu', 'Vv'};
-#+end_src
-
-*** Root Locus :ignore:
#+begin_src matlab :exports none
figure;
@@ -1760,7 +2217,7 @@ The rotating speed is set to $\Omega = 0.1 \omega_0$.
hold on;
set(gca,'ColorOrderIndex',1);
plot(real(pole(Giff)), imag(pole(Giff)), 'x', ...
- 'DisplayName', 'Pseudo Integrator');
+ 'DisplayName', 'IFF + LFP');
set(gca,'ColorOrderIndex',1);
plot(real(tzero(Giff)), imag(tzero(Giff)), 'o', ...
'HandleVisibility', 'off');
@@ -1774,7 +2231,7 @@ The rotating speed is set to $\Omega = 0.1 \omega_0$.
set(gca,'ColorOrderIndex',2);
plot(real(pole(Giff_kp)), imag(pole(Giff_kp)), 'x', ...
- 'DisplayName', 'Parallel Stiffness');
+ 'DisplayName', 'IFF + $k_p$');
set(gca,'ColorOrderIndex',2);
plot(real(tzero(Giff_kp)), imag(tzero(Giff_kp)), 'o', ...
'HandleVisibility', 'off');
@@ -1812,10 +2269,67 @@ The rotating speed is set to $\Omega = 0.1 \omega_0$.
#+end_src
#+name: fig:comp_root_locus
-#+caption:
+#+caption: Root Locus plot - Comparison of IFF with additional high pass filter, IFF with additional parallel stiffness and DVF
#+RESULTS:
[[file:figs/comp_root_locus.png]]
+#+begin_src matlab :exports none :tangle no
+ gains = logspace(-2, 2, 1000);
+
+ poles_iff_hpf = rootLocusPolesSorted(Giff, 1/(s + wi)*eye(2), gains, 'd_max', 1e-4);
+ poles_iff_kp = rootLocusPolesSorted(Giff_kp, 1/s*eye(2), gains, 'd_max', 1e-4);
+ poles_dvf = rootLocusPolesSorted(Gdvf, eye(2), gains, 'd_max', 1e-4);
+
+ figure;
+
+ hold on;
+ set(gca,'ColorOrderIndex',1);
+ plot(real(pole(Giff)), imag(pole(Giff)), 'x', ...
+ 'DisplayName', 'IFF + LFP');
+ set(gca,'ColorOrderIndex',1);
+ plot(real(tzero(Giff)), imag(tzero(Giff)), 'o', ...
+ 'HandleVisibility', 'off');
+ for p_i = 1:size(poles_iff_hpf, 2)
+ set(gca,'ColorOrderIndex',1);
+ plot(real(poles_iff_hpf(:, p_i)), imag(poles_iff_hpf(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+
+ set(gca,'ColorOrderIndex',2);
+ plot(real(pole(Giff_kp)), imag(pole(Giff_kp)), 'x', ...
+ 'DisplayName', 'IFF + $k_p$');
+ set(gca,'ColorOrderIndex',2);
+ plot(real(tzero(Giff_kp)), imag(tzero(Giff_kp)), 'o', ...
+ 'HandleVisibility', 'off');
+ for p_i = 1:size(poles_iff_kp, 2)
+ set(gca,'ColorOrderIndex',2);
+ plot(real(poles_iff_kp(:, p_i)), imag(poles_iff_kp(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+
+ set(gca,'ColorOrderIndex',3);
+ plot(real(pole(Gdvf)), imag(pole(Gdvf)), 'x', ...
+ 'DisplayName', 'DVF');
+ set(gca,'ColorOrderIndex',3);
+ plot(real(tzero(Gdvf)), imag(tzero(Gdvf)), 'o', ...
+ 'HandleVisibility', 'off');
+ for p_i = 1:size(poles_dvf, 2)
+ set(gca,'ColorOrderIndex',3);
+ plot(real(poles_dvf(:, p_i)), imag(poles_dvf(:, p_i)), '-', ...
+ 'HandleVisibility', 'off');
+ end
+ hold off;
+ axis square;
+ xlim([-1.2, 0.05]); ylim([0, 1.25]);
+
+ xlabel('Real Part'); ylabel('Imaginary Part');
+ legend('location', 'northwest');
+#+end_src
+
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/comp_root_locus.pdf', 'width', 'wide', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
** 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.
@@ -1883,11 +2397,16 @@ The obtained damping ratio and control are shown below.
| | Obtained $\xi$ | Control Gain |
|----------------+----------------+--------------|
| Modified IFF | 0.83 | 2.0 |
-| IFF with $k_p$ | 0.84 | 2.01 |
+| IFF with $k_p$ | 0.83 | 2.01 |
| DVF | 0.85 | 1.67 |
** Transmissibility
<>
+
+#+begin_src matlab
+ open('rotating_frame.slx');
+#+end_src
+
*** Open Loop :ignore:
#+begin_src matlab :exports none
Kdvf = tf(zeros(2));
@@ -2014,17 +2533,17 @@ The obtained damping ratio and control are shown below.
figure;
hold on;
plot(freqs, abs(squeeze(freqresp(Tiff(1,1), freqs))), ...
- 'DisplayName', 'IFF Pseudo int')
+ 'DisplayName', 'IFF + HPF')
plot(freqs, abs(squeeze(freqresp(Tiff_kp(1,1), freqs))), ...
- 'DisplayName', 'IFF Paral. stiff')
+ 'DisplayName', 'IFF + $k_p$')
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')
+ 'DisplayName', 'Open-Loop')
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [rad/s]'); ylabel('Transmissibility [m/m]');
- legend('location', 'northwest');
+ legend('location', 'southwest');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
@@ -2032,10 +2551,14 @@ The obtained damping ratio and control are shown below.
#+end_src
#+name: fig:comp_transmissibility
-#+caption:
+#+caption: Comparison of the transmissibility
#+RESULTS:
[[file:figs/comp_transmissibility.png]]
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/comp_transmissibility.pdf', 'width', 'wide', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
** Compliance
<>
*** Open Loop :ignore:
@@ -2134,13 +2657,13 @@ The obtained damping ratio and control are shown below.
figure;
hold on;
plot(freqs, abs(squeeze(freqresp(Ciff(1,1), freqs))), ...
- 'DisplayName', 'IFF Pseudo int')
+ 'DisplayName', 'IFF + LPF')
plot(freqs, abs(squeeze(freqresp(Ciff_kp(1,1), freqs))), ...
- 'DisplayName', 'IFF Paral. stiff')
+ 'DisplayName', 'IFF + $k_p$')
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')
+ 'DisplayName', 'Open-Loop')
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [rad/s]'); ylabel('Compliance [m/N]');
@@ -2156,6 +2679,10 @@ The obtained damping ratio and control are shown below.
#+RESULTS:
[[file:figs/comp_compliance.png]]
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/comp_compliance.pdf', 'width', 'wide', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
* Notations
<>
@@ -2185,3 +2712,150 @@ The obtained damping ratio and control are shown below.
| 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= | |
+
+* Function
+** Sort Poles for the Root Locus
+:PROPERTIES:
+:header-args:matlab+: :tangle src/rootLocusPolesSorted.m
+:header-args:matlab+: :comments none :mkdirp yes :eval no
+:END:
+<>
+
+This Matlab function is accessible [[file:src/rootLocusPolesSorted.m][here]].
+
+*** Function description
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+#+begin_src matlab
+ function [poles] = rootLocusPolesSorted(G, K, gains, args)
+ % rootLocusPolesSorted -
+ %
+ % Syntax: [poles] = rootLocusPolesSorted(G, K, gains, args)
+ %
+ % Inputs:
+ % - G, K, gains, args -
+ %
+ % Outputs:
+ % - poles -
+#+end_src
+
+*** Optional Parameters
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+#+begin_src matlab
+ arguments
+ G
+ K
+ gains
+ args.minreal double {mustBeNumericOrLogical} = false
+ args.p_half double {mustBeNumericOrLogical} = false
+ args.d_max double {mustBeNumeric} = -1
+ end
+#+end_src
+
+*** Function
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+#+begin_src matlab
+ if args.minreal
+ p1 = pole(minreal(feedback(G, gains(1)*K)));
+ [~, i_uniq] = uniquetol([real(p1), imag(p1)], 1e-10, 'ByRows', true);
+ p1 = p1(i_uniq);
+
+ poles = zeros(length(p1), length(gains));
+ poles(:, 1) = p1;
+ else
+ p1 = pole(feedback(G, gains(1)*K));
+ [~, i_uniq] = uniquetol([real(p1), imag(p1)], 1e-10, 'ByRows', true);
+ p1 = p1(i_uniq);
+
+ poles = zeros(length(p1), length(gains));
+ poles(:, 1) = p1;
+ end
+#+end_src
+
+#+begin_src matlab
+ if args.minreal
+ p2 = pole(minreal(feedback(G, gains(2)*K)));
+ [~, i_uniq] = uniquetol([real(p2), imag(p2)], 1e-10, 'ByRows', true);
+ p2 = p2(i_uniq);
+ poles(:, 2) = p2;
+ else
+ p2 = pole(feedback(G, gains(2)*K));
+ [~, i_uniq] = uniquetol([real(p2), imag(p2)], 1e-10, 'ByRows', true);
+ p2 = p2(i_uniq);
+ poles(:, 2) = p2;
+ end
+#+end_src
+
+#+begin_src matlab
+ for g_i = 3:length(gains)
+ % Estimated value of the poles
+ poles_est = poles(:, g_i-1) + (poles(:, g_i-1) - poles(:, g_i-2))*(gains(g_i) - gains(g_i-1))/(gains(g_i-1) - gains(g_i - 2));
+
+ % New values for the poles
+ poles_gi = pole(feedback(G, gains(g_i)*K));
+ [~, i_uniq] = uniquetol([real(poles_gi), imag(poles_gi)], 1e-10, 'ByRows', true);
+ poles_gi = poles_gi(i_uniq);
+
+ % Array of distances between all the poles
+ poles_dist = sqrt((poles_est-poles_gi.').*conj(poles_est-poles_gi.'));
+
+ % Get indices corresponding to distances from lowest to highest
+ [~, c] = sort(min(poles_dist));
+
+ as = 1:length(poles_gi);
+
+ % for each column of poles_dist corresponding to the i'th pole
+ % with closest previous poles
+ for p_i = c
+ % Get the indice a_i of the previous pole that is the closest
+ % to pole c(p_i)
+ [~, a_i] = min(poles_dist(:, p_i));
+
+ poles(as(a_i), g_i) = poles_gi(p_i);
+
+ % Remove old poles that are already matched
+ % poles_gi(as(a_i), :) = [];
+ poles_dist(a_i, :) = [];
+ as(a_i) = [];
+ end
+ end
+#+end_src
+
+*** Remove useless poles
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+
+#+begin_src matlab
+ if args.d_max > 0
+ poles = poles(max(abs(poles(:, 2:end) - poles(:, 1:end-1))') > args.d_max, :);
+ end
+
+ if args.p_half
+ poles = poles(1:round(end/2), :);
+ end
+#+end_src
+
+*** Sort poles
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+
+#+begin_src matlab
+ [~, s_p] = sort(imag(poles(:,1)), 'descend');
+ poles = poles(s_p, :);
+#+end_src
+
+*** Transpose for easy plotting
+:PROPERTIES:
+:UNNUMBERED: t
+:END:
+
+#+begin_src matlab
+ poles = poles.';
+#+end_src
diff --git a/paper/paper.org b/paper/paper.org
index e4ed2c9..628646a 100644
--- a/paper/paper.org
+++ b/paper/paper.org
@@ -65,7 +65,7 @@
- $k$: Actuator's Stiffness [N/m]
- $m$: Payload's mass [kg]
- $\omega_0 = \sqrt{\frac{k}{m}}$: Resonance of the (non-rotating) mass-spring system [rad/s]
-- $\omega_r = \dot{\theta}$: rotation speed [rad/s]
+- $\Omega = \dot{\theta}$: rotation speed [rad/s]
#+name: fig:rotating_xy_platform
@@ -141,7 +141,7 @@ Without the coupling terms, each equation is the equation of a one degree of fre
Thus, the term $- m\dot{\theta}^2$ acts like a negative stiffness (due to *centrifugal forces*).
** Constant Rotating Speed
-To simplify, let's consider a constant rotating speed $\dot{\theta} = \omega_r$ and thus $\ddot{\theta} = 0$.
+To simplify, let's consider a constant rotating speed $\dot{\theta} = \Omega$ and thus $\ddot{\theta} = 0$.
#+NAME: eq:coupledplant
\begin{equation}
@@ -154,13 +154,15 @@ To simplify, let's consider a constant rotating speed $\dot{\theta} = \omega_r$
\begin{bmatrix} F_u \\ F_v \end{bmatrix}
\end{equation}
+# Explain each term
+
#+NAME: eq:coupled_plant
\begin{equation}
\begin{bmatrix} d_u \\ d_v \end{bmatrix} =
-\frac{\frac{1}{k}}{\left( \frac{s^2}{{\omega_0}^2} + (1 - \frac{{\omega_r}^2}{{\omega_0}^2}) \right)^2 + \left( 2 \frac{{\omega_r} s}{{\omega_0}^2} \right)^2}
+\frac{\frac{1}{k}}{\left( \frac{s^2}{{\omega_0}^2} + (1 - \frac{{\Omega}^2}{{\omega_0}^2}) \right)^2 + \left( 2 \frac{{\Omega} s}{{\omega_0}^2} \right)^2}
\begin{bmatrix}
- \frac{s^2}{{\omega_0}^2} + 1 - \frac{{\omega_r}^2}{{\omega_0}^2} & 2 \frac{\omega_r s}{{\omega_0}^2} \\
- -2 \frac{\omega_r s}{{\omega_0}^2} & \frac{s^2}{{\omega_0}^2} + 1 - \frac{{\omega_r}^2}{{\omega_0}^2} \\
+ \frac{s^2}{{\omega_0}^2} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} & 2 \frac{\Omega s}{{\omega_0}^2} \\
+ -2 \frac{\Omega s}{{\omega_0}^2} & \frac{s^2}{{\omega_0}^2} + 1 - \frac{{\Omega}^2}{{\omega_0}^2} \\
\end{bmatrix}
\begin{bmatrix} F_u \\ F_v \end{bmatrix}
\end{equation}
6.5 Compliance
%% Name of the Simulink File
@@ -1459,7 +1741,7 @@ Cdvf.OutputName = {'Dx', 'Dy'};
+
-
Figure 21: Comparison of the obtained Compliance
@@ -1468,11 +1750,11 @@ Cdvf.OutputName = {'Dx', 'Dy'};
-
7 Notations
+ + +
+
8 Function
+
+
+
+
+8.1 Sort Poles for the Root Locus
+
+
+
+
+
++This Matlab function is accessible here. +
+
+
+
+Function description
+
+
+
+
+function [poles] = rootLocusPolesSorted(G, K, gains, args) +% rootLocusPolesSorted - +% +% Syntax: [poles] = rootLocusPolesSorted(G, K, gains, args) +% +% Inputs: +% - G, K, gains, args - +% +% Outputs: +% - poles - ++
+
+
+Optional Parameters
+
+
+
+
+arguments
+ G
+ K
+ gains
+ args.minreal double {mustBeNumericOrLogical} = false
+ args.p_half double {mustBeNumericOrLogical} = false
+ args.d_max double {mustBeNumeric} = -1
+end
+
+
+
+
+Function
+
+
+
+
+
+if args.minreal + p1 = pole(minreal(feedback(G, gains(1)*K))); + [~, i_uniq] = uniquetol([real(p1), imag(p1)], 1e-10, 'ByRows', true); + p1 = p1(i_uniq); + + poles = zeros(length(p1), length(gains)); + poles(:, 1) = p1; +else + p1 = pole(feedback(G, gains(1)*K)); + [~, i_uniq] = uniquetol([real(p1), imag(p1)], 1e-10, 'ByRows', true); + p1 = p1(i_uniq); + + poles = zeros(length(p1), length(gains)); + poles(:, 1) = p1; +end ++
+
+
+if args.minreal + p2 = pole(minreal(feedback(G, gains(2)*K))); + [~, i_uniq] = uniquetol([real(p2), imag(p2)], 1e-10, 'ByRows', true); + p2 = p2(i_uniq); + poles(:, 2) = p2; +else + p2 = pole(feedback(G, gains(2)*K)); + [~, i_uniq] = uniquetol([real(p2), imag(p2)], 1e-10, 'ByRows', true); + p2 = p2(i_uniq); + poles(:, 2) = p2; +end ++
+
+for g_i = 3:length(gains) + % Estimated value of the poles + poles_est = poles(:, g_i-1) + (poles(:, g_i-1) - poles(:, g_i-2))*(gains(g_i) - gains(g_i-1))/(gains(g_i-1) - gains(g_i - 2)); + + % New values for the poles + poles_gi = pole(feedback(G, gains(g_i)*K)); + [~, i_uniq] = uniquetol([real(poles_gi), imag(poles_gi)], 1e-10, 'ByRows', true); + poles_gi = poles_gi(i_uniq); + + % Array of distances between all the poles + poles_dist = sqrt((poles_est-poles_gi.').*conj(poles_est-poles_gi.')); + + % Get indices corresponding to distances from lowest to highest + [~, c] = sort(min(poles_dist)); + + as = 1:length(poles_gi); + + % for each column of poles_dist corresponding to the i'th pole + % with closest previous poles + for p_i = c + % Get the indice a_i of the previous pole that is the closest + % to pole c(p_i) + [~, a_i] = min(poles_dist(:, p_i)); + + poles(as(a_i), g_i) = poles_gi(p_i); + + % Remove old poles that are already matched + % poles_gi(as(a_i), :) = []; + poles_dist(a_i, :) = []; + as(a_i) = []; + end +end ++
+
+
+Remove useless poles
+
+
+
+
+if args.d_max > 0 + poles = poles(max(abs(poles(:, 2:end) - poles(:, 1:end-1))') > args.d_max, :); +end + +if args.p_half + poles = poles(1:round(end/2), :); +end ++
+
+Sort poles
+
+
+
+
+[~, s_p] = sort(imag(poles(:,1)), 'descend'); +poles = poles(s_p, :); ++
-
diff --git a/matlab/index.org b/matlab/index.org
index ba137be..d2ad407 100644
--- a/matlab/index.org
+++ b/matlab/index.org
@@ -49,17 +49,14 @@
#+begin_src matlab
addpath('./matlab/');
-#+end_src
-
-#+begin_src matlab
- open('rotating_frame.slx');
+ addpath('./src/');
#+end_src
** System description
The system consists of one 2 degree of freedom translation stage on top of a spindle (figure [[fig:rotating_xy_platform]]).
#+name: fig:rotating_xy_platform
-#+caption: Figure caption
+#+caption: Schematic of the studied system
[[file:figs-tikz/rotating_xy_platform.png]]
The control inputs are the forces applied by the actuators of the translation stage ($F_u$ and $F_v$).
@@ -72,16 +69,29 @@ Based on the Figure [[fig:rotating_xy_platform]], the equations of motions are:
#+begin_important
\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}
+\bm{G}_d
+\begin{bmatrix} F_u \\ F_v \end{bmatrix}
+\end{equation}
+
+\begin{equation}
+\begin{bmatrix} d_u \\ d_v \end{bmatrix} =
+\frac{1}{k} \frac{1}{G_{dp}}
\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} \\
+ G_{dz} & G_{dc} \\
+ -G_{dc} & G_{dz}
\end{bmatrix}
\begin{bmatrix} F_u \\ F_v \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}
#+end_important
Explain Coriolis and Centrifugal Forces (negative Stiffness)
+=> First write the equations in terms of $k$, $m$ and $c$ and explain the terms.
** Numerical Values
Let's define initial values for the model.
@@ -140,7 +150,7 @@ It is shown in Figure [[fig:campbell_diagram]], and one can see that the system
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
- exportFig('figs/campbell_diagram.pdf', 'width', 'full', 'height', 'tall');
+ exportFig('figs/campbell_diagram.pdf', 'width', 'full', 'height', 'normal');
#+end_src
#+name: fig:campbell_diagram
@@ -148,7 +158,7 @@ It is shown in Figure [[fig:campbell_diagram]], and one can see that the system
#+RESULTS:
[[file:figs/campbell_diagram.png]]
-#+begin_src matlab :exports none
+#+begin_src matlab :exports none :tangle no
figure;
ax1 = subplot(1,2,1);
@@ -159,10 +169,12 @@ It is shown in Figure [[fig:campbell_diagram]], and one can see that the system
plot(Ws, zeros(size(Ws)), 'k--')
hold off;
xlabel('Rotation Frequency $\Omega$'); ylabel('Real Part');
- xticks([0, w0, 2*w0])
- xticklabels({'$0$', '$\omega_0$', '$2\omega_0$'})
- yticks([-xi, 0])
- yticklabels({'$-\xi$', '$0$'})
+ xlim([0, 2*w0]);
+ xticks([0,w0/2,w0,3/2*w0,2*w0])
+ xticklabels({'$0$', '', '$\omega_0$', '', '$2 \omega_0$'})
+ ylim([-3*xi, xi]);
+ yticks([-3*xi, -2*xi, -xi, 0, xi])
+ yticklabels({'', '', '$-\xi$', '$0$', ''})
ax2 = subplot(1,2,2);
hold on;
@@ -173,10 +185,16 @@ It is shown in Figure [[fig:campbell_diagram]], and one can see that the system
plot(Ws, zeros(size(Ws)), 'k--')
hold off;
xlabel('Rotation Frequency $\Omega$'); ylabel('Imaginary Part');
- xticks([0, w0, 2*w0])
- xticklabels({'$0$', '$\omega_0$', '$2 \omega_0$'})
- yticks([-w0, 0, w0])
- yticklabels({'$-\omega_0$', '$0$', '$\omega_0$'})
+ xlim([0, 2*w0]);
+ xticks([0,w0/2,w0,3/2*w0,2*w0])
+ xticklabels({'$0$', '', '$\omega_0$', '', '$2 \omega_0$'})
+ ylim([-3*w0, 3*w0]);
+ yticks([-3*w0, -2*w0, -w0, 0, w0, 2*w0, 3*w0])
+ yticklabels({'', '', '$-\omega_0$', '$0$', '$\omega_0$', '', ''})
+#+end_src
+
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/campbell_diagram.pdf', 'width', 'full', 'height', 'normal', 'png', false, 'pdf', false, 'svg', true);
#+end_src
** Simscape Model
@@ -193,6 +211,10 @@ Define the rotating speed for the Simscape Model.
cp = 0; % Parallel Damping [N/(m/s)]
#+end_src
+#+begin_src matlab
+ open('rotating_frame.slx');
+#+end_src
+
The transfer function from $[F_u, F_v]$ to $[d_u, d_v]$ is identified from the Simscape model.
#+begin_src matlab
@@ -213,7 +235,6 @@ The transfer function from $[F_u, F_v]$ to $[d_u, d_v]$ is identified from the S
G.OutputName = {'du', 'dv'};
#+end_src
-** 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.
#+begin_src matlab
Gth = (1/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ...
@@ -285,7 +306,7 @@ Both transfer functions are compared in Figure [[fig:plant_simscape_analytical]]
** 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.
#+begin_src matlab
- Ws = [0, 0.1, 0.5, 0.8, 1.1]*w0; % [rad/s]
+ Ws = [0, 0.2, 0.7, 1.1]*w0; % Rotating Speeds [rad/s]
#+end_src
#+begin_src matlab
@@ -310,7 +331,7 @@ They are compared in Figure [[fig:plant_compare_rotating_speed]].
hold on;
for W_i = 1:length(Ws)
plot(freqs, abs(squeeze(freqresp(Gs{W_i}(1,1), freqs))), ...
- 'DisplayName', sprintf('$\\omega = %.1f \\omega_0 $', Ws(W_i)/w0))
+ 'DisplayName', sprintf('$\\Omega = %.1f \\omega_0 $', Ws(W_i)/w0))
end
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
@@ -332,8 +353,7 @@ They are compared in Figure [[fig:plant_compare_rotating_speed]].
ax2 = subplot(2, 2, 2);
hold on;
for W_i = 1:length(Ws)
- plot(freqs, abs(squeeze(freqresp(Gs{W_i}(2,1), freqs))), ...
- 'DisplayName', sprintf('$\\omega = %.2f \\omega_0 $', Ws(W_i)/w0))
+ plot(freqs, abs(squeeze(freqresp(Gs{W_i}(2,1), freqs))))
end
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
@@ -365,6 +385,10 @@ They are compared in Figure [[fig:plant_compare_rotating_speed]].
#+RESULTS:
[[file:figs/plant_compare_rotating_speed.png]]
+#+begin_src matlab :tangle no :exports none :results none
+ exportFig('figs-inkscape/plant_compare_rotating_speed.pdf', 'width', 'full', 'height', 'full', 'png', false, 'pdf', false, 'svg', true);
+#+end_src
+
* Problem with pure Integral Force Feedback
<Created: 2020-06-12 ven. 19:29
+Created: 2020-06-18 jeu. 11:51