Update Matlab code
This commit is contained in:
File diff suppressed because it is too large
Load Diff
463
matlab/index.org
463
matlab/index.org
@@ -1,4 +1,4 @@
|
||||
#+TITLE: Matlab Computation
|
||||
#+TITLE: Active Damping of Rotating Platforms using Integral Force Feedback - Matlab Computation
|
||||
:DRAWER:
|
||||
#+HTML_LINK_HOME: ../index.html
|
||||
#+HTML_LINK_UP: ../index.html
|
||||
@@ -24,13 +24,15 @@
|
||||
:END:
|
||||
|
||||
* Introduction :ignore:
|
||||
This document gathers the Matlab code used to for the paper cite:dehaeze20_activ_dampin_rotat_platf_integ_force_feedb.
|
||||
|
||||
- Section [[sec:system_description]]
|
||||
- Section [[sec:iff_pure_int]]
|
||||
- Section [[sec:iff_pseudo_int]]
|
||||
- Section [[sec:iff_parallel_stiffness]]
|
||||
- Section [[sec:comparison]]
|
||||
- Section [[sec:notations]]
|
||||
It is structured in several sections:
|
||||
- Section [[sec:system_description]]: presents a simple model of a rotating suspended platform that will be used throughout this study.
|
||||
- Section [[sec:iff_pure_int]]: explains how the unconditional stability of IFF is lost due to Gyroscopic effects induced by the rotation.
|
||||
- Section [[sec:iff_pseudo_int]]: suggests a simple modification of the control law such that damping can be added to the suspension modes in a robust way.
|
||||
- Section [[sec:iff_parallel_stiffness]]: proposes to add springs in parallel with the force sensors to regain the unconditional stability of IFF.
|
||||
- Section [[sec:comparison]]: compares both proposed modifications to the classical IFF in terms of damping authority and closed-loop system behavior.
|
||||
- Section [[sec:notations]]: contains the notations used for both the Matlab code and the paper
|
||||
|
||||
* System Description and Analysis
|
||||
:PROPERTIES:
|
||||
@@ -38,7 +40,6 @@
|
||||
:END:
|
||||
<<sec:system_description>>
|
||||
|
||||
** Introduction :ignore:
|
||||
** Matlab Init :noexport:ignore:
|
||||
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
|
||||
<<matlab-dir>>
|
||||
@@ -63,8 +64,6 @@ The system consists of one 2 degree of freedom translation stage on top of a spi
|
||||
The control inputs are the forces applied by the actuators of the translation stage ($F_u$ and $F_v$).
|
||||
As the translation stage is rotating around the Z axis due to the spindle, the forces are applied along $\vec{i}_u$ and $\vec{i}_v$.
|
||||
|
||||
The measurement is either the $x-y$ displacement of the object located on top of the translation stage or the $u-v$ displacement of the sample with respect to a fixed reference frame.
|
||||
|
||||
** Equations
|
||||
Based on the Figure [[fig:system]], the equations of motions are:
|
||||
#+begin_important
|
||||
@@ -90,9 +89,6 @@ With:
|
||||
\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.
|
||||
#+begin_src matlab
|
||||
@@ -207,7 +203,10 @@ It is shown in Figure [[fig:campbell_diagram]], and one can see that the system
|
||||
#+end_src
|
||||
|
||||
** Simscape Model
|
||||
Define the rotating speed for the Simscape Model.
|
||||
In order to validate all the equations of motion, a Simscape model of the same system has been developed.
|
||||
The dynamics of the system can be identified from the Simscape model and compare with the analytical model.
|
||||
|
||||
The rotating speed for the Simscape Model is defined.
|
||||
#+begin_src matlab
|
||||
W = 0.1; % Rotation Speed [rad/s]
|
||||
#+end_src
|
||||
@@ -404,8 +403,13 @@ They are compared in Figure [[fig:plant_compare_rotating_speed]].
|
||||
<<sec:iff_pure_int>>
|
||||
|
||||
** Introduction :ignore:
|
||||
- Diagram with the controller
|
||||
- Basic idea of IFF
|
||||
Force sensors are added in series with the two actuators (Figure [[fig:system_iff]]).
|
||||
|
||||
Two identical controllers $K_F$ are used to feedback each of the sensed force to its associated actuator.
|
||||
|
||||
#+name: fig:system_iff
|
||||
#+caption: System with added Force Sensor in series with the actuators
|
||||
[[file:figs-paper/system_iff.png]]
|
||||
|
||||
** Matlab Init :noexport:ignore:
|
||||
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
|
||||
@@ -421,12 +425,6 @@ They are compared in Figure [[fig:plant_compare_rotating_speed]].
|
||||
addpath('./src/');
|
||||
#+end_src
|
||||
|
||||
** Schematic
|
||||
|
||||
#+name: fig:system_iff
|
||||
#+caption: Figure caption
|
||||
[[file:figs-paper/system_iff.pdf]]
|
||||
|
||||
** Plant Parameters
|
||||
Let's define initial values for the model.
|
||||
#+begin_src matlab
|
||||
@@ -732,9 +730,6 @@ It is shown that for non-null rotating speed, one pole is bound to the right-hal
|
||||
:END:
|
||||
<<sec:iff_pseudo_int>>
|
||||
|
||||
** Introduction :ignore:
|
||||
- Classical modification of the IFF
|
||||
|
||||
** Matlab Init :noexport:ignore:
|
||||
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
|
||||
<<matlab-dir>>
|
||||
@@ -1118,7 +1113,6 @@ The gain at which the system becomes unstable is
|
||||
|
||||
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.
|
||||
|
||||
|
||||
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 [[fig:mod_iff_damping_wi]]).
|
||||
|
||||
@@ -1180,7 +1174,6 @@ To find the optimum, the gain that maximize the simultaneous damping of the mode
|
||||
:END:
|
||||
<<sec:iff_parallel_stiffness>>
|
||||
|
||||
** Introduction :ignore:
|
||||
** Matlab Init :noexport:ignore:
|
||||
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
|
||||
<<matlab-dir>>
|
||||
@@ -1196,11 +1189,17 @@ To find the optimum, the gain that maximize the simultaneous damping of the mode
|
||||
#+end_src
|
||||
|
||||
** Schematic
|
||||
In this section additional springs in parallel with the force sensors are added to counteract the negative stiffness induced by the rotation.
|
||||
|
||||
#+name: fig:system_parallel_springs
|
||||
#+caption: Figure caption
|
||||
#+caption: Studied system with additional springs in parallel with the actuators and force sensors
|
||||
[[file:figs-paper/system_parallel_springs.png]]
|
||||
|
||||
In order to keep the overall stiffness $k = k_a + k_p$ constant, a scalar parameter $\alpha$ ($0 \le \alpha < 1$) is defined to describe the fraction of the total stiffness in parallel with the actuator and force sensor
|
||||
\begin{equation}
|
||||
k_p = \alpha k, \quad k_a = (1 - \alpha) k
|
||||
\end{equation}
|
||||
|
||||
** Equations
|
||||
#+begin_important
|
||||
\begin{equation}
|
||||
@@ -1220,40 +1219,16 @@ To find the optimum, the gain that maximize the simultaneous damping of the mode
|
||||
\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)
|
||||
G_{kp} &= \left( \frac{s^2}{{\omega_0}^2} + 2\xi \frac{s}{{\omega_0}^2} + 1 - \frac{\Omega^2}{{\omega_0}^2} \right)^2 + \left( 2 \frac{\Omega}{\omega_0}\frac{s}{\omega_0} \right)^2 \\
|
||||
G_{kz} &= \left( \frac{s^2}{{\omega_0}^2} - \frac{\Omega^2}{{\omega_0}^2} + \alpha \right) \left( \frac{s^2}{{\omega_0}^2} + 2\xi \frac{s}{{\omega_0}^2} + 1 - \frac{\Omega^2}{{\omega_0}^2} \right) + \left( 2 \frac{\Omega}{\omega_0}\frac{s}{\omega_0} \right)^2 \\
|
||||
G_{kc} &= \left( 2 \xi \frac{s}{\omega_0} + 1 - \alpha \right) \left( 2 \frac{\Omega}{\omega_0}\frac{s}{\omega_0} \right)
|
||||
\end{align}
|
||||
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}$.
|
||||
If we compare $G_{kz}$ and $G_{fz}$, we see that the spring in parallel adds a term $\alpha$.
|
||||
In order to have two complex conjugate zeros (instead of real zeros):
|
||||
\begin{equation}
|
||||
\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
|
||||
\alpha > \frac{\Omega^2}{{\omega_0}^2} \quad \Leftrightarrow \quad k_p > m \Omega^2
|
||||
\end{equation}
|
||||
|
||||
** Plant Parameters
|
||||
@@ -1684,9 +1659,7 @@ Thus, decentralized IFF controller with pure integrators can be used if:
|
||||
#+end_src
|
||||
|
||||
** Effect of $k_p$ on the attainable damping
|
||||
However, having large values of $k_p$ may:
|
||||
- decrease the actuator force authority
|
||||
- decrease the attainable damping
|
||||
However, having large values of $k_p$ may decrease the attainable damping.
|
||||
|
||||
To study the second point, Root Locus plots for the following values of $k_p$ are shown in Figure [[fig:root_locus_iff_kps]].
|
||||
#+begin_src matlab
|
||||
@@ -1694,7 +1667,6 @@ To study the second point, Root Locus plots for the following values of $k_p$ ar
|
||||
#+end_src
|
||||
|
||||
It is shown that large values of $k_p$ decreases the attainable damping.
|
||||
|
||||
#+begin_src matlab :exports none
|
||||
figure;
|
||||
|
||||
@@ -1842,114 +1814,15 @@ It is shown that large values of $k_p$ decreases the attainable damping.
|
||||
#+RESULTS:
|
||||
[[file:figs/opt_damp_alpha.png]]
|
||||
|
||||
** TODO 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;
|
||||
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 ];
|
||||
#+end_src
|
||||
|
||||
#+begin_src matlab
|
||||
opt_xi = 0;
|
||||
opt_gain = 0;
|
||||
|
||||
gains = logspace(-2, 4, 1000);
|
||||
|
||||
for g = gains
|
||||
Kiff = (g/s)*eye(2);
|
||||
|
||||
[w, xi] = damp(minreal(feedback(Giff, Kiff)));
|
||||
|
||||
if min(xi) > opt_xi && all(xi > 0)
|
||||
opt_xi = min(xi);
|
||||
opt_gain = min(g);
|
||||
end
|
||||
end
|
||||
#+end_src
|
||||
|
||||
#+begin_src matlab :exports none
|
||||
figure;
|
||||
|
||||
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(minreal(feedback(Giff, (g/s)*eye(2))));
|
||||
plot(real(clpoles), imag(clpoles), 'k.');
|
||||
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 results :results file replace
|
||||
exportFig('figs/root_locus_opt_gain_iff_kp.pdf', 'width', 'wide', 'height', 'normal');
|
||||
#+end_src
|
||||
|
||||
#+name: fig:root_locus_opt_gain_iff_kp
|
||||
#+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-paper/root_locus_opt_gain_iff_kp.pdf', 'width', 'wide', 'height', 'tall', 'png', false, 'pdf', false, 'svg', true);
|
||||
#+end_src
|
||||
|
||||
* Comparison
|
||||
:PROPERTIES:
|
||||
:header-args:matlab+: :tangle matlab/s6_act_damp_comparison.m
|
||||
:header-args:matlab+: :tangle matlab/s5_act_damp_comparison.m
|
||||
:END:
|
||||
<<sec:comparison>>
|
||||
|
||||
** Introduction :ignore:
|
||||
Two modifications to adapt the IFF control strategy to rotating platforms have been proposed.
|
||||
These two methods are now compared in terms of added damping, closed-loop compliance and transmissibility.
|
||||
|
||||
** Matlab Init :noexport:ignore:
|
||||
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
|
||||
@@ -2315,262 +2188,6 @@ Critical Damping corresponds to to $\xi = 1$, and thus:
|
||||
exportFig('figs-paper/comp_compliance.pdf', 'width', 'half', 'height', 'normal', 'png', false, 'pdf', false, 'svg', true);
|
||||
#+end_src
|
||||
|
||||
** DC Compliance
|
||||
*** Pseudo Integrator IFF :ignore:
|
||||
#+begin_src matlab :exports none
|
||||
k = 1;
|
||||
m = 1;
|
||||
w0 = sqrt(k/m);
|
||||
#+end_src
|
||||
|
||||
#+begin_src matlab :exports none
|
||||
Gwi = tf(zeros(4,4));
|
||||
Gwi.InputName = {'Fx', 'Fy', 'Fu', 'Fv'};
|
||||
Gwi.OutputName = {'dx', 'dy', 'fu', 'fv'};
|
||||
|
||||
Gp = ((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2;
|
||||
|
||||
Gwi('dx', 'Fu') = (1/k)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))/Gp;
|
||||
Gwi('dy', 'Fv') = (1/k)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))/Gp;
|
||||
|
||||
Gwi('dx', 'Fx') = (1/k)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))/Gp;
|
||||
Gwi('dy', 'Fy') = (1/k)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))/Gp;
|
||||
|
||||
Gwi('dx', 'Fv') = (1/k)*(2*W*s/(w0^2))/Gp;
|
||||
Gwi('dy', 'Fu') = -(1/k)*(2*W*s/(w0^2))/Gp;
|
||||
|
||||
Gwi('dx', 'Fy') = (1/k)*(2*W*s/(w0^2))/Gp;
|
||||
Gwi('dy', 'Fx') = -(1/k)*(2*W*s/(w0^2))/Gp;
|
||||
|
||||
Gwi('fu', 'Fu') = ((s^2/w0^2 - W^2/w0^2)*(s^2/w0^2 + 2*xi*s/w0 + 1 - W^2/w0^2) + (2*(s/w0)*(W/w0))^2)/Gp;
|
||||
Gwi('fv', 'Fv') = ((s^2/w0^2 - W^2/w0^2)*(s^2/w0^2 + 2*xi*s/w0 + 1 - W^2/w0^2) + (2*(s/w0)*(W/w0))^2)/Gp;
|
||||
|
||||
Gwi('fu', 'Fv') = -(2*xi*s/w0 + 1)*(2*(s/w0)*(W/w0))/Gp;
|
||||
Gwi('fv', 'Fu') = (2*xi*s/w0 + 1)*(2*(s/w0)*(W/w0))/Gp;
|
||||
|
||||
Gwi('fu', 'Fx') = -(c*s + k)*(1/k)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))/Gp;
|
||||
Gwi('fv', 'Fy') = -(c*s + k)*(1/k)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))/Gp;
|
||||
|
||||
Gwi('fu', 'Fy') = -(c*s + k)*(1/k)*(2*W*s/(w0^2))/Gp;
|
||||
Gwi('fv', 'Fx') = (c*s + k)*(1/k)*(2*W*s/(w0^2))/Gp;
|
||||
#+end_src
|
||||
|
||||
#+begin_src matlab
|
||||
wis = logspace(-2, 1, 100)*w0; % [rad/s]
|
||||
|
||||
opt_xi_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_xi_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
|
||||
#+end_src
|
||||
|
||||
*** IFF With parallel Stiffness :ignore:
|
||||
#+begin_src matlab
|
||||
alphas = logspace(-2, 0, 100);
|
||||
|
||||
opt_xi_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_xi_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
|
||||
#+end_src
|
||||
|
||||
*** Comparison :ignore:
|
||||
#+begin_src matlab :exports none
|
||||
figure;
|
||||
yyaxis left
|
||||
plot(alphas, C_dc_gain_kp, '-', 'DisplayName', '$|T(0)|$');
|
||||
set(gca, 'YScale', 'log');
|
||||
ylim([0, 1e3]);
|
||||
ylabel('DC value of the Compliance');
|
||||
|
||||
yyaxis right
|
||||
hold on;
|
||||
plot(alphas, opt_xi_kp, '-', 'DisplayName', '$\xi_{opt}$');
|
||||
set(gca, 'YScale', 'lin');
|
||||
ylim([0,1]);
|
||||
ylabel('Optimal Damping Ratio');
|
||||
|
||||
xlabel('$\alpha$');
|
||||
set(gca, 'XScale', 'log');
|
||||
legend('location', 'northeast');
|
||||
#+end_src
|
||||
|
||||
#+begin_src matlab :tangle no :exports results :results file replace
|
||||
exportFig('figs/compliance_dc_gain_wi.pdf', 'width', 'wide', 'height', 'normal');
|
||||
#+end_src
|
||||
|
||||
#+name: fig:compliance_dc_gain_wi
|
||||
#+caption:
|
||||
#+RESULTS:
|
||||
[[file:figs/compliance_dc_gain_wi.png]]
|
||||
|
||||
#+begin_src matlab :exports none
|
||||
figure;
|
||||
yyaxis left
|
||||
plot(wis, abs(C_dc_gain_wi), '-', 'DisplayName', '$|T(0)|$');
|
||||
set(gca, 'YScale', 'log');
|
||||
ylim([0, 1e3]);
|
||||
ylabel('DC value of the Compliance');
|
||||
|
||||
yyaxis right
|
||||
hold on;
|
||||
plot(wis, opt_xi_wi, '-', 'DisplayName', '$\xi_{opt}$');
|
||||
set(gca, 'YScale', 'lin');
|
||||
ylim([0,1]);
|
||||
ylabel('Optimal Damping Ratio');
|
||||
|
||||
xlabel('$\omega_i$');
|
||||
set(gca, 'XScale', 'log');
|
||||
legend('location', 'northeast');
|
||||
#+end_src
|
||||
|
||||
#+begin_src matlab :tangle no :exports results :results file replace
|
||||
exportFig('figs/compliance_dc_gain_kp.pdf', 'width', 'wide', 'height', 'normal');
|
||||
#+end_src
|
||||
|
||||
#+name: fig:compliance_dc_gain_kp
|
||||
#+caption:
|
||||
#+RESULTS:
|
||||
[[file:figs/compliance_dc_gain_kp.png]]
|
||||
|
||||
#+begin_src matlab :exports none
|
||||
figure;
|
||||
hold on;
|
||||
plot(opt_xi_wi, C_dc_gain_wi, '-', 'DisplayName', '$\omega_i$');
|
||||
plot(opt_xi_kp, C_dc_gain_kp, '-', 'DisplayName', '$k_p$');
|
||||
hold off
|
||||
set(gca, 'YScale', 'log');
|
||||
ylim([0, 1e3]);
|
||||
ylabel('DC value of the Compliance');
|
||||
xlabel('Optimal Damping Ratio');
|
||||
legend('location', 'northwest');
|
||||
#+end_src
|
||||
|
||||
#+begin_src matlab :tangle no :exports results :results file replace
|
||||
exportFig('figs/opt_damp_vs_dc_comp.pdf', 'width', 'wide', 'height', 'normal');
|
||||
#+end_src
|
||||
|
||||
#+name: fig:opt_damp_vs_dc_comp
|
||||
#+caption:
|
||||
#+RESULTS:
|
||||
[[file:figs/opt_damp_vs_dc_comp.png]]
|
||||
|
||||
*** Comparison :ignore:
|
||||
#+begin_src matlab :exports none
|
||||
figure;
|
||||
yyaxis left
|
||||
plot(wis, opt_xi_wi, '-', 'DisplayName', '$\xi_{cl}$');
|
||||
set(gca, 'YScale', 'lin');
|
||||
ylim([0,1]);
|
||||
ylabel('Attainable Damping Ratio $\xi$');
|
||||
|
||||
yyaxis right
|
||||
hold on;
|
||||
plot(wis, opt_gain_wi, '-', 'DisplayName', '$g_{opt}$');
|
||||
plot(wis, wis*((w0/W)^2 - 1), '--', 'DisplayName', '$g_{max}$');
|
||||
set(gca, 'YScale', 'lin');
|
||||
ylim([0,10]);
|
||||
ylabel('Controller gain $g$');
|
||||
|
||||
xlabel('$\omega_i/\omega_0$');
|
||||
set(gca, 'XScale', 'log');
|
||||
legend('location', 'northeast');
|
||||
#+end_src
|
||||
|
||||
#+begin_src matlab :exports none
|
||||
figure;
|
||||
yyaxis left
|
||||
plot(alphas, opt_xi_kp, '-', 'DisplayName', '$\xi_{cl}$');
|
||||
set(gca, 'YScale', 'lin');
|
||||
ylim([0,1]);
|
||||
ylabel('Attainable Damping Ratio $\xi$');
|
||||
|
||||
yyaxis right
|
||||
hold on;
|
||||
plot(alphas, opt_gain_kp, '-', 'DisplayName', '$g_{opt}$');
|
||||
set(gca, 'YScale', 'lin');
|
||||
ylim([0,2.5]);
|
||||
ylabel('Controller gain $g$');
|
||||
|
||||
xlabel('$\alpha$');
|
||||
set(gca, 'XScale', 'log');
|
||||
legend('location', 'northeast');
|
||||
#+end_src
|
||||
|
||||
#+begin_src matlab :tangle no :exports results :results file replace
|
||||
exportFig('figs/mod_iff_damping_kp.pdf', 'width', 'normal', 'height', 'normal');
|
||||
#+end_src
|
||||
|
||||
#+name: fig:mod_iff_damping_kp
|
||||
#+caption:
|
||||
#+RESULTS:
|
||||
[[file:figs/mod_iff_damping_kp.png]]
|
||||
|
||||
#+begin_src matlab :tangle no :exports none :results none
|
||||
exportFig('figs-paper/mod_iff_damping_kp.pdf', 'width', 'half', 'height', '650', 'png', false, 'pdf', false, 'svg', true);
|
||||
#+end_src
|
||||
|
||||
* Notations
|
||||
<<sec:notations>>
|
||||
|
||||
@@ -2583,7 +2200,6 @@ exportFig('figs/opt_damp_vs_dc_comp.pdf', 'width', 'wide', 'height', 'normal');
|
||||
| Actuator Force | $\bm{F}, F_u, F_v$ | =F= =Fu= =Fv= | N |
|
||||
| Force Sensor signal | $\bm{f}, f_u, f_v$ | =f= =fu= =fv= | N |
|
||||
| Relative Displacement | $\bm{d}, d_u, d_v$ | =d= =du= =dv= | m |
|
||||
| Relative Velocity | $\bm{v}, v_u, v_v$ | =v= =vu= =vv= | m/s |
|
||||
| Resonance freq. when $\Omega = 0$ | $\omega_0$ | =w0= | rad/s |
|
||||
| Rotation Speed | $\Omega = \dot{\theta}$ | =W= | rad/s |
|
||||
| Low Pass Filter corner frequency | $\omega_i$ | =wi= | rad/s |
|
||||
@@ -2594,13 +2210,6 @@ exportFig('figs/opt_damp_vs_dc_comp.pdf', 'width', 'wide', 'height', 'normal');
|
||||
| Complex number | $j$ | =j= | |
|
||||
| Frequency | $\omega$ | =w= | [rad/s] |
|
||||
|
||||
| | Mathematical Notation | Matlab | Unit |
|
||||
|----------------+------------------------------------------------+--------+---------|
|
||||
| IFF Plant | $\bm{G}_\text{IFF}(s) = \frac{\bm{f}}{\bm{F}}$ | =Giff= | N/N |
|
||||
| DVF Plant | $\bm{G}_\text{DVF}(s) = \frac{\bm{v}}{\bm{F}}$ | =Gdvf= | (m/s)/N |
|
||||
| IFF Controller | $\bm{K}_\text{IFF}(s)$ | =Kiff= | |
|
||||
| DVF Controller | $\bm{K}_\text{DVF}(s)$ | =Kdvf= | |
|
||||
|
||||
* Functions :noexport:
|
||||
** Sort Poles for the Root Locus
|
||||
:PROPERTIES:
|
||||
|
||||
@@ -43,24 +43,25 @@ for p_i = 1:size(p_ws, 1)
|
||||
end
|
||||
plot(Ws, zeros(size(Ws)), 'k--')
|
||||
hold off;
|
||||
xlabel('Rotation Frequency [rad/s]'); ylabel('Real Part');
|
||||
xlabel('Rotational Speed [rad/s]'); ylabel('Real Part');
|
||||
|
||||
ax2 = subplot(1,2,2);
|
||||
hold on;
|
||||
for p_i = 1:size(p_ws, 1)
|
||||
plot(Ws, imag(p_ws(p_i, :)), 'k-')
|
||||
plot(Ws, -imag(p_ws(p_i, :)), 'k-')
|
||||
plot(Ws, imag(p_ws(p_i, :)), 'k-')
|
||||
end
|
||||
hold off;
|
||||
xlabel('Rotation Frequency [rad/s]'); ylabel('Imaginary Part');
|
||||
xlabel('Rotational Speed [rad/s]'); ylabel('Imaginary Part');
|
||||
|
||||
% Simscape Model
|
||||
% Define the rotating speed for the Simscape Model.
|
||||
% In order to validate all the equations of motion, a Simscape model of the same system has been developed.
|
||||
% The dynamics of the system can be identified from the Simscape model and compare with the analytical model.
|
||||
|
||||
% The rotating speed for the Simscape Model is defined.
|
||||
|
||||
W = 0.1; % Rotation Speed [rad/s]
|
||||
|
||||
Kiff = tf(zeros(2));
|
||||
Kdvf = tf(zeros(2));
|
||||
|
||||
kp = 0; % Parallel Stiffness [N/m]
|
||||
cp = 0; % Parallel Damping [N/(m/s)]
|
||||
@@ -206,7 +207,7 @@ title('$d_u/F_v$, $d_v/F_u$');
|
||||
ax4 = subplot(2, 2, 4);
|
||||
hold on;
|
||||
for W_i = 1:length(Ws)
|
||||
plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{W_i}(1,1), freqs))))
|
||||
plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{W_i}(2,1), freqs))))
|
||||
end
|
||||
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
||||
xlabel('Frequency [rad/s]'); ylabel('Phase [deg]');
|
||||
|
||||
@@ -17,13 +17,12 @@ w0 = sqrt(k/m); % [rad/s]
|
||||
kp = 0; % [N/m]
|
||||
cp = 0; % [N/(m/s)]
|
||||
|
||||
% Simscape Model
|
||||
% Comparison of the Analytical Model and the Simscape Model
|
||||
% The rotation speed is set to $\Omega = 0.1 \omega_0$.
|
||||
|
||||
W = 0.1*w0; % [rad/s]
|
||||
|
||||
Kiff = tf(zeros(2));
|
||||
Kdvf = tf(zeros(2));
|
||||
|
||||
open('rotating_frame.slx');
|
||||
|
||||
@@ -45,7 +44,8 @@ Giff = linearize(mdl, io, 0);
|
||||
Giff.InputName = {'Fu', 'Fv'};
|
||||
Giff.OutputName = {'fu', 'fv'};
|
||||
|
||||
% Comparison of the Analytical Model and the Simscape Model
|
||||
|
||||
|
||||
% The same transfer function from $[F_u, F_v]$ to $[f_u, f_v]$ is written down from the analytical model.
|
||||
|
||||
Giff_th = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ...
|
||||
@@ -56,7 +56,6 @@ Giff_th = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2)
|
||||
|
||||
% The two are compared in Figure [[fig:plant_iff_comp_simscape_analytical]] and found to perfectly match.
|
||||
|
||||
|
||||
freqs = logspace(-1, 1, 1000);
|
||||
|
||||
figure;
|
||||
|
||||
@@ -206,41 +206,36 @@ xlabel('Real Part'); ylabel('Imaginary Part');
|
||||
|
||||
% 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.
|
||||
|
||||
|
||||
% 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 [[fig:mod_iff_damping_wi]]).
|
||||
|
||||
wis = logspace(-2, 1, 31)*w0; % [rad/s]
|
||||
|
||||
opt_zeta = zeros(1, length(wis)); % Optimal simultaneous damping
|
||||
wis = logspace(-2, 1, 100)*w0; % [rad/s]
|
||||
|
||||
opt_xi = zeros(1, length(wis)); % Optimal simultaneous damping
|
||||
opt_gain = zeros(1, length(wis)); % Corresponding optimal gain
|
||||
|
||||
for wi_i = 1:length(wis)
|
||||
wi = wis(wi_i);
|
||||
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_xi(wi_i) = 1/xi_opt;
|
||||
opt_gain(wi_i) = g_opt;
|
||||
end
|
||||
|
||||
figure;
|
||||
yyaxis left
|
||||
plot(wis, opt_zeta, '-o', 'DisplayName', '$\xi_{cl}$');
|
||||
plot(wis, opt_xi, '-', 'DisplayName', '$\xi_{cl}$');
|
||||
set(gca, 'YScale', 'lin');
|
||||
ylim([0,1]);
|
||||
ylabel('Attainable Damping Ratio $\xi$');
|
||||
|
||||
yyaxis right
|
||||
hold on;
|
||||
plot(wis, opt_gain, '-x', 'DisplayName', '$g_{opt}$');
|
||||
plot(wis, opt_gain, '-', 'DisplayName', '$g_{opt}$');
|
||||
plot(wis, wis*((w0/W)^2 - 1), '--', 'DisplayName', '$g_{max}$');
|
||||
set(gca, 'YScale', 'lin');
|
||||
ylim([0,10]);
|
||||
|
||||
@@ -1,61 +1,339 @@
|
||||
% Attainable Damping as a function of $k_p$
|
||||
%% Clear Workspace and Close figures
|
||||
clear; close all; clc;
|
||||
|
||||
tic;
|
||||
alphas = logspace(-2, 0, 10);
|
||||
gains = linspace(0.5, 2.5, 100);
|
||||
%% Intialize Laplace variable
|
||||
s = zpk('s');
|
||||
|
||||
opt_zeta = zeros(1, length(alphas)); % Optimal simultaneous damping
|
||||
opt_gain = zeros(1, length(alphas)); % Corresponding optimal gain
|
||||
% Plant Parameters
|
||||
% Let's define initial values for the model.
|
||||
|
||||
for alpha_i = 1:length(alphas)
|
||||
kp = alphas(alpha_i);
|
||||
k = 1 - alphas(alpha_i);
|
||||
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]
|
||||
|
||||
kp = 0; % [N/m]
|
||||
cp = 0; % [N/(m/s)]
|
||||
|
||||
% 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.
|
||||
|
||||
W = 0.1*w0; % [rad/s]
|
||||
|
||||
kp = 1.5*m*W^2;
|
||||
cp = 0;
|
||||
|
||||
Kiff = tf(zeros(2));
|
||||
|
||||
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, '/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'};
|
||||
|
||||
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'};
|
||||
|
||||
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');
|
||||
|
||||
% Effect of the parallel stiffness on the IFF plant
|
||||
% The rotation speed is set to $\Omega = 0.1 \omega_0$.
|
||||
|
||||
W = 0.1*w0; % [rad/s]
|
||||
|
||||
|
||||
|
||||
% And the IFF plant (transfer function from $[F_u, F_v]$ to $[f_u, f_v]$) is identified in three different cases:
|
||||
% - without parallel stiffness
|
||||
% - with a small parallel stiffness $k_p < m \Omega^2$
|
||||
% - with a large parallel stiffness $k_p > m \Omega^2$
|
||||
|
||||
% The results are shown in Figure [[fig:plant_iff_kp]].
|
||||
|
||||
% One can see that for $k_p > m \Omega^2$, the systems shows alternating complex conjugate poles and zeros.
|
||||
|
||||
|
||||
kp = 0;
|
||||
|
||||
w0p = sqrt((k + kp)/m);
|
||||
xip = c/(2*sqrt((k+kp)*m));
|
||||
|
||||
Giff = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ...
|
||||
(s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p));
|
||||
(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2];
|
||||
|
||||
kp = 0.5*m*W^2;
|
||||
k = 1 - kp;
|
||||
|
||||
w0p = sqrt((k + kp)/m);
|
||||
xip = c/(2*sqrt((k+kp)*m));
|
||||
|
||||
Giff_s = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ...
|
||||
(s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p));
|
||||
(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2];
|
||||
|
||||
kp = 1.5*m*W^2;
|
||||
k = 1 - kp;
|
||||
|
||||
w0p = sqrt((k + kp)/m);
|
||||
xip = c/(2*sqrt((k+kp)*m));
|
||||
|
||||
Giff_l = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ...
|
||||
(s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p));
|
||||
(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2];
|
||||
|
||||
freqs = logspace(-2, 1, 1000);
|
||||
|
||||
figure;
|
||||
|
||||
ax1 = subplot(2, 1, 1);
|
||||
hold on;
|
||||
plot(freqs, abs(squeeze(freqresp(Giff(1,1), freqs))), 'k-')
|
||||
plot(freqs, abs(squeeze(freqresp(Giff_s(1,1), freqs))), 'k--')
|
||||
plot(freqs, abs(squeeze(freqresp(Giff_l(1,1), freqs))), 'k:')
|
||||
hold off;
|
||||
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
||||
set(gca, 'XTickLabel',[]); ylabel('Magnitude [N/N]');
|
||||
ylim([1e-5, 2e1]);
|
||||
|
||||
ax2 = subplot(2, 1, 2);
|
||||
hold on;
|
||||
plot(freqs, 180/pi*angle(squeeze(freqresp(Giff(1,1), freqs))), 'k-', ...
|
||||
'DisplayName', '$k_p = 0$')
|
||||
plot(freqs, 180/pi*angle(squeeze(freqresp(Giff_s(1,1), freqs))), 'k--', ...
|
||||
'DisplayName', '$k_p < m\Omega^2$')
|
||||
plot(freqs, 180/pi*angle(squeeze(freqresp(Giff_l(1,1), freqs))), 'k:', ...
|
||||
'DisplayName', '$k_p > m\Omega^2$')
|
||||
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', 'southwest');
|
||||
|
||||
linkaxes([ax1,ax2],'x');
|
||||
xlim([freqs(1), freqs(end)]);
|
||||
|
||||
% 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}
|
||||
% K_{\text{IFF}} = \frac{g}{s} \begin{bmatrix}
|
||||
% 1 & 0 \\
|
||||
% 0 & 1
|
||||
% \end{bmatrix}
|
||||
% \end{equation}
|
||||
|
||||
% One can see that for $k_p > m \Omega^2$, the root locus stays in the left half of the complex plane and thus the control system is unconditionally stable.
|
||||
|
||||
% Thus, decentralized IFF controller with pure integrators can be used if:
|
||||
% \begin{equation}
|
||||
% k_{p} > m \Omega^2
|
||||
% \end{equation}
|
||||
|
||||
|
||||
figure;
|
||||
|
||||
gains = logspace(-2, 2, 100);
|
||||
|
||||
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 g = gains
|
||||
cl_poles = pole(feedback(Giff, (g/s)*eye(2)));
|
||||
set(gca,'ColorOrderIndex',1);
|
||||
plot(real(cl_poles), imag(cl_poles), '.', ...
|
||||
'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 g = gains
|
||||
cl_poles = pole(feedback(Giff_s, (g/s)*eye(2)));
|
||||
set(gca,'ColorOrderIndex',2);
|
||||
plot(real(cl_poles), imag(cl_poles), '.', ...
|
||||
'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 g = gains
|
||||
set(gca,'ColorOrderIndex',3);
|
||||
cl_poles = pole(feedback(Giff_l, (g/s)*eye(2)));
|
||||
plot(real(cl_poles), imag(cl_poles), '.', ...
|
||||
'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 g = gains
|
||||
cl_poles = pole(feedback(Giff, (g/s)*eye(2)));
|
||||
set(gca,'ColorOrderIndex',1);
|
||||
plot(real(cl_poles), imag(cl_poles), '.');
|
||||
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 g = gains
|
||||
cl_poles = pole(feedback(Giff_s, (g/s)*eye(2)));
|
||||
set(gca,'ColorOrderIndex',2);
|
||||
plot(real(cl_poles), imag(cl_poles), '.');
|
||||
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 g = gains
|
||||
set(gca,'ColorOrderIndex',3);
|
||||
cl_poles = pole(feedback(Giff_l, (g/s)*eye(2)));
|
||||
plot(real(cl_poles), imag(cl_poles), '.');
|
||||
end
|
||||
hold off;
|
||||
axis square;
|
||||
xlim([-0.04, 0.06]); ylim([0, 0.1]);
|
||||
|
||||
xlabel('Real Part'); ylabel('Imaginary Part');
|
||||
|
||||
% Effect of $k_p$ on the attainable damping
|
||||
% However, having large values of $k_p$ may decrease the attainable damping.
|
||||
|
||||
% To study the second point, Root Locus plots for the following values of $k_p$ are shown in Figure [[fig:root_locus_iff_kps]].
|
||||
|
||||
kps = [2, 20, 40]*m*W^2;
|
||||
|
||||
|
||||
|
||||
% It is shown that large values of $k_p$ decreases the attainable damping.
|
||||
|
||||
figure;
|
||||
|
||||
gains = logspace(-2, 4, 500);
|
||||
|
||||
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];
|
||||
(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', ...
|
||||
'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 g = gains
|
||||
Kiff = g/s*eye(2);
|
||||
|
||||
[w, zeta] = damp(minreal(feedback(Giff, Kiff)));
|
||||
|
||||
if min(zeta) > opt_zeta(alpha_i) && all(zeta > 0)
|
||||
opt_zeta(alpha_i) = min(zeta);
|
||||
opt_gain(alpha_i) = g;
|
||||
end
|
||||
Kiffa = (g/s)*eye(2);
|
||||
cl_poles = pole(feedback(Giff, Kiffa));
|
||||
set(gca,'ColorOrderIndex',kp_i);
|
||||
plot(real(cl_poles), imag(cl_poles), '.', ...
|
||||
'HandleVisibility', 'off');
|
||||
end
|
||||
end
|
||||
toc
|
||||
hold off;
|
||||
axis square;
|
||||
xlim([-1.2, 0.2]); ylim([0, 1.4]);
|
||||
|
||||
figure;
|
||||
yyaxis left
|
||||
plot(alphas, opt_zeta, '-o', 'DisplayName', '$\xi_{cl}$');
|
||||
set(gca, 'YScale', 'lin');
|
||||
ylim([0,1]);
|
||||
ylabel('Attainable Damping Ratio $\xi$');
|
||||
|
||||
yyaxis right
|
||||
hold on;
|
||||
plot(alphas, opt_gain, '-x', 'DisplayName', '$g_{opt}$');
|
||||
set(gca, 'YScale', 'lin');
|
||||
ylim([0,3]);
|
||||
ylabel('Controller gain $g$');
|
||||
|
||||
xlabel('$\alpha$');
|
||||
set(gca, 'XScale', 'log');
|
||||
legend('location', 'northeast');
|
||||
|
||||
|
||||
|
||||
% Alternative using fminserach
|
||||
xlabel('Real Part'); ylabel('Imaginary Part');
|
||||
legend('location', 'northwest');
|
||||
|
||||
alphas = logspace(-2, 0, 100);
|
||||
|
||||
opt_zeta = zeros(1, length(alphas)); % Optimal simultaneous damping
|
||||
opt_xi = zeros(1, length(alphas)); % Optimal simultaneous damping
|
||||
opt_gain = zeros(1, length(alphas)); % Corresponding optimal gain
|
||||
|
||||
Kiff = 1/s*eye(2);
|
||||
@@ -74,6 +352,24 @@ for alpha_i = 1:length(alphas)
|
||||
fun = @(g)computeSimultaneousDamping(g, Giff, Kiff);
|
||||
|
||||
[g_opt, xi_opt] = fminsearch(fun, 2);
|
||||
opt_zeta(alpha_i) = 1/xi_opt;
|
||||
opt_xi(alpha_i) = 1/xi_opt;
|
||||
opt_gain(alpha_i) = g_opt;
|
||||
end
|
||||
|
||||
figure;
|
||||
yyaxis left
|
||||
plot(alphas, opt_xi, '-', 'DisplayName', '$\xi_{cl}$');
|
||||
set(gca, 'YScale', 'lin');
|
||||
ylim([0,1]);
|
||||
ylabel('Attainable Damping Ratio $\xi$');
|
||||
|
||||
yyaxis right
|
||||
hold on;
|
||||
plot(alphas, opt_gain, '-', 'DisplayName', '$g_{opt}$');
|
||||
set(gca, 'YScale', 'lin');
|
||||
ylim([0,2.5]);
|
||||
ylabel('Controller gain $g$');
|
||||
|
||||
xlabel('$\alpha$');
|
||||
set(gca, 'XScale', 'log');
|
||||
legend('location', 'northeast');
|
||||
|
||||
227
matlab/matlab/s5_act_damp_comparison.m
Normal file
227
matlab/matlab/s5_act_damp_comparison.m
Normal file
@@ -0,0 +1,227 @@
|
||||
%% Clear Workspace and Close figures
|
||||
clear; close all; clc;
|
||||
|
||||
%% Intialize Laplace variable
|
||||
s = zpk('s');
|
||||
|
||||
% 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]
|
||||
|
||||
kp = 0; % [N/m]
|
||||
cp = 0; % [N/(m/s)]
|
||||
|
||||
|
||||
|
||||
% The rotating speed is set to $\Omega = 0.1 \omega_0$.
|
||||
|
||||
W = 0.1*w0;
|
||||
|
||||
% Root Locus
|
||||
% IFF with High Pass Filter
|
||||
|
||||
wi = 0.1*w0; % [rad/s]
|
||||
|
||||
Giff = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ...
|
||||
[(s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)) + (2*W*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*W*s/(w0^2) ; ...
|
||||
(2*xi*s/w0 + 1)*2*W*s/(w0^2), (s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))+ (2*W*s/(w0^2))^2];
|
||||
|
||||
|
||||
|
||||
% IFF With parallel Stiffness
|
||||
|
||||
kp = 5*m*W^2;
|
||||
k = k - kp;
|
||||
|
||||
w0p = sqrt((k + kp)/m);
|
||||
xip = c/(2*sqrt((k+kp)*m));
|
||||
|
||||
Giff_kp = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ...
|
||||
(s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p));
|
||||
(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2 ];
|
||||
|
||||
k = k + kp;
|
||||
|
||||
figure;
|
||||
|
||||
gains = logspace(-2, 2, 100);
|
||||
|
||||
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 g = gains
|
||||
Kiff = (g/(wi + s))*eye(2);
|
||||
cl_poles = pole(feedback(Giff, Kiff));
|
||||
set(gca,'ColorOrderIndex',1);
|
||||
plot(real(cl_poles), imag(cl_poles), '.', ...
|
||||
'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 g = gains
|
||||
Kiffa = (g/s)*eye(2);
|
||||
cl_poles = pole(feedback(Giff_kp, Kiffa));
|
||||
set(gca,'ColorOrderIndex',2);
|
||||
plot(real(cl_poles), imag(cl_poles), '.', ...
|
||||
'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');
|
||||
|
||||
% 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.
|
||||
|
||||
|
||||
fun = @(g)computeSimultaneousDamping(g, Giff, (1/(wi+s))*eye(2));
|
||||
|
||||
[opt_gain_iff, opt_xi_iff] = fminsearch(fun, 0.5*(w0^2/W^2 - 1)*wi);
|
||||
opt_xi_iff = 1/opt_xi_iff;
|
||||
|
||||
fun = @(g)computeSimultaneousDamping(g, Giff_kp, 1/s*eye(2));
|
||||
|
||||
[opt_gain_kp, opt_xi_kp] = fminsearch(fun, 2);
|
||||
opt_xi_kp = 1/opt_xi_kp;
|
||||
|
||||
% Passive Damping - Critical Damping
|
||||
% \begin{equation}
|
||||
% \xi = \frac{c}{2 \sqrt{km}}
|
||||
% \end{equation}
|
||||
|
||||
% Critical Damping corresponds to to $\xi = 1$, and thus:
|
||||
% \begin{equation}
|
||||
% c_{\text{crit}} = 2 \sqrt{km}
|
||||
% \end{equation}
|
||||
|
||||
|
||||
c_opt = 2*sqrt(k*m);
|
||||
|
||||
% Transmissibility And Compliance
|
||||
% <<sec:comp_transmissibilty>>
|
||||
|
||||
|
||||
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, '/fd'], 1, 'input'); io_i = io_i + 1;
|
||||
io(io_i) = linio([mdl, '/Meas'], 1, 'output'); io_i = io_i + 1;
|
||||
|
||||
% Open Loop :ignore:
|
||||
|
||||
Kiff = tf(zeros(2));
|
||||
|
||||
kp = 0;
|
||||
cp = 0;
|
||||
|
||||
G_ol = linearize(mdl, io, 0);
|
||||
|
||||
%% Input/Output definition
|
||||
G_ol.InputName = {'Dwx', 'Dwy', 'Fdx', 'Fdy'};
|
||||
G_ol.OutputName = {'Dx', 'Dy'};
|
||||
|
||||
% Passive Damping
|
||||
|
||||
kp = 0;
|
||||
cp = 0;
|
||||
|
||||
c_old = c;
|
||||
c = c_opt;
|
||||
|
||||
Kiff = tf(zeros(2));
|
||||
|
||||
G_pas = linearize(mdl, io, 0);
|
||||
|
||||
%% Input/Output definition
|
||||
G_pas.InputName = {'Dwx', 'Dwy', 'Fdx', 'Fdy'};
|
||||
G_pas.OutputName = {'Dx', 'Dy'};
|
||||
|
||||
c = c_old;
|
||||
|
||||
% Pseudo Integrator IFF :ignore:
|
||||
|
||||
kp = 0;
|
||||
cp = 0;
|
||||
|
||||
Kiff = opt_gain_iff/(wi + s)*tf(eye(2));
|
||||
|
||||
G_iff = linearize(mdl, io, 0);
|
||||
|
||||
%% Input/Output definition
|
||||
G_iff.InputName = {'Dwx', 'Dwy', 'Fdx', 'Fdy'};
|
||||
G_iff.OutputName = {'Dx', 'Dy'};
|
||||
|
||||
% IFF With parallel Stiffness :ignore:
|
||||
|
||||
kp = 5*m*W^2;
|
||||
cp = 0.01;
|
||||
|
||||
Kiff = opt_gain_kp/s*tf(eye(2));
|
||||
|
||||
G_kp = linearize(mdl, io, 0);
|
||||
|
||||
%% Input/Output definition
|
||||
G_kp.InputName = {'Dwx', 'Dwy', 'Fdx', 'Fdy'};
|
||||
G_kp.OutputName = {'Dx', 'Dy'};
|
||||
|
||||
% Transmissibility :ignore:
|
||||
|
||||
freqs = logspace(-2, 1, 1000);
|
||||
|
||||
figure;
|
||||
hold on;
|
||||
plot(freqs, abs(squeeze(freqresp(G_iff({'Dx'}, {'Dwx'}), freqs))), ...
|
||||
'DisplayName', 'IFF + HPF')
|
||||
plot(freqs, abs(squeeze(freqresp(G_kp( {'Dx'}, {'Dwx'}), freqs))), ...
|
||||
'DisplayName', 'IFF + $k_p$')
|
||||
plot(freqs, abs(squeeze(freqresp(G_pas({'Dx'}, {'Dwx'}), freqs))), ...
|
||||
'DisplayName', 'Passive')
|
||||
plot(freqs, abs(squeeze(freqresp(G_ol( {'Dx'}, {'Dwx'}), freqs))), 'k-', ...
|
||||
'DisplayName', 'Open-Loop')
|
||||
hold off;
|
||||
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
||||
ylim([1e-2, 3e1]);
|
||||
xlabel('Frequency [rad/s]'); ylabel('Transmissibility [m/m]');
|
||||
legend('location', 'southwest');
|
||||
|
||||
% Compliance :ignore:
|
||||
|
||||
freqs = logspace(-2, 1, 1000);
|
||||
|
||||
figure;
|
||||
hold on;
|
||||
plot(freqs, abs(squeeze(freqresp(G_iff({'Dx'}, {'Fdx'}), freqs))), ...
|
||||
'DisplayName', 'IFF + HPF')
|
||||
plot(freqs, abs(squeeze(freqresp(G_kp( {'Dx'}, {'Fdx'}), freqs))), ...
|
||||
'DisplayName', 'IFF + $k_p$')
|
||||
plot(freqs, abs(squeeze(freqresp(G_pas({'Dx'}, {'Fdx'}), freqs))), ...
|
||||
'DisplayName', 'Passive')
|
||||
plot(freqs, abs(squeeze(freqresp(G_ol( {'Dx'}, {'Fdx'}), freqs))), 'k-', ...
|
||||
'DisplayName', 'Open-Loop')
|
||||
hold off;
|
||||
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
||||
ylim([1e-2, 3e1]);
|
||||
xlabel('Frequency [rad/s]'); ylabel('Compliance [m/N]');
|
||||
legend('location', 'southwest');
|
||||
@@ -1,159 +0,0 @@
|
||||
%% Clear Workspace and Close figures
|
||||
clear; close all; clc;
|
||||
|
||||
%% Intialize Laplace variable
|
||||
s = zpk('s');
|
||||
|
||||
% 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]
|
||||
|
||||
kp = 0; % [N/m]
|
||||
cp = 0; % [N/(m/s)]
|
||||
|
||||
% Comparison of the Analytical Model and the Simscape Model
|
||||
% The rotating speed is set to $\Omega = 0.1 \omega_0$.
|
||||
|
||||
W = 0.1*w0;
|
||||
|
||||
Kiff = tf(zeros(2));
|
||||
Kdvf = tf(zeros(2));
|
||||
|
||||
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 [[fig:plant_iff_comp_simscape_analytical]] and found to perfectly match.
|
||||
|
||||
|
||||
freqs = logspace(-1, 1, 1000);
|
||||
|
||||
figure;
|
||||
ax1 = subplot(2, 2, 1);
|
||||
hold on;
|
||||
plot(freqs, abs(squeeze(freqresp(Gdvf(1,1), freqs))), '-')
|
||||
plot(freqs, abs(squeeze(freqresp(Gdvf_th(1,1), freqs))), '--')
|
||||
hold off;
|
||||
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
||||
set(gca, 'XTickLabel',[]); ylabel('Magnitude [$\frac{m/s}{N}$]');
|
||||
title('$v_u/F_u$, $v_v/F_v$');
|
||||
|
||||
ax3 = subplot(2, 2, 3);
|
||||
hold on;
|
||||
plot(freqs, 180/pi*angle(squeeze(freqresp(Gdvf(1,1), freqs))), '-')
|
||||
plot(freqs, 180/pi*angle(squeeze(freqresp(Gdvf_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(Gdvf(1,2), freqs))), '-')
|
||||
plot(freqs, abs(squeeze(freqresp(Gdvf_th(1,2), freqs))), '--')
|
||||
hold off;
|
||||
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
||||
set(gca, 'XTickLabel',[]); ylabel('Magnitude [$\frac{m/s}{N}$]');
|
||||
title('$v_u/F_v$, $v_v/F_u$');
|
||||
|
||||
ax4 = subplot(2, 2, 4);
|
||||
hold on;
|
||||
plot(freqs, 180/pi*angle(squeeze(freqresp(Gdvf(1,2), freqs))), '-')
|
||||
plot(freqs, 180/pi*angle(squeeze(freqresp(Gdvf_th(1,2), 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;
|
||||
|
||||
linkaxes([ax1,ax2,ax3,ax4],'x');
|
||||
xlim([freqs(1), freqs(end)]);
|
||||
|
||||
linkaxes([ax1,ax2],'y');
|
||||
|
||||
% 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 [[fig:root_locus_dvf]].
|
||||
|
||||
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.
|
||||
|
||||
gains = logspace(-2, 1, 100);
|
||||
|
||||
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');
|
||||
|
||||
for g = gains
|
||||
set(gca,'ColorOrderIndex',W_i);
|
||||
cl_poles = pole(feedback(Gdvf, g*eye(2)));
|
||||
|
||||
plot(real(cl_poles), imag(cl_poles), '.', ...
|
||||
'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');
|
||||
@@ -1,364 +0,0 @@
|
||||
%% Clear Workspace and Close figures
|
||||
clear; close all; clc;
|
||||
|
||||
%% Intialize Laplace variable
|
||||
s = zpk('s');
|
||||
|
||||
% 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]
|
||||
|
||||
kp = 0; % [N/m]
|
||||
cp = 0; % [N/(m/s)]
|
||||
|
||||
|
||||
|
||||
% The rotating speed is set to $\Omega = 0.1 \omega_0$.
|
||||
|
||||
W = 0.1*w0;
|
||||
|
||||
% Root Locus
|
||||
% IFF with High Pass Filter
|
||||
|
||||
wi = 0.1*w0; % [rad/s]
|
||||
|
||||
Giff = 1/(((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))^2 + (2*W*s/(w0^2))^2) * ...
|
||||
[(s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2)) + (2*W*s/(w0^2))^2, - (2*xi*s/w0 + 1)*2*W*s/(w0^2) ; ...
|
||||
(2*xi*s/w0 + 1)*2*W*s/(w0^2), (s^2/w0^2 - W^2/w0^2)*((s^2)/(w0^2) + 2*xi*s/w0 + 1 - (W^2)/(w0^2))+ (2*W*s/(w0^2))^2];
|
||||
|
||||
|
||||
|
||||
% IFF With parallel Stiffness
|
||||
|
||||
kp = 5*m*W^2;
|
||||
k = k - kp;
|
||||
|
||||
w0p = sqrt((k + kp)/m);
|
||||
xip = c/(2*sqrt((k+kp)*m));
|
||||
|
||||
Giff_kp = 1/( (s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2)^2 + (2*(s/w0p)*(W/w0p))^2 ) * [ ...
|
||||
(s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2, -(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p));
|
||||
(2*xip*s/w0p + k/(k + kp))*(2*(s/w0p)*(W/w0p)), (s^2/w0p^2 + kp/(k + kp) - W^2/w0p^2)*(s^2/w0p^2 + 2*xip*s/w0p + 1 - W^2/w0p^2) + (2*(s/w0p)*(W/w0p))^2 ];
|
||||
|
||||
k = k + kp;
|
||||
|
||||
|
||||
|
||||
% 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;
|
||||
|
||||
gains = logspace(-2, 2, 100);
|
||||
|
||||
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 g = gains
|
||||
Kiff = (g/(wi + s))*eye(2);
|
||||
cl_poles = pole(feedback(Giff, Kiff));
|
||||
set(gca,'ColorOrderIndex',1);
|
||||
plot(real(cl_poles), imag(cl_poles), '.', ...
|
||||
'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 g = gains
|
||||
Kiffa = (g/s)*eye(2);
|
||||
cl_poles = pole(feedback(Giff_kp, Kiffa));
|
||||
set(gca,'ColorOrderIndex',2);
|
||||
plot(real(cl_poles), imag(cl_poles), '.', ...
|
||||
'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 g = gains
|
||||
Kdvf = g*eye(2);
|
||||
cl_poles = pole(feedback(Gdvf, Kdvf));
|
||||
set(gca,'ColorOrderIndex',3);
|
||||
plot(real(cl_poles), imag(cl_poles), '.', ...
|
||||
'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');
|
||||
|
||||
% 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.
|
||||
|
||||
|
||||
%% IFF with pseudo integrators
|
||||
gains = linspace(0, (w0^2/W^2 - 1)*wi, 100);
|
||||
opt_zeta_iff = 0;
|
||||
opt_gain_iff = 0;
|
||||
|
||||
for g = gains
|
||||
Kiff = (g/(wi+s))*eye(2);
|
||||
|
||||
[w, zeta] = damp(minreal(feedback(Giff, Kiff)));
|
||||
|
||||
if min(zeta) > opt_zeta_iff && all(zeta > 0)
|
||||
opt_zeta_iff = min(zeta);
|
||||
opt_gain_iff = g;
|
||||
end
|
||||
end
|
||||
|
||||
%% IFF with Parallel Stiffness
|
||||
gains = logspace(-2, 4, 100);
|
||||
opt_zeta_kp = 0;
|
||||
opt_gain_kp = 0;
|
||||
|
||||
for g = gains
|
||||
Kiff = g/s*eye(2);
|
||||
|
||||
[w, zeta] = damp(minreal(feedback(Giff_kp, Kiff)));
|
||||
|
||||
if min(zeta) > opt_zeta_kp && all(zeta > 0)
|
||||
opt_zeta_kp = min(zeta);
|
||||
opt_gain_kp = g;
|
||||
end
|
||||
end
|
||||
|
||||
%% Direct Velocity Feedback
|
||||
gains = logspace(0, 2, 100);
|
||||
opt_zeta_dvf = 0;
|
||||
opt_gain_dvf = 0;
|
||||
|
||||
for g = gains
|
||||
Kdvf = g*eye(2);
|
||||
|
||||
[w, zeta] = damp(minreal(feedback(Gdvf, Kdvf)));
|
||||
|
||||
if min(zeta) > opt_zeta_dvf && all(zeta > 0) && min(zeta) < 0.85
|
||||
opt_zeta_dvf = min(zeta);
|
||||
opt_gain_dvf = g;
|
||||
end
|
||||
end
|
||||
|
||||
% Transmissibility
|
||||
% <<sec:comp_transmissibilty>>
|
||||
|
||||
|
||||
open('rotating_frame.slx');
|
||||
|
||||
% Open Loop :ignore:
|
||||
|
||||
Kdvf = tf(zeros(2));
|
||||
Kiff = tf(zeros(2));
|
||||
|
||||
kp = 0;
|
||||
cp = 0;
|
||||
|
||||
%% Name of the Simulink File
|
||||
mdl = 'rotating_frame';
|
||||
|
||||
%% Input/Output definition
|
||||
clear io; io_i = 1;
|
||||
io(io_i) = linio([mdl, '/dw'], 1, 'input'); io_i = io_i + 1;
|
||||
io(io_i) = linio([mdl, '/Meas'], 1, 'output'); io_i = io_i + 1;
|
||||
|
||||
Tol = linearize(mdl, io, 0);
|
||||
|
||||
%% Input/Output definition
|
||||
Tol.InputName = {'Dwx', 'Dwy'};
|
||||
Tol.OutputName = {'Dx', 'Dy'};
|
||||
|
||||
% Pseudo Integrator IFF :ignore:
|
||||
|
||||
kp = 0;
|
||||
cp = 0;
|
||||
|
||||
Kdvf = tf(zeros(2));
|
||||
|
||||
Kiff = opt_gain_iff/(wi + s)*tf(eye(2));
|
||||
|
||||
%% Name of the Simulink File
|
||||
mdl = 'rotating_frame';
|
||||
|
||||
%% Input/Output definition
|
||||
clear io; io_i = 1;
|
||||
io(io_i) = linio([mdl, '/dw'], 1, 'input'); io_i = io_i + 1;
|
||||
io(io_i) = linio([mdl, '/Meas'], 1, 'output'); io_i = io_i + 1;
|
||||
|
||||
Tiff = linearize(mdl, io, 0);
|
||||
|
||||
%% Input/Output definition
|
||||
Tiff.InputName = {'Dwx', 'Dwy'};
|
||||
Tiff.OutputName = {'Dx', 'Dy'};
|
||||
|
||||
% IFF With parallel Stiffness :ignore:
|
||||
|
||||
kp = 5*m*W^2;
|
||||
cp = 0.01;
|
||||
|
||||
Kiff = opt_gain_kp/s*tf(eye(2));
|
||||
|
||||
Kdvf = tf(zeros(2));
|
||||
|
||||
%% Name of the Simulink File
|
||||
mdl = 'rotating_frame';
|
||||
|
||||
%% Input/Output definition
|
||||
clear io; io_i = 1;
|
||||
io(io_i) = linio([mdl, '/dw'], 1, 'input'); io_i = io_i + 1;
|
||||
io(io_i) = linio([mdl, '/Meas'], 1, 'output'); io_i = io_i + 1;
|
||||
|
||||
Tiff_kp = linearize(mdl, io, 0);
|
||||
|
||||
%% Input/Output definition
|
||||
Tiff_kp.InputName = {'Dwx', 'Dwy'};
|
||||
Tiff_kp.OutputName = {'Dx', 'Dy'};
|
||||
|
||||
% DVF :ignore:
|
||||
|
||||
kp = 0;
|
||||
cp = 0;
|
||||
|
||||
Kiff = tf(zeros(2));
|
||||
|
||||
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'};
|
||||
|
||||
% Transmissibility :ignore:
|
||||
|
||||
freqs = logspace(-2, 1, 1000);
|
||||
|
||||
figure;
|
||||
hold on;
|
||||
plot(freqs, abs(squeeze(freqresp(Tiff(1,1), freqs))), ...
|
||||
'DisplayName', 'IFF + HPF')
|
||||
plot(freqs, abs(squeeze(freqresp(Tiff_kp(1,1), freqs))), ...
|
||||
'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', 'Open-Loop')
|
||||
hold off;
|
||||
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
||||
xlabel('Frequency [rad/s]'); ylabel('Transmissibility [m/m]');
|
||||
legend('location', 'southwest');
|
||||
|
||||
% Open Loop :ignore:
|
||||
|
||||
Kdvf = tf(zeros(2));
|
||||
Kiff = tf(zeros(2));
|
||||
|
||||
kp = 0;
|
||||
cp = 0;
|
||||
|
||||
%% Name of the Simulink File
|
||||
mdl = 'rotating_frame';
|
||||
|
||||
%% Input/Output definition
|
||||
clear io; io_i = 1;
|
||||
io(io_i) = linio([mdl, '/fd'], 1, 'input'); io_i = io_i + 1;
|
||||
io(io_i) = linio([mdl, '/Meas'], 1, 'output'); io_i = io_i + 1;
|
||||
|
||||
Col = linearize(mdl, io, 0);
|
||||
|
||||
%% Input/Output definition
|
||||
Col.InputName = {'Fdx', 'Fdy'};
|
||||
Col.OutputName = {'Dx', 'Dy'};
|
||||
|
||||
% Pseudo Integrator IFF :ignore:
|
||||
|
||||
kp = 0;
|
||||
cp = 0;
|
||||
|
||||
Kdvf = tf(zeros(2));
|
||||
|
||||
Kiff = opt_gain_iff/(wi + s)*tf(eye(2));
|
||||
|
||||
Ciff = linearize(mdl, io, 0);
|
||||
|
||||
%% Input/Output definition
|
||||
Ciff.InputName = {'Fdx', 'Fdy'};
|
||||
Ciff.OutputName = {'Dx', 'Dy'};
|
||||
|
||||
% IFF With parallel Stiffness :ignore:
|
||||
|
||||
kp = 5*m*W^2;
|
||||
cp = 0.01;
|
||||
|
||||
Kiff = opt_gain_kp/s*tf(eye(2));
|
||||
|
||||
Kdvf = tf(zeros(2));
|
||||
|
||||
Ciff_kp = linearize(mdl, io, 0);
|
||||
|
||||
%% Input/Output definition
|
||||
Ciff_kp.InputName = {'Fdx', 'Fdy'};
|
||||
Ciff_kp.OutputName = {'Dx', 'Dy'};
|
||||
|
||||
% DVF :ignore:
|
||||
|
||||
kp = 0;
|
||||
cp = 0;
|
||||
|
||||
Kiff = tf(zeros(2));
|
||||
|
||||
Kdvf = opt_gain_kp*tf(eye(2));
|
||||
|
||||
Cdvf = linearize(mdl, io, 0);
|
||||
|
||||
%% Input/Output definition
|
||||
Cdvf.InputName = {'Fdx', 'Fdy'};
|
||||
Cdvf.OutputName = {'Dx', 'Dy'};
|
||||
|
||||
% Compliance :ignore:
|
||||
|
||||
freqs = logspace(-2, 1, 1000);
|
||||
|
||||
figure;
|
||||
hold on;
|
||||
plot(freqs, abs(squeeze(freqresp(Ciff(1,1), freqs))), ...
|
||||
'DisplayName', 'IFF + HPF')
|
||||
plot(freqs, abs(squeeze(freqresp(Ciff_kp(1,1), freqs))), ...
|
||||
'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', 'Open-Loop')
|
||||
hold off;
|
||||
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
||||
xlabel('Frequency [rad/s]'); ylabel('Compliance [m/N]');
|
||||
legend('location', 'southwest');
|
||||
@@ -64,4 +64,13 @@
|
||||
H-infinity control},
|
||||
year = 1992,
|
||||
keywords = {},
|
||||
}
|
||||
}
|
||||
|
||||
@inproceedings{dehaeze20_activ_dampin_rotat_platf_integ_force_feedb,
|
||||
author = {Dehaeze, T. and Collette, C.},
|
||||
title = {Active Damping of Rotating Platforms using Integral Force
|
||||
Feedback},
|
||||
booktitle = {Proceedings of the International Conference on Modal
|
||||
Analysis Noise and Vibration Engineering (ISMA)},
|
||||
year = 2020,
|
||||
}
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
% Compute Damping
|
||||
% =computeSimultaneousDamping=
|
||||
|
||||
function [xi_min] = computeSimultaneousDamping(g, G, K)
|
||||
[w, xi] = damp(minreal(feedback(G, g*K)));
|
||||
|
||||
Reference in New Issue
Block a user