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: Schematic of the studied system
@@ -142,11 +136,11 @@ The measurement is either the \(x-y\) displacement of the object located on top1.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:
1.3 Numerical Values
+1.3 Numerical Values
Let’s define initial values for the model. @@ -204,19 +198,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
@@ -224,8 +218,8 @@ It is shown in Figure 2, and one can see that the syst1.5 Simscape Model
+1.5 Simscape Model
Define the rotating speed for the Simscape Model. @@ -275,11 +269,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
@@ -287,8 +281,8 @@ Both transfer functions are compared in Figure 3 and a1.6 Effect of the rotation speed
+1.6 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. @@ -312,11 +306,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
@@ -325,11 +319,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 @@ -337,8 +331,8 @@ They are compared in Figure 4.
2.1 Schematic
+2.1 Schematic
figs-tikz/system_iff.pdf @@ -346,8 +340,8 @@ They are compared in Figure 4.
2.2 Plant Parameters
+2.2 Plant Parameters
Let’s define initial values for the model. @@ -367,8 +361,8 @@ w0 = sqrt(k/m); % [rad/s]
2.3 Equations
+2.3 Equations
The sensed forces are equal to: @@ -413,8 +407,8 @@ Which then gives:
2.4 Simscape Model
+2.4 Comparison of the Analytical Model and the Simscape Model
The rotation speed is set to \(\Omega = 0.1 \omega_0\). @@ -451,12 +445,7 @@ Giff.InputName = {'Fu', 'Fv'}; Giff.OutputName = {'fu', 'fv'};
2.5 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.
@@ -468,11 +457,10 @@ 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
@@ -480,9 +468,9 @@ The two are compared in Figure 5 and found to perfectl2.6 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.
@@ -505,10 +493,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
@@ -516,9 +504,9 @@ The obtained transfer functions are shown in Figure 6.2.7 Decentralized Integral Force Feedback
-2.6 Decentralized Integral Force Feedback
+The decentralized IFF controller consists of pure integrators:
@@ -530,12 +518,12 @@ The decentralized IFF controller consists of pure integrators: \end{equation}-The Root Locus (evolution of the poles of the closed loop system in the complex plane as a function of \(g\)) is shown in Figure 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.
@@ -544,19 +532,19 @@ It is shown that for non-null rotating speed, one pole is bound to the right-hal3 Integral Force Feedback with an High Pass Filter
+3 Integral Force Feedback with an High Pass Filter
-3.1 Plant Parameters
+3.1 Plant Parameters
Let’s define initial values for the model. @@ -576,8 +564,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: @@ -612,10 +600,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
@@ -623,15 +611,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
@@ -639,11 +627,11 @@ As shown in the Root Locus plot (Figure 9), for some v3.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] @@ -651,7 +639,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)
@@ -674,40 +662,36 @@ The gain at which the system becomes unstable is \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]
+wis = logspace(-2, 1, 100)*w0; % [rad/s]
opt_zeta = zeros(1, length(wis)); % Optimal simultaneous damping
opt_gain = zeros(1, length(wis)); % Corresponding optimal gain
for wi_i = 1:length(wis)
wi = wis(wi_i);
- gains = linspace(0, (w0^2/W^2 - 1)*wi, 100);
+ Kiff = 1/(s + wi)*eye(2);
- for g = gains
- Kiff = (g/(wi+s))*eye(2);
+ fun = @(g)computeSimultaneousDamping(g, Giff, Kiff);
- [w, zeta] = damp(minreal(feedback(Giff, Kiff)));
-
- if min(zeta) > opt_zeta(wi_i) && all(zeta > 0)
- opt_zeta(wi_i) = min(zeta);
- opt_gain(wi_i) = g;
- end
- end
+ [g_opt, xi_opt] = fminsearch(fun, 0.5*wi*((w0/W)^2 - 1));
+ opt_zeta(wi_i) = 1/xi_opt;
+ opt_gain(wi_i) = g_opt;
end
Figure 11: Simultaneous attainable damping of the closed loop poles as a function of \(\omega_i\)
@@ -716,19 +700,19 @@ end4 IFF with a stiffness in parallel with the force sensor
+4 IFF with a stiffness in parallel with the force sensor
-4.1 Schematic
+4.1 Schematic
Figure 12: Figure caption
@@ -736,8 +720,8 @@ end4.2 Equations
+4.2 Equations
4.3 Physical Explanation
+4.3 Physical Explanation
- Negative stiffness induced by gyroscopic effects @@ -816,8 +800,8 @@ And thus, the stiffness in parallel should be such that:
4.4 Plant Parameters
+4.4 Plant Parameters
Let’s define initial values for the model. @@ -837,8 +821,8 @@ w0 = sqrt(k/m); % [rad/s]
4.5 Comparison of the Analytical Model and the Simscape Model
+4.5 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. @@ -886,7 +870,7 @@ Giff_th.OutputName = {'fu', 'fv'};
Figure 13: Comparison of the transfer functions from \([F_u, F_v]\) to \([f_u, f_v]\) between the Simscape model and the analytical one
@@ -894,8 +878,8 @@ Giff_th.OutputName = {'fu', 'fv'};4.6 Effect of the parallel stiffness on the IFF plant
+4.6 Effect of the parallel stiffness on the IFF plant
The rotation speed is set to \(\Omega = 0.1 \omega_0\). @@ -915,7 +899,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 14. +The results are shown in Figure 14.
@@ -924,7 +908,6 @@ One can see that for \(k_p > m \Omega^2\), the systems shows alternating complex
kp = 0; -cp = 0; w0p = sqrt((k + kp)/m); xip = c/(2*sqrt((k+kp)*m)); @@ -937,7 +920,7 @@ Giff = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 )kp = 0.5*m*W^2; -cp = 0; +k = 1 - kp; w0p = sqrt((k + kp)/m); xip = c/(2*sqrt((k+kp)*m)); @@ -950,7 +933,7 @@ Giff_s = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2-kp = 1.5*m*W^2; -cp = 0; +k = 1 - kp; w0p = sqrt((k + kp)/m); xip = c/(2*sqrt((k+kp)*m)); @@ -962,7 +945,7 @@ Giff_l = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2+-
Figure 14: 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\)
@@ -970,11 +953,11 @@ Giff_l = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2-4.7 IFF when adding a spring in parallel
++4.7 IFF when adding a spring in parallel
-In Figure 15 is displayed the Root Locus in the three considered cases with +In Figure 15 is displayed the Root Locus in the three considered cases with
\begin{equation} K_{\text{IFF}} = \frac{g}{s} \begin{bmatrix} @@ -995,7 +978,7 @@ Thus, decentralized IFF controller with pure integrators can be used if: \end{equation} -+-
Figure 15: Root Locus
@@ -1003,8 +986,8 @@ Thus, decentralized IFF controller with pure integrators can be used if:-4.8 Effect of \(k_p\) on the attainable damping
++4.8 Effect of \(k_p\) on the attainable damping
However, having large values of \(k_p\) may: @@ -1015,7 +998,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 16. +To study the second point, Root Locus plots for the following values of \(k_p\) are shown in Figure 16.
-kps = [2, 20, 40]*m*W^2; @@ -1027,16 +1010,50 @@ It is shown that large values of \(k_p\) decreases the attainable damping. -++ +
Figure 16: Root Locus plot
++ + +alphas = logspace(-2, 0, 100); + +opt_zeta = zeros(1, length(alphas)); % Optimal simultaneous damping +opt_gain = zeros(1, length(alphas)); % Corresponding optimal gain + +Kiff = 1/s*eye(2); + +for alpha_i = 1:length(alphas) + kp = alphas(alpha_i); + k = 1 - alphas(alpha_i); + + w0p = sqrt((k + kp)/m); + xip = c/(2*sqrt((k+kp)*m)); + + Giff = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ... + (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)); + (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2]; + + fun = @(g)computeSimultaneousDamping(g, Giff, Kiff); + + [g_opt, xi_opt] = fminsearch(fun, 2); + opt_zeta(alpha_i) = 1/xi_opt; + opt_gain(alpha_i) = g_opt; +end ++++
+
-4.9 Optimal Gain
++-4.9 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. @@ -1044,7 +1061,7 @@ Let’s take \(k_p = 5 m \Omega^2\) and find the optimal IFF control gain \(
-kp = 5*m*W^2; -cp = 0.01; +k = 1 - kp; w0p = sqrt((k + kp)/m); xip = c/(2*sqrt((k+kp)*m)); @@ -1075,231 +1092,56 @@ end+-
Figure 17: Root Locus for \(k_p = 5 m \Omega^2\) and the poles corresponding to the identified optimal gain
+Figure 18: Root Locus for \(k_p = 5 m \Omega^2\) and the poles corresponding to the identified optimal gain
-5 Direct Velocity Feedback
++5 Comparison
--5.1 Schematic
++5.1 Plant Parameters
- ---
+
+Let’s define initial values for the model.
-Figure 18: Figure caption
+++ +k = 1; % Actuator Stiffness [N/m] +c = 0.05; % Actuator Damping [N/(m/s)] +m = 1; % Payload mass [kg] ++++ +xi = c/(2*sqrt(k*m)); +w0 = sqrt(k/m); % [rad/s] +++The rotating speed is set to \(\Omega = 0.1 \omega_0\). +
++W = 0.1*w0; +-- -5.2 Equations
++- -5.2 Root Locus
--The sensed relative velocity are equal to: -
--\begin{equation} -\begin{bmatrix} \dot{d}_u \\ \dot{d}_v \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.3 Plant Parameters
----Let’s define initial values for the model. -
--- -k = 1; % Actuator Stiffness [N/m] -c = 0.05; % Actuator Damping [N/(m/s)] -m = 1; % Payload mass [kg] ----xi = c/(2*sqrt(k*m)); -w0 = sqrt(k/m); % [rad/s] ---- -5.4 Comparison of the Analytical Model and the Simscape Model
----The rotating speed is set to \(\Omega = 0.1 \omega_0\). -
--- -W = 0.1*w0; ---- -open('rotating_frame.slx'); ---And the transfer function from \([F_u, F_v]\) to \([v_u, v_v]\) is identified using the Simscape model. -
--- -%% Name of the Simulink File -mdl = 'rotating_frame'; - -%% Input/Output definition -clear io; io_i = 1; -io(io_i) = linio([mdl, '/K'], 1, 'openinput'); io_i = io_i + 1; -io(io_i) = linio([mdl, '/G'], 1, 'openoutput'); io_i = io_i + 1; ---- -Gdvf = linearize(mdl, io, 0); - -%% Input/Output definition -Gdvf.InputName = {'Fu', 'Fv'}; -Gdvf.OutputName = {'Vu', 'Vv'}; ---The same transfer function from \([F_u, F_v]\) to \([v_u, v_v]\) is written down from the analytical model. -
--- -Gdvf_th = (s/k)/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ... - [(s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2), 2*W*s/(w0^2) ; ... - -2*W*s/(w0^2), (s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)]; - -Gdvf_th.InputName = {'Fu', 'Fv'}; -Gdvf_th.OutputName = {'vu', 'vv'}; ---The two are compared in Figure 5 and found to perfectly match. -
- - ----
-
Figure 19: Comparison of the transfer functions from \([F_u, F_v]\) to \([v_u, v_v]\) between the Simscape model and the analytical one
---5.5 Root Locus
----The Decentralized Direct Velocity Feedback controller consist of a pure gain on the diagonal: -
-\begin{equation} - K_{\text{DVF}}(s) = g \begin{bmatrix} - 1 & 0 \\ - 0 & 1 -\end{bmatrix} -\end{equation} - --The corresponding Root Locus plots for the following rotating speeds are shown in Figure 20. -
--- -Ws = [0, 0.2, 0.7, 1.1]*w0; % Rotating Speeds [rad/s] ---It is shown that for rotating speed \(\Omega < \omega_0\), the closed loop system is unconditionally stable and arbitrary damping can be added to the poles. -
- ----
-
Figure 20: Root Locus for the Decentralized Direct Velocity Feedback controller. Several rotating speed are shown.
--6 Comparison
- - --- -6.1 Plant Parameters
----Let’s define initial values for the model. -
--- -k = 1; % Actuator Stiffness [N/m] -c = 0.05; % Actuator Damping [N/(m/s)] -m = 1; % Payload mass [kg] ---- -xi = c/(2*sqrt(k*m)); -w0 = sqrt(k/m); % [rad/s] ---The rotating speed is set to \(\Omega = 0.1 \omega_0\). -
---W = 0.1*w0; ---6.2 Root Locus
---IFF with High Pass Filter
@@ -1329,28 +1171,18 @@ k = k + kp;--DVF -
--- -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)]; --+-
Figure 21: Root Locus plot - Comparison of IFF with additional high pass filter, IFF with additional parallel stiffness and DVF
+Figure 19: Root Locus plot - Comparison of IFF with additional high pass filter, IFF with additional parallel stiffness
-6.3 Controllers - Optimal Gains
-++-5.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.
@@ -1388,51 +1220,43 @@ The obtained damping ratio and control are shown below.0.83 2.01 - -- DVF -0.85 -1.67 --6.4 Transmissibility
-+++ +5.4 Passive Damping - Critical Damping
+ ++- -5.5 Transmissibility And Compliance
++--open('rotating_frame.slx');-- -%% Name of the Simulink File -mdl = 'rotating_frame'; - -%% Input/Output definition -clear io; io_i = 1; -io(io_i) = linio([mdl, '/dw'], 1, 'input'); io_i = io_i + 1; -io(io_i) = linio([mdl, '/Meas'], 1, 'output'); io_i = io_i + 1; ----Tol = linearize(mdl, io, 0); - -%% Input/Output definition -Tol.InputName = {'Dwx', 'Dwy'}; -Tol.OutputName = {'Dx', 'Dy'}; ---Kiff = opt_gain_iff/(wi + s)*tf(eye(2)); --- -%% Name of the Simulink File @@ -1441,119 +1265,68 @@ mdl = 'rotating_frame'; %% Input/Output definition clear io; io_i = 1; io(io_i) = linio([mdl, '/dw'], 1, 'input'); io_i = io_i + 1; -io(io_i) = linio([mdl, '/Meas'], 1, 'output'); io_i = io_i + 1; ----Tiff = linearize(mdl, io, 0); - -%% Input/Output definition -Tiff.InputName = {'Dwx', 'Dwy'}; -Tiff.OutputName = {'Dx', 'Dy'}; ---- -kp = 5*m*W^2; -cp = 0.01; ---- -Kiff = opt_gain_kp/s*tf(eye(2)); ---- -%% Name of the Simulink File -mdl = 'rotating_frame'; - -%% Input/Output definition -clear io; io_i = 1; -io(io_i) = linio([mdl, '/dw'], 1, 'input'); io_i = io_i + 1; -io(io_i) = linio([mdl, '/Meas'], 1, 'output'); io_i = io_i + 1; ----Tiff_kp = linearize(mdl, io, 0); - -%% Input/Output definition -Tiff_kp.InputName = {'Dwx', 'Dwy'}; -Tiff_kp.OutputName = {'Dx', 'Dy'}; ---- -Kdvf = opt_gain_kp*tf(eye(2)); ---- -%% Name of the Simulink File -mdl = 'rotating_frame'; - -%% Input/Output definition -clear io; io_i = 1; -io(io_i) = linio([mdl, '/dw'], 1, 'input'); io_i = io_i + 1; -io(io_i) = linio([mdl, '/Meas'], 1, 'output'); io_i = io_i + 1; ---- -Tdvf = linearize(mdl, io, 0); - -%% Input/Output definition -Tdvf.InputName = {'Dwx', 'Dwy'}; -Tdvf.OutputName = {'Dx', 'Dy'}; -----
-
Figure 22: Comparison of the transmissibility
--6.5 Compliance
-- -+ +--%% Name of the Simulink File -mdl = 'rotating_frame'; - -%% Input/Output definition -clear io; io_i = 1; io(io_i) = linio([mdl, '/fd'], 1, 'input'); io_i = io_i + 1; io(io_i) = linio([mdl, '/Meas'], 1, 'output'); io_i = io_i + 1;-+Col = linearize(mdl, io, 0); +G_ol = linearize(mdl, io, 0); %% Input/Output definition -Col.InputName = {'Fdx', 'Fdy'}; -Col.OutputName = {'Dx', 'Dy'}; +G_ol.InputName = {'Dwx', 'Dwy', 'Fdx', 'Fdy'}; +G_ol.OutputName = {'Dx', 'Dy'};++ +5.5.1 Passive Damping
+++++ +kp = 0; +cp = 0; ++++ +c_old = c; +c = c_opt; ++++ +G_pas = linearize(mdl, io, 0); + +%% Input/Output definition +G_pas.InputName = {'Dwx', 'Dwy', 'Fdx', 'Fdy'}; +G_pas.OutputName = {'Dx', 'Dy'}; ++++c = c_old; +++ +Kiff = opt_gain_iff/(wi + s)*tf(eye(2));-+Ciff = linearize(mdl, io, 0); +G_iff = linearize(mdl, io, 0); %% Input/Output definition -Ciff.InputName = {'Fdx', 'Fdy'}; -Ciff.OutputName = {'Dx', 'Dy'}; +G_iff.InputName = {'Dwx', 'Dwy', 'Fdx', 'Fdy'}; +G_iff.OutputName = {'Dx', 'Dy'};kp = 5*m*W^2; cp = 0.01; @@ -1566,41 +1339,142 @@ cp = 0.01;--Ciff_kp = linearize(mdl, io, 0); +G_kp = linearize(mdl, io, 0); %% Input/Output definition -Ciff_kp.InputName = {'Fdx', 'Fdy'}; -Ciff_kp.OutputName = {'Dx', 'Dy'}; +G_kp.InputName = {'Dwx', 'Dwy', 'Fdx', 'Fdy'}; +G_kp.OutputName = {'Dx', 'Dy'};--Kdvf = opt_gain_kp*tf(eye(2)); ---Cdvf = linearize(mdl, io, 0); +-%% Input/Output definition -Cdvf.InputName = {'Fdx', 'Fdy'}; -Cdvf.OutputName = {'Dx', 'Dy'}; - +-++
Figure 20: Comparison of the transmissibility
+++ ++ +++-
Figure 23: Comparison of the obtained Compliance
+Figure 21: Comparison of the obtained Compliance
++5.6 DC Compliance
++++wis = logspace(-2, 1, 100)*w0; % [rad/s] + +opt_zeta_wi = zeros(1, length(wis)); % Optimal simultaneous damping +opt_gain_wi = zeros(1, length(wis)); % Corresponding optimal gain +C_dc_gain_wi = zeros(1, length(wis)); % Compliance DC value + +for wi_i = 1:length(wis) + wi = wis(wi_i); + Kiff = 1/(s + wi)*eye(2); + + fun = @(g)computeSimultaneousDamping(g, Gwi({'fu', 'fv'}, {'Fu', 'Fv'}), Kiff); + + [g_opt, xi_opt] = fminsearch(fun, 0.5*wi*((w0/W)^2 - 1)); + opt_zeta_wi(wi_i) = 1/xi_opt; + opt_gain_wi(wi_i) = g_opt; + + G_dc_gain = dcgain(lft(Gwi, -g_opt/(s + wi)*eye(2), 2, 2)); + C_dc_gain_wi(wi_i) = G_dc_gain(1,1); +end ++++ +alphas = logspace(-2, 0, 100); + +opt_zeta_kp = zeros(1, length(alphas)); % Optimal simultaneous damping +opt_gain_kp = zeros(1, length(alphas)); % Corresponding optimal gain +C_dc_gain_kp = zeros(1, length(alphas)); % DC gain of the compliance + +Kiff = 1/s*eye(2); + +for alpha_i = 1:length(alphas) + kp = alphas(alpha_i); + k = 1 - alphas(alpha_i); + + w0p = sqrt((k + kp)/m); + xip = c/(2*sqrt((k+kp)*m)); + + Gkp = tf(zeros(4,4)); + Gkp.InputName = {'Fx', 'Fy', 'Fu', 'Fv'}; + Gkp.OutputName = {'dx', 'dy', 'fu', 'fv'}; + + Gp = ((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))^2 + (2*W*s/(w0p^2))^2; + + Gkp('dx', 'Fu') = (1/(k + kp))*((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))/Gp; + Gkp('dy', 'Fv') = (1/(k + kp))*((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))/Gp; + + Gkp('dx', 'Fx') = (1/(k + kp))*((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))/Gp; + Gkp('dy', 'Fy') = (1/(k + kp))*((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))/Gp; + + Gkp('dx', 'Fv') = (1/(k + kp))*(2*W*s/(w0p^2))/Gp; + Gkp('dy', 'Fu') = -(1/(k + kp))*(2*W*s/(w0p^2))/Gp; + + Gkp('dx', 'Fy') = (1/(k + kp))*(2*W*s/(w0p^2))/Gp; + Gkp('dy', 'Fx') = -(1/(k + kp))*(2*W*s/(w0p^2))/Gp; + + Gkp('fu', 'Fu') = ((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)/Gp; + Gkp('fv', 'Fv') = ((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)/Gp; + + Gkp('fu', 'Fv') = -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p))/Gp; + Gkp('fv', 'Fu') = (2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p))/Gp; + + Gkp('fu', 'Fx') = -(c*s + k)*(1/(k + kp))*((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))/Gp; + Gkp('fv', 'Fy') = -(c*s + k)*(1/(k + kp))*((s^2)/(w0p^2) + 2*xip*s/w0p + 1 - (W^2)/(w0p^2))/Gp; + + Gkp('fu', 'Fy') = -(c*s + k)*(1/(k + kp))*(2*W*s/(w0p^2))/Gp; + Gkp('fv', 'Fx') = (c*s + k)*(1/(k + kp))*(2*W*s/(w0p^2))/Gp; + + fun = @(g)computeSimultaneousDamping(g, Gkp({'fu', 'fv'}, {'Fu', 'Fv'}), Kiff); + + [g_opt, xi_opt] = fminsearch(fun, 2); + opt_zeta_kp(alpha_i) = 1/xi_opt; + opt_gain_kp(alpha_i) = g_opt; + + G_dc_gain = dcgain(lft(Gkp, -g_opt/s*eye(2), 2, 2)); + C_dc_gain_kp(alpha_i) = G_dc_gain(1,1); +end ++++ + ++
+
++ + ++
+
+
+
-7 Notations
-++6 Notations
+@@ -1802,7 +1676,7 @@ Cdvf.OutputName = {'Dx', 'Dy'};
-diff --git a/matlab/src/computeSimultaneousDamping.m b/matlab/src/computeSimultaneousDamping.m new file mode 100644 index 0000000..2cfb68d --- /dev/null +++ b/matlab/src/computeSimultaneousDamping.m @@ -0,0 +1,10 @@ +% Compute Damping + +function [xi_min] = computeSimultaneousDamping(g, G, K) + [w, xi] = damp(minreal(feedback(G, g*K))); + xi_min = 1/min(xi); + + if xi_min < 0 + xi_min = 1e8; + end +endCreated: 2020-06-24 mer. 08:44
+Created: 2020-07-06 lun. 23:23