1758 lines
58 KiB
Org Mode
1758 lines
58 KiB
Org Mode
#+TITLE: Control in a rotating frame
|
|
:DRAWER:
|
|
#+STARTUP: overview
|
|
|
|
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="css/htmlize.css"/>
|
|
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="css/readtheorg.css"/>
|
|
#+HTML_HEAD: <script src="js/jquery.min.js"></script>
|
|
#+HTML_HEAD: <script src="js/bootstrap.min.js"></script>
|
|
#+HTML_HEAD: <script type="text/javascript" src="js/jquery.stickytableheaders.min.js"></script>
|
|
#+HTML_HEAD: <script type="text/javascript" src="js/readtheorg.js"></script>
|
|
|
|
#+LATEX_CLASS: cleanreport
|
|
#+LaTeX_CLASS_OPTIONS: [tocnp, secbreak, minted]
|
|
#+LaTeX_HEADER: \usepackage{svg}
|
|
#+LaTeX_HEADER: \newcommand{\authorFirstName}{Thomas}
|
|
#+LaTeX_HEADER: \newcommand{\authorLastName}{Dehaeze}
|
|
#+LaTeX_HEADER: \newcommand{\authorEmail}{dehaeze.thomas@gmail.com}
|
|
|
|
#+PROPERTY: header-args:matlab :session *MATLAB*
|
|
#+PROPERTY: header-args:matlab+ :comments org
|
|
#+PROPERTY: header-args:matlab+ :exports both
|
|
#+PROPERTY: header-args:matlab+ :eval no-export
|
|
#+PROPERTY: header-args:matlab+ :output-dir figs
|
|
#+PROPERTY: header-args:matlab+ :mkdirp yes
|
|
:END:
|
|
|
|
* Introduction
|
|
The objective of this note it to highlight some control problems that arises when controlling the position of an object using actuators that are rotating with respect to a fixed reference frame.
|
|
|
|
In section ref:sec:system, a simple system composed of a spindle and a translation stage is defined and the equations of motion are written.
|
|
The rotation induces some coupling between the actuators and their displacement, and modifies the dynamics of the system.
|
|
This is studied using the equations, and some numerical computations are used to compare the use of voice coil and piezoelectric actuators.
|
|
|
|
Then, in section ref:sec:control_strategies, two different control approach are compared where:
|
|
- the measurement is made in the fixed frame
|
|
- the measurement is made in the rotating frame
|
|
|
|
In section ref:sec:simscape, the analytical study will be validated using a multi body model of the studied system.
|
|
|
|
Finally, in section ref:sec:control, the control strategies are implemented using Simulink and Simscape and compared.
|
|
|
|
* System Description and Analysis
|
|
:PROPERTIES:
|
|
:HEADER-ARGS:matlab+: :tangle matlab/system_numerical_analysis.m
|
|
:END:
|
|
<<sec:system>>
|
|
** System description
|
|
The system consists of one 2 degree of freedom translation stage on top of a spindle (figure ref:fig:rotating_frame_2dof).
|
|
|
|
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 $u$ and $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.
|
|
|
|
#+name: fig:rotating_frame_2dof
|
|
#+caption: Schematic of the mecanical system
|
|
[[./figs/rotating_frame_2dof.png]]
|
|
|
|
In the following block diagram:
|
|
- $G$ is the transfer function from the forces applied in the actuators to the measurement
|
|
- $K$ is the controller to design
|
|
- $J$ is a Jacobian matrix usually used to change the reference frame
|
|
|
|
Indices $x$ and $y$ corresponds to signals in the fixed reference frame (along $\vec{i}_x$ and $\vec{i}_y$):
|
|
- $D_x$ is the measured position of the sample
|
|
- $r_x$ is the reference signal which corresponds to the wanted $D_x$
|
|
- $\epsilon_x$ is the position error
|
|
|
|
Indices $u$ and $v$ corresponds to signals in the rotating reference frame ($\vec{i}_u$ and $\vec{i}_v$):
|
|
- $F_u$ and $F_v$ are forces applied by the actuators
|
|
- $\epsilon_u$ and $\epsilon_v$ are position error of the sample along $\vec{i}_u$ and $\vec{i}_v$
|
|
|
|
** Equations
|
|
<<sec:equations>>
|
|
Based on the figure ref:fig:rotating_frame_2dof, we can write the equations of motion of the system.
|
|
|
|
Let's express the kinetic energy $T$ and the potential energy $V$ of the mass $m$:
|
|
#+name: eq:energy_inertial_frame
|
|
\begin{align}
|
|
T & = \frac{1}{2} m \left( \dot{x}^2 + \dot{y}^2 \right) \\
|
|
V & = \frac{1}{2} k \left( x^2 + y^2 \right)
|
|
\end{align}
|
|
|
|
The Lagrangian is the kinetic energy minus the potential energy.
|
|
#+name: eq:lagrangian_inertial_frame
|
|
\begin{equation}
|
|
L = T-V = \frac{1}{2} m \left( \dot{x}^2 + \dot{y}^2 \right) - \frac{1}{2} k \left( x^2 + y^2 \right)
|
|
\end{equation}
|
|
|
|
|
|
The partial derivatives of the Lagrangian with respect to the variables $(x, y)$ are:
|
|
#+name: eq:inertial_frame_deriv
|
|
\begin{align*}
|
|
\frac{\partial L}{\partial x} & = -kx \\
|
|
\frac{\partial L}{\partial y} & = -ky \\
|
|
\frac{d}{dt}\frac{\partial L}{\partial \dot{x}} & = m\ddot{x} \\
|
|
\frac{d}{dt}\frac{\partial L}{\partial \dot{y}} & = m\ddot{y}
|
|
\end{align*}
|
|
|
|
The external forces applied to the mass are:
|
|
\begin{align*}
|
|
F_{\text{ext}, x} &= F_u \cos{\theta} - F_v \sin{\theta}\\
|
|
F_{\text{ext}, y} &= F_u \sin{\theta} + F_v \cos{\theta}
|
|
\end{align*}
|
|
|
|
By appling the Lagrangian equations, we obtain:
|
|
\begin{align}
|
|
m\ddot{x} + kx = F_u \cos{\theta} - F_v \sin{\theta}\\
|
|
m\ddot{y} + ky = F_u \sin{\theta} + F_v \cos{\theta}
|
|
\end{align}
|
|
|
|
We then change coordinates from $(x, y)$ to $(d_x, d_y, \theta)$.
|
|
\begin{align*}
|
|
x & = d_u \cos{\theta} - d_v \sin{\theta}\\
|
|
y & = d_u \sin{\theta} + d_v \cos{\theta}
|
|
\end{align*}
|
|
|
|
We obtain:
|
|
\begin{align*}
|
|
\ddot{x} & = \ddot{d_u} \cos{\theta} - 2\dot{d_u}\dot{\theta}\sin{\theta} - d_u\ddot{\theta}\sin{\theta} - d_u\dot{\theta}^2 \cos{\theta}
|
|
- \ddot{d_v} \sin{\theta} - 2\dot{d_v}\dot{\theta}\cos{\theta} - d_v\ddot{\theta}\cos{\theta} + d_v\dot{\theta}^2 \sin{\theta} \\
|
|
\ddot{y} & = \ddot{d_u} \sin{\theta} + 2\dot{d_u}\dot{\theta}\cos{\theta} + d_u\ddot{\theta}\cos{\theta} - d_u\dot{\theta}^2 \sin{\theta}
|
|
+ \ddot{d_v} \cos{\theta} - 2\dot{d_v}\dot{\theta}\sin{\theta} - d_v\ddot{\theta}\sin{\theta} - d_v\dot{\theta}^2 \cos{\theta} \\
|
|
\end{align*}
|
|
|
|
By injecting the previous result into the Lagrangian equation, we obtain:
|
|
\begin{align*}
|
|
m \ddot{d_u} \cos{\theta} - 2m\dot{d_u}\dot{\theta}\sin{\theta} - m d_u\ddot{\theta}\sin{\theta} - m d_u\dot{\theta}^2 \cos{\theta}
|
|
-m \ddot{d_v} \sin{\theta} - 2m\dot{d_v}\dot{\theta}\cos{\theta} - m d_v\ddot{\theta}\cos{\theta} + m d_v\dot{\theta}^2 \sin{\theta}
|
|
+ k d_u \cos{\theta} - k d_v \sin{\theta} = F_u \cos{\theta} - F_v \sin{\theta} \\
|
|
m \ddot{d_u} \sin{\theta} + 2m\dot{d_u}\dot{\theta}\cos{\theta} + m d_u\ddot{\theta}\cos{\theta} - m d_u\dot{\theta}^2 \sin{\theta}
|
|
+ m \ddot{d_v} \cos{\theta} - 2m\dot{d_v}\dot{\theta}\sin{\theta} - m d_v\ddot{\theta}\sin{\theta} - m d_v\dot{\theta}^2 \cos{\theta}
|
|
+ k d_u \sin{\theta} + k d_v \cos{\theta} = F_u \sin{\theta} + F_v \cos{\theta}
|
|
\end{align*}
|
|
|
|
Which is equivalent to:
|
|
\begin{align*}
|
|
m \ddot{d_u} - 2m\dot{d_u}\dot{\theta}\frac{\sin{\theta}}{\cos{\theta}} - m d_u\ddot{\theta}\frac{\sin{\theta}}{\cos{\theta}} - m d_u\dot{\theta}^2
|
|
-m \ddot{d_v} \frac{\sin{\theta}}{\cos{\theta}} - 2m\dot{d_v}\dot{\theta} - m d_v\ddot{\theta} + m d_v\dot{\theta}^2 \frac{\sin{\theta}}{\cos{\theta}}
|
|
+ k d_u - k d_v \frac{\sin{\theta}}{\cos{\theta}} = F_u - F_v \frac{\sin{\theta}}{\cos{\theta}} \\
|
|
m \ddot{d_u} + 2m\dot{d_u}\dot{\theta}\frac{\cos{\theta}}{\sin{\theta}} + m d_u\ddot{\theta}\frac{\cos{\theta}}{\sin{\theta}} - m d_u\dot{\theta}^2
|
|
+ m \ddot{d_v} \frac{\cos{\theta}}{\sin{\theta}} - 2m\dot{d_v}\dot{\theta} - m d_v\ddot{\theta} - m d_v\dot{\theta}^2 \frac{\cos{\theta}}{\sin{\theta}}
|
|
+ k d_u + k d_v \frac{\cos{\theta}}{\sin{\theta}} = F_u + F_v \frac{\cos{\theta}}{\sin{\theta}}
|
|
\end{align*}
|
|
|
|
We can then subtract and add the previous equations to obtain the following equations:
|
|
#+begin_important
|
|
#+NAME: eq:du_coupled
|
|
\begin{equation}
|
|
m \ddot{d_u} + (k - m\dot{\theta}^2) d_u = F_u + 2 m\dot{d_v}\dot{\theta} + m d_v\ddot{\theta}
|
|
\end{equation}
|
|
#+NAME: eq:dv_coupled
|
|
\begin{equation}
|
|
m \ddot{d_v} + (k - m\dot{\theta}^2) d_v = F_v - 2 m\dot{d_u}\dot{\theta} - m d_u\ddot{\theta}
|
|
\end{equation}
|
|
#+end_important
|
|
|
|
We obtain two differential equations that are coupled through:
|
|
- *Euler forces*: $m d_v \ddot{\theta}$
|
|
- *Coriolis forces*: $2 m \dot{d_v} \dot{\theta}$
|
|
|
|
Without the coupling terms, each equation is the equation of a one degree of freedom mass-spring system with mass $m$ and stiffness $k- m\dot{\theta}^2$.
|
|
Thus, the term $- m\dot{\theta}^2$ acts like a negative stiffness (due to *centrifugal forces*).
|
|
|
|
The forces induced by the rotating reference frame are independent of the stiffness of the actuator.
|
|
The resulting effect of those forces should then be higher when using softer actuators.
|
|
|
|
** Numerical Values for the NASS
|
|
#+begin_src matlab :exports none :results silent :noweb yes
|
|
<<matlab-init>>
|
|
#+end_src
|
|
|
|
Let's define the parameters for the NASS.
|
|
#+begin_src matlab :exports none :results silent
|
|
mlight = 35; % Mass for light sample [kg]
|
|
mheavy = 85; % Mass for heavy sample [kg]
|
|
|
|
wlight = 2*pi; % Max rot. speed for light sample [rad/s]
|
|
wheavy = 2*pi/60; % Max rot. speed for heavy sample [rad/s]
|
|
|
|
kvc = 1e3; % Voice Coil Stiffness [N/m]
|
|
kpz = 1e8; % Piezo Stiffness [N/m]
|
|
|
|
wdot = 1; % Maximum rotation acceleration [rad/s2]
|
|
|
|
d = 0.01; % Maximum excentricity from rotational axis [m]
|
|
ddot = 0.2; % Maximum Horizontal speed [m/s]
|
|
|
|
save('./mat/parameters.mat');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :results table :exports results
|
|
labels = {'Light sample mass [kg]', ...
|
|
'Heavy sample mass [kg]', ...
|
|
'Max rot. speed - light [rpm]', ...
|
|
'Max rot. speed - heavy [rpm]', ...
|
|
'Voice Coil Stiffness [N/m]', ...
|
|
'Piezo Stiffness [N/m]', ...
|
|
'Max rot. acceleration [rad/s2]', ...
|
|
'Max mass excentricity [m]', ...
|
|
'Max Horizontal speed [m/s]'};
|
|
data = [mlight, mheavy, 60*wlight/2/pi, 60*wheavy/2/pi, kvc, kpz, wdot, d, ddot];
|
|
data2orgtable(data', labels, {}, ' %.1e ')
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
| Light sample mass [kg] | 3.5e+01 |
|
|
| Heavy sample mass [kg] | 8.5e+01 |
|
|
| Max rot. speed - light [rpm] | 6.0e+01 |
|
|
| Max rot. speed - heavy [rpm] | 1.0e+00 |
|
|
| Voice Coil Stiffness [N/m] | 1.0e+03 |
|
|
| Piezo Stiffness [N/m] | 1.0e+08 |
|
|
| Max rot. acceleration [rad/s2] | 1.0e+00 |
|
|
| Max mass excentricity [m] | 1.0e-02 |
|
|
| Max Horizontal speed [m/s] | 2.0e-01 |
|
|
|
|
** Euler and Coriolis forces - Numerical Result
|
|
First we will determine the value for Euler and Coriolis forces during regular experiment.
|
|
- *Euler forces*: $m d_v \ddot{\theta}$
|
|
- *Coriolis forces*: $2 m \dot{d_v} \dot{\theta}$
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
Felight = mlight*d*wdot;
|
|
Feheavy = mheavy*d*wdot;
|
|
|
|
Fclight = 2*mlight*ddot*wlight;
|
|
Fcheavy = 2*mheavy*ddot*wheavy;
|
|
#+end_src
|
|
|
|
The obtained values are displayed in table ref:tab:euler_coriolis.
|
|
|
|
#+begin_src matlab :results value table :exports results :post addhdr(*this*)
|
|
data = [Fclight, Fcheavy ;
|
|
Felight, Feheavy];
|
|
data2orgtable(data, {'Coriolis', 'Euler'}, {'Light', 'Heavy'}, ' %.1fN ')
|
|
#+end_src
|
|
|
|
#+NAME: tab:euler_coriolis
|
|
#+CAPTION: Euler and Coriolis forces for the NASS
|
|
#+RESULTS:
|
|
| | Light | Heavy |
|
|
|----------+-------+-------|
|
|
| Coriolis | 88.0N | 3.6N |
|
|
| Euler | 0.4N | 0.8N |
|
|
|
|
** Negative Spring Effect - Numerical Result
|
|
The negative stiffness due to the rotation is equal to $-m{\omega_0}^2$.
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
Klight = mlight*wlight^2;
|
|
Kheavy = mheavy*wheavy^2;
|
|
#+end_src
|
|
|
|
The values for the negative spring effect are displayed in table ref:tab:negative_spring.
|
|
|
|
This is definitely negligible when using piezoelectric actuators. It may not be the case when using voice coil actuators.
|
|
|
|
#+begin_src matlab :results value table :exports results :post addhdr(*this*)
|
|
data = [Klight, Kheavy];
|
|
data2orgtable(data, {'Neg. Spring'}, {'Light', 'Heavy'}, ' %.1f[N/m] ')
|
|
#+end_src
|
|
|
|
#+NAME: tab:negative_spring
|
|
#+CAPTION: Negative Spring effect
|
|
#+RESULTS:
|
|
| | Light | Heavy |
|
|
|-------------+-------------+----------|
|
|
| Neg. Spring | 1381.7[N/m] | 0.9[N/m] |
|
|
|
|
** Limitations due to coupling
|
|
To simplify, we consider a constant rotating speed $\dot{\theta} = \omega_0$ and thus $\ddot{\theta} = 0$.
|
|
|
|
From equations eqref:eq:du_coupled and eqref:eq:dv_coupled, we obtain:
|
|
\begin{align*}
|
|
(m s^2 + (k - m{\omega_0}^2)) d_u &= F_u + 2 m {\omega_0} s d_v \\
|
|
(m s^2 + (k - m{\omega_0}^2)) d_v &= F_v - 2 m {\omega_0} s d_u \\
|
|
\end{align*}
|
|
|
|
From second equation:
|
|
\[ d_v = \frac{1}{m s^2 + (k - m{\omega_0}^2)} F_v - \frac{2 m {\omega_0} s}{m s^2 + (k - m{\omega_0}^2)} d_u \]
|
|
|
|
And we re-inject $d_v$ into the first equation:
|
|
\begin{equation*}
|
|
(m s^2 + (k - m{\omega_0}^2)) d_u = F_u + \frac{2 m {\omega_0} s}{m s^2 + (k - m{\omega_0}^2)} F_v - \frac{(2 m {\omega_0} s)^2}{m s^2 + (k - m{\omega_0}^2)} d_u
|
|
\end{equation*}
|
|
|
|
\begin{equation*}
|
|
\frac{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2}{m s^2 + (k - m{\omega_0}^2)} d_u = F_u + \frac{2 m {\omega_0} s}{m s^2 + (k - m{\omega_0}^2)} F_v
|
|
\end{equation*}
|
|
|
|
Finally we obtain $d_u$ function of $F_u$ and $F_v$.
|
|
\[ d_u = \frac{m s^2 + (k - m{\omega_0}^2)}{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2} F_u + \frac{2 m {\omega_0} s}{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2} F_v \]
|
|
|
|
Similarly we can obtain $d_v$ function of $F_u$ and $F_v$:
|
|
\[ d_v = \frac{m s^2 + (k - m{\omega_0}^2)}{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2} F_v - \frac{2 m {\omega_0} s}{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2} F_u \]
|
|
|
|
The two previous equations can be written in a matrix form:
|
|
#+begin_important
|
|
#+NAME: eq:coupledplant
|
|
\begin{equation}
|
|
\begin{bmatrix} d_u \\ d_v \end{bmatrix} =
|
|
\frac{1}{(m s^2 + (k - m{\omega_0}^2))^2 + (2 m {\omega_0} s)^2}
|
|
\begin{bmatrix}
|
|
ms^2 + (k-m{\omega_0}^2) & 2 m \omega_0 s \\
|
|
-2 m \omega_0 s & ms^2 + (k-m{\omega_0}^2) \\
|
|
\end{bmatrix}
|
|
\begin{bmatrix} F_u \\ F_v \end{bmatrix}
|
|
\end{equation}
|
|
#+end_important
|
|
|
|
Then, coupling is negligible if $|-m \omega^2 + (k - m{\omega_0}^2)| \gg |2 m {\omega_0} \omega|$.
|
|
|
|
*** Numerical Analysis
|
|
We plot on the same graph $\frac{|-m \omega^2 + (k - m {\omega_0}^2)|}{|2 m \omega_0 \omega|}$ for the voice coil and the piezo:
|
|
- with the light sample (figure ref:fig:coupling_light).
|
|
- with the heavy sample (figure ref:fig:coupling_heavy).
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
f = logspace(-1, 3, 1000);
|
|
|
|
figure;
|
|
hold on;
|
|
plot(f, abs(-mlight*(2*pi*f).^2 + kvc - mlight * wlight^2)./abs(2*mlight*wlight*2*pi*f), 'DisplayName', 'Voice Coil')
|
|
plot(f, abs(-mlight*(2*pi*f).^2 + kpz - mlight * wlight^2)./abs(2*mlight*wlight*2*pi*f), 'DisplayName', 'Piezo')
|
|
plot(f, ones(1, length(f)), 'k--', 'HandleVisibility', 'off')
|
|
hold off;
|
|
xlim([f(1), f(end)]);
|
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
|
xlabel('Frequency [Hz]');
|
|
legend('Location', 'northeast');
|
|
#+end_src
|
|
|
|
#+NAME: fig:coupling_light
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/coupling_light.pdf" :var figsize="normal-normal" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:coupling_light
|
|
#+CAPTION: Relative Coupling for light mass and high rotation speed
|
|
#+RESULTS: fig:coupling_light
|
|
[[file:figs/coupling_light.png]]
|
|
#+begin_src matlab :exports none :results silent
|
|
figure;
|
|
hold on;
|
|
plot(f, abs(-mheavy*(2*pi*f).^2 + kvc - mheavy * wheavy^2)./abs(2*mheavy*wheavy*2*pi*f), 'DisplayName', 'Voice Coil')
|
|
plot(f, abs(-mheavy*(2*pi*f).^2 + kpz - mheavy * wheavy^2)./abs(2*mheavy*wheavy*2*pi*f), 'DisplayName', 'Piezo')
|
|
plot(f, ones(1, length(f)), 'k--', 'HandleVisibility', 'off')
|
|
hold off;
|
|
xlim([f(1), f(end)]);
|
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
|
xlabel('Frequency [Hz]');
|
|
legend('Location', 'northeast');
|
|
#+end_src
|
|
|
|
#+NAME: fig:coupling_heavy
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/coupling_heavy.pdf" :var figsize="normal-normal" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:coupling_heavy
|
|
#+CAPTION: Relative Coupling for heavy mass and low rotation speed
|
|
#+RESULTS: fig:coupling_heavy
|
|
[[file:figs/coupling_heavy.png]]
|
|
|
|
#+begin_important
|
|
Coupling is higher for actuators with small stiffness.
|
|
#+end_important
|
|
|
|
** Limitations due to negative stiffness effect
|
|
If $\max{\dot{\theta}} \ll \sqrt{\frac{k}{m}}$, then the negative spring effect is negligible and $k - m\dot{\theta}^2 \approx k$.
|
|
|
|
Let's estimate what is the maximum rotation speed for which the negative stiffness effect is still negligible ($\omega_\text{max} = 0.1 \sqrt{\frac{k}{m}}$). Results are shown table ref:tab:negative_stiffness.
|
|
#+begin_src matlab :results table :exports results :post addhdr(*this*)
|
|
data = 0.1*60*(1/2/pi)*[sqrt(kvc/mlight), sqrt(kpz/mlight);
|
|
sqrt(kvc/mheavy), sqrt(kpz/mheavy)];
|
|
data2orgtable(data, {'Light', 'Heavy'}, {'Voice Coil', 'Piezo'}, ' %.0f[rpm] ')
|
|
#+end_src
|
|
|
|
#+NAME: tab:negative_stiffness
|
|
#+CAPTION: Maximum rotation speed at which negative stiffness is negligible ($0.1\sqrt{\frac{k}{m}}$)
|
|
#+RESULTS:
|
|
| | Voice Coil | Piezo |
|
|
|-------+------------+-----------|
|
|
| Light | 5[rpm] | 1614[rpm] |
|
|
| Heavy | 3[rpm] | 1036[rpm] |
|
|
|
|
The negative spring effect is proportional to the rotational speed $\omega$.
|
|
The system dynamics will be much more affected when using soft actuator.
|
|
|
|
#+begin_important
|
|
Negative stiffness effect has very important effect when using soft actuators.
|
|
#+end_important
|
|
|
|
The system can even goes unstable when $m \omega^2 > k$, that is when the centrifugal forces are higher than the forces due to stiffness.
|
|
|
|
From this analysis, we can determine the lowest practical stiffness that is possible to use: $k_\text{min} = 10 m \omega^2$ (table sec:tab:min_k)
|
|
|
|
#+begin_src matlab :results table :exports results :post addhdr(*this*)
|
|
data = 10*[mlight*2*pi, mheavy*2*pi/60]
|
|
data2orgtable(data, {'k min [N/m]'}, {'Light', 'Heavy'}, ' %.0f ')
|
|
#+end_src
|
|
|
|
#+NAME: tab:min_k
|
|
#+CAPTION: Minimum possible stiffness
|
|
#+RESULTS:
|
|
| | Light | Heavy |
|
|
|-------------+-------+-------|
|
|
| k min [N/m] | 2199 | 89 |
|
|
|
|
** Effect of rotation speed on the plant
|
|
As shown in equation eqref:eq:coupledplant, the plant changes with the rotation speed $\omega_0$.
|
|
|
|
Then, we compute the bode plot of the direct term and coupling term for multiple rotating speed.
|
|
|
|
Then we compare the result between voice coil and piezoelectric actuators.
|
|
|
|
*** Voice coil actuator
|
|
#+begin_src matlab :exports none :results silent
|
|
m = mlight;
|
|
k = kvc;
|
|
|
|
ws = linspace(0, 2*pi, 5); % Rotation speed vector [rad/s]
|
|
|
|
Gs = {zeros(1, length(ws))};
|
|
Gcs = {zeros(1, length(ws))};
|
|
|
|
for i = 1:length(ws)
|
|
w = ws(i);
|
|
Gs(i) = {(m*s^2 + (k-m*w^2))/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)};
|
|
Gcs(i) = {(2*m*w*s)/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)};
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
freqs = logspace(-2, 1, 1000);
|
|
|
|
figure;
|
|
ax1 = subaxis(2,1,1);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, abs(squeeze(freqresp(Gs{i}, freqs, 'Hz'))));
|
|
end
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
set(gca, 'XTickLabel',[]);
|
|
ylabel('Magnitude [m/N]');
|
|
|
|
ax2 = subaxis(2,1,2);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
|
|
end
|
|
hold off;
|
|
yticks(-180:90:180);
|
|
ylim([-180 180]);
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
legend('Location', 'southwest');
|
|
linkaxes([ax1,ax2],'x');
|
|
#+end_src
|
|
|
|
#+NAME: fig:G_ws_vc
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/G_ws_vc.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:G_ws_vc
|
|
#+CAPTION: Bode plot of the direct transfer function term (from $F_u$ to $D_u$) for multiple rotation speed - Voice coil
|
|
#+RESULTS: fig:G_ws_vc
|
|
[[file:figs/G_ws_vc.png]]
|
|
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
freqs = logspace(-2, 1, 1000);
|
|
|
|
figure;
|
|
ax1 = subaxis(2,1,1);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, abs(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))));
|
|
end
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
set(gca, 'XTickLabel',[]);
|
|
ylabel('Magnitude [m/N]');
|
|
|
|
ax2 = subaxis(2,1,2);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
|
|
end
|
|
hold off;
|
|
yticks(-180:90:180);
|
|
ylim([-180 180]);
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
legend('Location', 'southwest');
|
|
linkaxes([ax1,ax2],'x');
|
|
#+end_src
|
|
|
|
#+NAME: fig:Gc_ws_vc
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/Gc_ws_vc.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:Gc_ws_vc
|
|
#+CAPTION: caption
|
|
#+RESULTS: fig:Gc_ws_vc
|
|
[[file:figs/Gc_ws_vc.png]]
|
|
|
|
*** Piezoelectric actuator
|
|
#+begin_src matlab :exports none :results silent
|
|
m = mlight;
|
|
k = kpz;
|
|
|
|
ws = linspace(0, 2*pi, 5); % Rotation speed vector [rad/s]
|
|
|
|
Gs = {zeros(1, length(ws))};
|
|
Gcs = {zeros(1, length(ws))};
|
|
|
|
for i = 1:length(ws)
|
|
w = ws(i);
|
|
Gs(i) = {(m*s^2 + (k-m*w^2))/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)};
|
|
Gcs(i) = {(2*m*w*s)/((m*s^2 + (k - m*w^2))^2 + (2*m*w*s)^2)};
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
freqs = logspace(2, 3, 1000);
|
|
|
|
figure;
|
|
ax1 = subaxis(2,1,1);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, abs(squeeze(freqresp(Gs{i}, freqs, 'Hz'))));
|
|
end
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
set(gca, 'XTickLabel',[]);
|
|
ylabel('Magnitude [m/N]');
|
|
|
|
ax2 = subaxis(2,1,2);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
|
|
end
|
|
hold off;
|
|
yticks(-180:90:180);
|
|
ylim([-180 180]);
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
legend('Location', 'southwest');
|
|
linkaxes([ax1,ax2],'x');
|
|
#+end_src
|
|
|
|
#+NAME: fig:G_ws_pz
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/G_ws_pz.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:G_ws_pz
|
|
#+CAPTION: Bode plot of the direct transfer function term (from $F_u$ to $D_u$) for multiple rotation speed - Piezoelectric actuator
|
|
#+RESULTS: fig:G_ws_pz
|
|
[[file:figs/G_ws_pz.png]]
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
figure;
|
|
ax1 = subaxis(2,1,1);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, abs(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))));
|
|
end
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
set(gca, 'XTickLabel',[]);
|
|
ylabel('Magnitude [m/N]');
|
|
|
|
ax2 = subaxis(2,1,2);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(Gcs{i}, freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
|
|
end
|
|
hold off;
|
|
yticks(-180:90:180);
|
|
ylim([-180 180]);
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
legend('Location', 'southwest');
|
|
linkaxes([ax1,ax2],'x');
|
|
#+end_src
|
|
|
|
#+NAME: fig:Gc_ws_pz
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/Gc_ws_pz.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:Gc_ws_pz
|
|
#+CAPTION: Bode plot of the coupling transfer function term (from $F_u$ to $D_v$) for multiple rotation speed - Piezoelectric actuator
|
|
#+RESULTS: fig:Gc_ws_pz
|
|
[[file:figs/Gc_ws_pz.png]]
|
|
|
|
*** Analysis
|
|
When the rotation speed is null, the coupling terms are equal to zero and the diagonal terms corresponds to one degree of freedom mass spring system.
|
|
|
|
When the rotation speed in not null, the resonance frequency is duplicated into two pairs of complex conjugate poles.
|
|
|
|
As the rotation speed increases, one of the two resonant frequency goes to lower frequencies as the other one goes to higher frequencies.
|
|
|
|
The poles of the coupling terms are the same as the poles of the diagonal terms. The magnitude of the coupling terms are increasing with the rotation speed.
|
|
|
|
#+begin_important
|
|
As shown in the previous figures, the system with voice coil is much more sensitive to rotation speed.
|
|
#+end_important
|
|
|
|
*** Campbell diagram
|
|
The poles of the system are computed for multiple values of the rotation frequency. To simplify the computation of the poles, we add some damping to the system.
|
|
|
|
#+begin_src matlab :results silent :exports code
|
|
m = mlight;
|
|
k = kvc;
|
|
c = 0.1*sqrt(k*m);
|
|
|
|
wsvc = linspace(0, 10, 100); % [rad/s]
|
|
|
|
polesvc = zeros(2, length(wsvc));
|
|
|
|
for i = 1:length(wsvc)
|
|
polei = pole(1/((m*s^2 + c*s + (k - m*wsvc(i)^2))^2 + (2*m*wsvc(i)*s)^2));
|
|
polesvc(:, i) = sort(polei(imag(polei) > 0));
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :results silent :exports code
|
|
m = mlight;
|
|
k = kpz;
|
|
c = 0.1*sqrt(k*m);
|
|
|
|
wspz = linspace(0, 1000, 100); % [rad/s]
|
|
|
|
polespz = zeros(2, length(wspz));
|
|
|
|
for i = 1:length(wspz)
|
|
polei = pole(1/((m*s^2 + c*s + (k - m*wspz(i)^2))^2 + (2*m*wspz(i)*s)^2));
|
|
polespz(:, i) = sort(polei(imag(polei) > 0));
|
|
end
|
|
#+end_src
|
|
|
|
We then plot the real and imaginary part of the poles as a function of the rotation frequency (figures ref:fig:poles_w_vc and ref:fig:poles_w_pz).
|
|
|
|
When the real part of one pole becomes positive, the system goes unstable.
|
|
|
|
For the voice coil (figure ref:fig:poles_w_vc), the system is unstable when the rotation speed is above 5 rad/s. The real and imaginary part of the poles of the system with piezoelectric actuators are changing much less (figure ref:fig:poles_w_pz).
|
|
|
|
#+begin_src matlab :results silent :exports none
|
|
figure;
|
|
% Amplitude
|
|
ax1 = subplot(1,2,1);
|
|
hold on;
|
|
plot(wsvc, real(polesvc(1, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wsvc, real(polesvc(2, :)))
|
|
plot(wsvc, zeros(size(wsvc)), 'k--')
|
|
hold off;
|
|
xlabel('Rotation Frequency [rad/s]');
|
|
ylabel('Pole Real Part');
|
|
ax2 = subplot(1,2,2);
|
|
hold on;
|
|
plot(wsvc, imag(polesvc(1, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wsvc, -imag(polesvc(1, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wsvc, imag(polesvc(2, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wsvc, -imag(polesvc(2, :)))
|
|
hold off;
|
|
xlabel('Rotation Frequency [rad/s]');
|
|
ylabel('Pole Imaginary Part');
|
|
#+end_src
|
|
|
|
#+NAME: fig:poles_w_vc
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes :mkdirp yes
|
|
#+begin_src matlab :var filepath="figs/poles_w_vc.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:poles_w_vc
|
|
#+CAPTION: Real and Imaginary part of the poles of the system as a function of the rotation speed - Voice Coil and light sample
|
|
#+RESULTS: fig:poles_w_vc
|
|
[[file:figs/poles_w_vc.png]]
|
|
|
|
|
|
#+begin_src matlab :results silent :exports none
|
|
figure;
|
|
% Amplitude
|
|
ax1 = subplot(1,2,1);
|
|
hold on;
|
|
plot(wspz, real(polespz(1, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wspz, real(polespz(2, :)))
|
|
plot(wspz, zeros(size(wspz)), 'k--')
|
|
hold off;
|
|
xlabel('Rotation Frequency [rad/s]');
|
|
ylabel('Pole Real Part');
|
|
ax2 = subplot(1,2,2);
|
|
hold on;
|
|
plot(wspz, imag(polespz(1, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wspz, -imag(polespz(1, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wspz, imag(polespz(2, :)))
|
|
set(gca,'ColorOrderIndex',1)
|
|
plot(wspz, -imag(polespz(2, :)))
|
|
hold off;
|
|
xlabel('Rotation Frequency [rad/s]');
|
|
ylabel('Pole Imaginary Part');
|
|
#+end_src
|
|
|
|
#+NAME: fig:poles_w_pz
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/poles_w_pz.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:poles_w_pz
|
|
#+CAPTION: Real and Imaginary part of the poles of the system as a function of the rotation speed - Piezoelectric actuator and light sample
|
|
#+RESULTS: fig:poles_w_pz
|
|
[[file:figs/poles_w_pz.png]]
|
|
|
|
* Control Strategies
|
|
<<sec:control_strategies>>
|
|
** Measurement in the fixed reference frame
|
|
First, let's consider a measurement in the fixed referenced frame.
|
|
|
|
The transfer function from actuator $[F_u, F_v]$ to sensor $[D_x, D_y]$ is then $G(\theta)$.
|
|
|
|
Then the measurement is subtracted to the reference signal $[r_x, r_y]$ to obtain the position error in the fixed reference frame $[\epsilon_x, \epsilon_y]$.
|
|
|
|
The position error $[\epsilon_x, \epsilon_y]$ is then express in the rotating frame corresponding to the actuators $[\epsilon_u, \epsilon_v]$.
|
|
|
|
Finally, the control low $K$ links the position errors $[\epsilon_u, \epsilon_v]$ to the actuator forces $[F_u, F_v]$.
|
|
|
|
The block diagram is shown on figure ref:fig:control_measure_fixed_2dof.
|
|
|
|
#+name: fig:control_measure_fixed_2dof
|
|
#+caption: Control with a measure from fixed frame
|
|
[[./figs/control_measure_fixed_2dof.png]]
|
|
|
|
The loop gain is then $L = G(\theta) K J(\theta)$.
|
|
|
|
One question we wish to answer is: is $G(\theta) J(\theta) = G(\theta_0) J(\theta_0)$?
|
|
|
|
** Measurement in the rotating frame
|
|
Let's consider that the measurement is made in the rotating reference frame.
|
|
|
|
The corresponding block diagram is shown figure ref:fig:control_measure_rotating_2dof
|
|
|
|
#+name: fig:control_measure_rotating_2dof
|
|
#+caption: Control with a measure from rotating frame
|
|
[[./figs/control_measure_rotating_2dof.png]]
|
|
|
|
The loop gain is $L = G K$.
|
|
|
|
* Multi Body Model - Simscape
|
|
:PROPERTIES:
|
|
:HEADER-ARGS:matlab+: :tangle matlab/simscape_analysis.m
|
|
:END:
|
|
<<sec:simscape>>
|
|
|
|
** Initialization
|
|
#+begin_src matlab :exports none :results silent :noweb yes
|
|
<<matlab-init>>
|
|
|
|
load('./mat/parameters.mat');
|
|
|
|
bode_opts = bodeoptions;
|
|
bode_opts.FreqUnits = 'Hz';
|
|
bode_opts.MagUnits = 'abs';
|
|
bode_opts.MagScale = 'log';
|
|
bode_opts.Grid = 'on';
|
|
bode_opts.PhaseVisible = 'off';
|
|
bode_opts.Title.FontSize = 10;
|
|
bode_opts.XLabel.FontSize = 10;
|
|
bode_opts.YLabel.FontSize = 10;
|
|
bode_opts.TickLabel.FontSize = 10;
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
open rotating_frame.slx
|
|
#+end_src
|
|
|
|
First we define the parameters that must be defined in order to run the Simscape simulation.
|
|
#+begin_src matlab :exports code :results silent
|
|
w = 2*pi; % Rotation speed [rad/s]
|
|
|
|
theta_e = 0; % Static measurement error on the angle theta [rad]
|
|
|
|
m = 5; % mass of the sample [kg]
|
|
|
|
mTuv = 30;% Mass of the moving part of the Tuv stage [kg]
|
|
kTuv = 1e8; % Stiffness of the Tuv stage [N/m]
|
|
cTuv = 0; % Damping of the Tuv stage [N/(m/s)]
|
|
#+end_src
|
|
|
|
Then, we defined parameters that will be used in the following analysis.
|
|
#+begin_src matlab :exports code :results silent
|
|
mlight = 5; % Mass for light sample [kg]
|
|
mheavy = 55; % Mass for heavy sample [kg]
|
|
|
|
wlight = 2*pi; % Max rot. speed for light sample [rad/s]
|
|
wheavy = 2*pi/60; % Max rot. speed for heavy sample [rad/s]
|
|
|
|
kvc = 1e3; % Voice Coil Stiffness [N/m]
|
|
kpz = 1e8; % Piezo Stiffness [N/m]
|
|
|
|
d = 0.01; % Maximum excentricity from rotational axis [m]
|
|
|
|
freqs = logspace(-2, 3, 1000); % Frequency vector for analysis [Hz]
|
|
#+end_src
|
|
|
|
** Identification in the rotating referenced frame
|
|
We initialize the inputs and outputs of the system to identify:
|
|
- Inputs: $f_u$ and $f_v$
|
|
- Outputs: $d_u$ and $d_v$
|
|
|
|
#+begin_src matlab :exports code :results silent
|
|
%% Options for Linearized
|
|
options = linearizeOptions;
|
|
options.SampleTime = 0;
|
|
|
|
%% Name of the Simulink File
|
|
mdl = 'rotating_frame';
|
|
|
|
%% Input/Output definition
|
|
io(1) = linio([mdl, '/fu'], 1, 'input');
|
|
io(2) = linio([mdl, '/fv'], 1, 'input');
|
|
|
|
io(3) = linio([mdl, '/du'], 1, 'output');
|
|
io(4) = linio([mdl, '/dv'], 1, 'output');
|
|
#+end_src
|
|
|
|
We start we identify the transfer functions at high speed with the light sample.
|
|
#+begin_src matlab :exports code :results silent
|
|
w = wlight; % Rotation speed [rad/s]
|
|
m = mlight; % mass of the sample [kg]
|
|
|
|
kTuv = kpz;
|
|
Gpz_light = linearize(mdl, io, 0.1);
|
|
Gpz_light.InputName = {'Fu', 'Fv'};
|
|
Gpz_light.OutputName = {'Du', 'Dv'};
|
|
|
|
kTuv = kvc;
|
|
Gvc_light = linearize(mdl, io, 0.1);
|
|
Gvc_light.InputName = {'Fu', 'Fv'};
|
|
Gvc_light.OutputName = {'Du', 'Dv'};
|
|
#+end_src
|
|
|
|
Then we identify the system with an heavy mass and low speed.
|
|
#+begin_src matlab :exports code :results silent
|
|
w = wheavy; % Rotation speed [rad/s]
|
|
m = mheavy; % mass of the sample [kg]
|
|
|
|
kTuv = kpz;
|
|
Gpz_heavy = linearize(mdl, io, 0.1);
|
|
Gpz_heavy.InputName = {'Fu', 'Fv'};
|
|
Gpz_heavy.OutputName = {'Du', 'Dv'};
|
|
|
|
kTuv = kvc;
|
|
Gvc_heavy = linearize(mdl, io, 0.1);
|
|
Gvc_heavy.InputName = {'Fu', 'Fv'};
|
|
Gvc_heavy.OutputName = {'Du', 'Dv'};
|
|
#+end_src
|
|
|
|
** Coupling ratio between $f_{uv}$ and $d_{uv}$
|
|
In order to validate the equations written, we can compute the coupling ratio using the simscape model and compare with the equations.
|
|
|
|
From the previous identification, we plot the coupling ratio in both case (figure ref:fig:coupling_ratio_light_heavy).
|
|
|
|
We obtain the same result than the analytical case (figures ref:fig:coupling_light and ref:fig:coupling_heavy).
|
|
#+begin_src matlab :results silent :exports none
|
|
figure;
|
|
hold on;
|
|
plot(freqs, abs(squeeze(freqresp(Gvc_light('Du', 'Fu'), freqs, 'Hz')))./abs(squeeze(freqresp(Gvc_light('Dv', 'Fu'), freqs, 'Hz'))));
|
|
plot(freqs, abs(squeeze(freqresp(Gpz_light('Du', 'Fu'), freqs, 'Hz')))./abs(squeeze(freqresp(Gpz_light('Dv', 'Fu'), freqs, 'Hz'))));
|
|
set(gca,'ColorOrderIndex',1);
|
|
plot(freqs, abs(squeeze(freqresp(Gvc_heavy('Du', 'Fu'), freqs, 'Hz')))./abs(squeeze(freqresp(Gvc_heavy('Dv', 'Fu'), freqs, 'Hz'))), '--');
|
|
plot(freqs, abs(squeeze(freqresp(Gpz_heavy('Du', 'Fu'), freqs, 'Hz')))./abs(squeeze(freqresp(Gpz_heavy('Dv', 'Fu'), freqs, 'Hz'))), '--');
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Coupling ratio');
|
|
legend({'light - VC', 'light - PZ', 'heavy - VC', 'heavy - PZ'}, 'Location', 'northeast')
|
|
#+end_src
|
|
|
|
#+NAME: fig:coupling_ratio_light_heavy
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes :mkdirp yes
|
|
#+begin_src matlab :var filepath="figs/coupling_ratio_light_heavy.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:coupling_ratio_light_heavy
|
|
#+CAPTION: Coupling ratio obtained with the Simscape model
|
|
#+RESULTS: fig:coupling_ratio_light_heavy
|
|
[[file:figs/coupling_ratio_light_heavy.png]]
|
|
|
|
** Plant Control - SISO approach
|
|
*** Plant identification
|
|
The goal is to study the control problems due to the coupling that appears because of the rotation.
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
%% Options for Linearized
|
|
options = linearizeOptions;
|
|
options.SampleTime = 0;
|
|
|
|
%% Name of the Simulink File
|
|
mdl = 'rotating_frame';
|
|
|
|
%% Input/Output definition
|
|
io(1) = linio([mdl, '/fu'], 1, 'input');
|
|
io(2) = linio([mdl, '/fv'], 1, 'input');
|
|
|
|
io(3) = linio([mdl, '/du'], 1, 'output');
|
|
io(4) = linio([mdl, '/dv'], 1, 'output');
|
|
#+end_src
|
|
|
|
First, we identify the system when the rotation speed is null and then when the rotation speed is equal to 60rpm.
|
|
|
|
The actuators are voice coil with some damping added.
|
|
|
|
The bode plot of the system not rotating and rotating at 60rpm is shown figure ref:fig:Gvc_speed.
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
w = 0; % Rotation speed [rad/s]
|
|
m = mlight; % mass of the sample [kg]
|
|
kTuv = kvc;
|
|
cTuv = 0.1*sqrt(kTuv*m);
|
|
|
|
Gvc = linearize(mdl, io, 0.1);
|
|
Gvc.InputName = {'Fu', 'Fv'};
|
|
Gvc.OutputName = {'Du', 'Dv'};
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
w = wlight; % Rotation speed [rad/s]
|
|
m = mlight; % mass of the sample [kg]
|
|
kTuv = kvc;
|
|
cTuv = 0.1*sqrt(kTuv*m);
|
|
|
|
Gtvc = linearize(mdl, io, 0.1);
|
|
Gtvc.InputName = {'Fu', 'Fv'};
|
|
Gtvc.OutputName = {'Du', 'Dv'};
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
figure;
|
|
bode(Gvc, Gtvc, bode_opts)
|
|
legend({'Gvc - $\omega = 0$', 'Gvc - $\omega = 60$rpm'}, 'Location', 'southwest');
|
|
title('');
|
|
#+end_src
|
|
|
|
#+NAME: fig:Gvc_speed
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/Gvc_speed.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:Gvc_speed
|
|
#+CAPTION: Change of transfer functions due to rotating speed
|
|
#+RESULTS: fig:Gvc_speed
|
|
[[file:figs/Gvc_speed.png]]
|
|
|
|
*** Effect of rotation speed
|
|
We first identify the system (voice coil and light mass) for multiple rotation speed.
|
|
Then we compute the bode plot of the diagonal element (figure ref:fig:Guu_ws) and of the coupling element (figure ref:fig:Guv_ws).
|
|
|
|
As the rotation frequency increases:
|
|
- one pole goes to lower frequencies while the other goes to higher frequencies
|
|
- one zero appears between the two poles
|
|
- the zero disappears when $\omega > \sqrt{\frac{k}{m}}$ and the low frequency pole becomes unstable (positive real part)
|
|
|
|
To stabilize the unstable pole, we need a control bandwidth of at least twice of frequency of the unstable pole.
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
ws = linspace(0, 2*pi, 5); % Rotation speed vector [rad/s]
|
|
m = mlight; % mass of the sample [kg]
|
|
|
|
kTuv = kvc;
|
|
cTuv = 0.1*sqrt(kTuv*m);
|
|
|
|
Gs = {zeros(1, length(ws))};
|
|
|
|
for i = 1:length(ws)
|
|
w = ws(i);
|
|
Gs{i} = linearize(mdl, io, 0.1);
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
freqs = logspace(-2, 2, 1000);
|
|
|
|
figure;
|
|
ax1 = subaxis(2,1,1);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, abs(squeeze(freqresp(Gs{i}(1, 1), freqs, 'Hz'))));
|
|
end
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
set(gca, 'XTickLabel',[]);
|
|
ylabel('Magnitude [m/N]');
|
|
|
|
ax2 = subaxis(2,1,2);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}(1, 1), freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
|
|
end
|
|
hold off;
|
|
yticks(-180:90:180);
|
|
ylim([-180 180]);
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
legend('Location', 'northeast');
|
|
linkaxes([ax1,ax2],'x');
|
|
#+end_src
|
|
|
|
#+NAME: fig:Guu_ws
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/Guu_ws.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:Guu_ws
|
|
#+CAPTION: Diagonal term as a function of the rotation frequency
|
|
#+RESULTS: fig:Guu_ws
|
|
[[file:figs/Guu_ws.png]]
|
|
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
freqs = logspace(-2, 2, 1000);
|
|
|
|
figure;
|
|
ax1 = subaxis(2,1,1);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, abs(squeeze(freqresp(Gs{i}(1, 2), freqs, 'Hz'))));
|
|
end
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
set(gca, 'XTickLabel',[]);
|
|
ylabel('Magnitude [m/N]');
|
|
|
|
ax2 = subaxis(2,1,2);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}(1, 2), freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
|
|
end
|
|
hold off;
|
|
yticks(-180:90:180);
|
|
ylim([-180 180]);
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
legend('Location', 'northeast');
|
|
linkaxes([ax1,ax2],'x');
|
|
#+end_src
|
|
|
|
#+NAME: fig:Guv_ws
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/Guv_ws.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:Guv_ws
|
|
#+CAPTION: Couplin term as a function of the rotation frequency
|
|
#+RESULTS: fig:Guv_ws
|
|
[[file:figs/Guv_ws.png]]
|
|
|
|
Then, we can look at the same plots for the piezoelectric actuator (figure ref:fig:Guu_ws_pz). The effect of the rotation frequency has very little effect on the dynamics of the system to control.
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
ws = linspace(0, 2*pi, 5); % Rotation speed vector [rad/s]
|
|
m = mlight; % mass of the sample [kg]
|
|
|
|
kTuv = kpz;
|
|
cTuv = 0.1*sqrt(kTuv*m);
|
|
|
|
Gs = {zeros(1, length(ws))};
|
|
|
|
for i = 1:length(ws)
|
|
w = ws(i);
|
|
Gs{i} = linearize(mdl, io, 0.1);
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
freqs = logspace(2, 3, 1000);
|
|
|
|
figure;
|
|
ax1 = subaxis(2,1,1);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, abs(squeeze(freqresp(Gs{i}(1, 1), freqs, 'Hz'))));
|
|
end
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
set(gca, 'XTickLabel',[]);
|
|
ylabel('Magnitude [m/N]');
|
|
|
|
ax2 = subaxis(2,1,2);
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(Gs{i}(1, 1), freqs, 'Hz'))), 'DisplayName', sprintf('w = %.0f [rpm]', ws(i)*60/2/pi));
|
|
end
|
|
hold off;
|
|
yticks(-180:90:180);
|
|
ylim([-180 180]);
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
legend('Location', 'northeast');
|
|
linkaxes([ax1,ax2],'x');
|
|
#+end_src
|
|
|
|
#+NAME: fig:Guu_ws_pz
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/Guu_ws_pz.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:Guu_ws_pz
|
|
#+CAPTION: Diagonal term as a function of the rotation frequency
|
|
#+RESULTS: fig:Guu_ws_pz
|
|
[[file:figs/Guu_ws_pz.png]]
|
|
|
|
*** Controller design
|
|
We design a controller based on the identification when the system is not rotating.
|
|
|
|
The obtained controller is a lead-lag controller with the following transfer function.
|
|
#+begin_src matlab :results none :exports code
|
|
Kll = 2.0698e09*(s+40.45)*(s+1.181)/((s+0.01)*(s+198.4)*(s+2790));
|
|
K = [Kll 0;
|
|
0 Kll];
|
|
#+end_src
|
|
|
|
The loop gain is displayed figure ref:fig:Gvc_loop_gain.
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
freqs = logspace(-2, 2, 1000);
|
|
|
|
figure;
|
|
% Amplitude
|
|
ax1 = subaxis(2,1,1);
|
|
hold on;
|
|
plot(freqs, abs(squeeze(freqresp(Gvc('Du', 'fu')*Kll, freqs, 'Hz'))), '-');
|
|
set(gca,'xscale','log'); set(gca,'yscale','log');
|
|
ylabel('Amplitude [m/N]');
|
|
set(gca, 'XTickLabel',[]);
|
|
hold off;
|
|
% Phase
|
|
ax2 = subaxis(2,1,2);
|
|
hold on;
|
|
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gvc('Du', 'fu')*Kll, freqs, 'Hz')))), '-');
|
|
set(gca,'xscale','log');
|
|
yticks(-180:180:180);
|
|
ylim([-180 180]);
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
hold off;
|
|
linkaxes([ax1,ax2],'x');
|
|
#+end_src
|
|
|
|
#+NAME: fig:Gvc_loop_gain
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/Gvc_loop_gain.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:Gvc_loop_gain
|
|
#+CAPTION: Loop gain obtained for a lead-lag controller on the system with a voice coil
|
|
#+RESULTS: fig:Gvc_loop_gain
|
|
[[file:figs/Gvc_loop_gain.png]]
|
|
|
|
*** Controlling the rotating system
|
|
We here want to see if the system is robust with respect to the rotation speed. We then use the controller based on the non-rotating system, and see if the system is stable and its dynamics.
|
|
|
|
We can then plot the same loop gain with the rotating system using the same controller (figure ref:fig:Gtvc_loop_gain). The result obtained is unstable.
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
freqs = logspace(-2, 2, 1000);
|
|
|
|
figure;
|
|
% Amplitude
|
|
ax1 = subaxis(2,1,1);
|
|
hold on;
|
|
plot(freqs, abs(squeeze(freqresp(Gtvc('Du', 'fu')*Kll, freqs, 'Hz'))), '-');
|
|
set(gca,'xscale','log'); set(gca,'yscale','log');
|
|
ylabel('Amplitude [m/N]');
|
|
set(gca, 'XTickLabel',[]);
|
|
hold off;
|
|
% Phase
|
|
ax2 = subaxis(2,1,2);
|
|
hold on;
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(Gtvc('Du', 'fu')*Kll, freqs, 'Hz'))), '-');
|
|
set(gca,'xscale','log');
|
|
yticks(-180:180:180);
|
|
ylim([-180 180]);
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
hold off;
|
|
linkaxes([ax1,ax2],'x');
|
|
#+end_src
|
|
|
|
#+NAME: fig:Gtvc_loop_gain
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/Gtvc_loop_gain.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:Gtvc_loop_gain
|
|
#+CAPTION: Loop gain with the rotating system
|
|
#+RESULTS: fig:Gtvc_loop_gain
|
|
[[file:figs/Gtvc_loop_gain.png]]
|
|
|
|
We can look at the poles of the system where we control only one direction ($u$ for instance). We obtain a pole with a positive real part.
|
|
|
|
#+begin_important
|
|
The system is then unstable when controlling only one direction.
|
|
#+end_important
|
|
|
|
#+begin_src matlab :results table :exports both
|
|
pole(feedback(Gtvc, blkdiag(Kll, 0)))
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
| -2798 |
|
|
| -58.916+94.248i |
|
|
| -58.916-94.248i |
|
|
| -71.644 |
|
|
| 3.1647 |
|
|
| -3.3034 |
|
|
| -1.1901 |
|
|
|
|
However, when we look at the poles of the closed loop with a diagonal controller, all the poles have negative real part and the system is stable.
|
|
#+begin_src matlab :results table :exports both
|
|
pole(feedback(Gtvc, blkdiag(Kll, Kll)))
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
| -2798+0.035765i |
|
|
| -2798-0.035765i |
|
|
| -56.414+105.34i |
|
|
| -56.414-105.34i |
|
|
| -64.495+79.314i |
|
|
| -64.495-79.314i |
|
|
| -68.509+13.499i |
|
|
| -68.509-13.499i |
|
|
| -1.1837+0.0041422i |
|
|
| -1.1837-0.0041422i |
|
|
|
|
Check stability of MIMO system.
|
|
#+begin_src matlab :result silent
|
|
isstable(1/(1+K*Gtvc))
|
|
isstable(Gtvc/(1+K*Gtvc))
|
|
isstable(Gtvc/(1+K*Gtvc))
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
: 0
|
|
|
|
*** Close loop performance
|
|
First, we create the closed loop systems. Then, we plot the transfer function from the reference signals $[\epsilon_u, \epsilon_v]$ to the output $[d_u, d_v]$ (figure ref:fig:perfcomp).
|
|
|
|
#+begin_src matlab :results none :exports code
|
|
S = eye(2)/(eye(2) + Gvc*K);
|
|
T = Gvc*K /(eye(2) + Gvc*K);
|
|
|
|
St = eye(2)/(eye(2) + Gtvc*K);
|
|
Tt = Gtvc*K/(eye(2) + Gtvc*K);
|
|
|
|
freqs = logspace(-3, 3, 1000);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :results none :exports code
|
|
figure;
|
|
bode(S, St, 2*pi*freqs, bode_opts)
|
|
#+end_src
|
|
|
|
#+begin_src matlab :results none :exports code
|
|
figure;
|
|
bode(T, Tt, 2*pi*freqs, bode_opts)
|
|
#+end_src
|
|
|
|
#+begin_src matlab :results none :exports code
|
|
freqs = logspace(-2, 2, 1000);
|
|
|
|
figure;
|
|
ax1 = subplot(1,2,1);
|
|
hold on;
|
|
plot(freqs, abs(squeeze(freqresp(Tvc(1, 1), freqs, 'Hz'))));
|
|
plot(freqs, abs(squeeze(freqresp(Ttvc(1, 1), freqs, 'Hz'))));
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Magnitude [m/N]');
|
|
legend({'w = 0 [rpm]', 'w = 60 [rpm]'}, 'Location', 'southwest')
|
|
title('$G_{r_u \to d_u}$')
|
|
|
|
ax2 = subplot(1,2,2);
|
|
hold on;
|
|
plot(freqs, abs(squeeze(freqresp(Tvc(1, 2), freqs, 'Hz'))));
|
|
plot(freqs, abs(squeeze(freqresp(Ttvc(1, 2), freqs, 'Hz'))));
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
ylim([1e-5, 1]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]');
|
|
title('$G_{r_u \to d_v}$')
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
#+end_src
|
|
|
|
#+NAME: fig:perfconp
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/perfconp.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:perfconp
|
|
#+CAPTION: Close loop performance for $\omega = 0$ and $\omega = 60 rpm$
|
|
#+RESULTS: fig:perfconp
|
|
[[file:figs/perfconp.png]]
|
|
|
|
** Plant Control - MIMO approach
|
|
*** TODO Analysis - SVD
|
|
\[ G = U \Sigma V^H \]
|
|
|
|
With:
|
|
- $\Sigma$ is an $2 \times 2$ matrix with 2 non-negative *singular values* $\sigma_i$, arranged in descending order along its main diagonal
|
|
- $U$ is an $2 \times 2$ unitary matrix. The columns vectors of $U$, denoted $u_i$, represent the *output directions* of the plant. They are orthonomal
|
|
- $V$ is an $2 \times 2$ unitary matrix. The columns vectors of $V$, denoted $v_i$, represent the *input directions* of the plant. They are orthonomal
|
|
|
|
We first look at the evolution of the singular values as a function of frequency (figure ref:fig:G_sigma).
|
|
#+begin_src matlab :exports none :results silent
|
|
ws = linspace(0, 2*pi, 5); % Rotation speed vector [rad/s]
|
|
m = mlight; % mass of the sample [kg]
|
|
|
|
kTuv = kvc;
|
|
cTuv = 0.1*sqrt(kTuv*m);
|
|
|
|
Gs = {zeros(1, length(ws))};
|
|
|
|
for i = 1:length(ws)
|
|
w = ws(i);
|
|
Gs{i} = linearize(mdl, io, 0.1);
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :results silent :exports none
|
|
freqs = logspace(-2, 2, 1000);
|
|
|
|
figure;
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
plot(freqs, abs(squeeze(freqresp(Gs{i}(1,1), freqs, 'Hz'))), '-');
|
|
end
|
|
hold off;
|
|
set(gca,'xscale','log'); set(gca,'yscale','log');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :results silent :exports none
|
|
freqs = logspace(-2, 1, 1000);
|
|
|
|
figure;
|
|
hold on;
|
|
for i = 1:length(ws)
|
|
sv = sigma(Gs{i}, 2*pi*freqs);
|
|
set(gca,'ColorOrderIndex',i)
|
|
plot(freqs, sv(1, :), 'DisplayName', sprintf('w = %.0f rpm', ws(i)*60/2/pi));
|
|
set(gca,'ColorOrderIndex',i)
|
|
plot(freqs, sv(2, :), '--', 'HandleVisibility', 'off');
|
|
end
|
|
hold off;
|
|
set(gca,'xscale','log'); set(gca,'yscale','log');
|
|
legend('location', 'southwest');
|
|
#+end_src
|
|
|
|
#+NAME: fig:G_sigma
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/G_sigma.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:G_sigma
|
|
#+CAPTION: caption
|
|
#+RESULTS: fig:G_sigma
|
|
#+CAPTION: Evolution of the singular values with frequency
|
|
[[file:figs/G_sigma.png]]
|
|
|
|
We compute
|
|
#+begin_src matlab :results silent :exports code
|
|
[U,S,V] = svd(freqresp(Gtvc, 2*pi*10));
|
|
#+end_src
|
|
|
|
#+begin_src matlab :results output :exports results
|
|
U, S, V
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
#+begin_example
|
|
U, S, V
|
|
U =
|
|
-0.707101109012986 - 0.00283224868340902i -0.707104254409621 - 0.00189034277692295i
|
|
0.00283224868340845 - 0.707101109012987i -0.00189034277692242 + 0.70710425440962i
|
|
S =
|
|
9.01532756059351e-06 0
|
|
0 6.01714794171208e-06
|
|
V =
|
|
0.707106781186547 + 0i 0.707106781186548 + 0i
|
|
-1.57009245868378e-16 + 0.707106781186548i 1.57009245868377e-16 - 0.707106781186547i
|
|
#+end_example
|
|
|
|
The input and output directions are related through the singular values
|
|
\[ G v_i = \sigma_i u_i \]
|
|
|
|
So, if we consider an input in the direction $v_i$, then the output is in the direction $u_i$. Furthermore, since $\normtwo{v_i}=1$ and $\normtwo{u_i}=1$, we see that *the singular value $\sigma_i$ directly gives the gain of the matrix $G$ in this direction*.
|
|
|
|
#+begin_src matlab
|
|
freqresp(Gtvc, 2*pi*10)*V(:, 1)
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
| -6.3747e-06-2.5534e-08i |
|
|
| 2.5534e-08-6.3747e-06i |
|
|
|
|
#+begin_src matlab
|
|
S(1)*U(:, 1)
|
|
#+end_src
|
|
|
|
#+RESULTS:
|
|
| -6.3747e-06-2.5534e-08i |
|
|
| 2.5534e-08-6.3747e-06i |
|
|
|
|
*** Closed loop SVD
|
|
#+begin_src matlab :results silent :exports none
|
|
figure;
|
|
sigma(Tvc, Ttvc)
|
|
#+end_src
|
|
|
|
** test
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
figure;
|
|
% Amplitude
|
|
ax1 = subaxis(2,1,1);
|
|
hold on;
|
|
plot(freqs, abs(squeeze(freqresp(Gvc_light('du', 'fu'), freqs, 'Hz'))), '-');
|
|
set(gca,'xscale','log'); set(gca,'yscale','log');
|
|
ylabel('Amplitude [m/N]');
|
|
set(gca, 'XTickLabel',[]);
|
|
hold off;
|
|
% Phase
|
|
ax2 = subaxis(2,1,2);
|
|
hold on;
|
|
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gvc_light('du', 'fu'), freqs, 'Hz')))), '-');
|
|
set(gca,'xscale','log');
|
|
yticks(-180:180:180);
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
hold off;
|
|
linkaxes([ax1,ax2],'x');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
figure;
|
|
% Amplitude
|
|
ax1 = subaxis(2,1,1);
|
|
hold on;
|
|
plot(freqs, abs(squeeze(freqresp(-Gvc_light('dv', 'fv'), freqs, 'Hz'))), '-');
|
|
set(gca,'xscale','log'); set(gca,'yscale','log');
|
|
ylabel('Amplitude [m/N]');
|
|
set(gca, 'XTickLabel',[]);
|
|
hold off;
|
|
% Phase
|
|
ax2 = subaxis(2,1,2);
|
|
hold on;
|
|
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(-Gvc_light('dv', 'fv'), freqs, 'Hz')))), '-');
|
|
set(gca,'xscale','log');
|
|
yticks(-180:180:180);
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
hold off;
|
|
linkaxes([ax1,ax2],'x');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
figure;
|
|
bode(Gpz_light, Gvc_light);
|
|
#+end_src
|
|
|
|
#+NAME: fig:coupling_simscape
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/coupling_simscape.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:coupling_simscape
|
|
#+CAPTION: caption
|
|
#+RESULTS: fig:coupling_simscape
|
|
[[file:figs/coupling_simscape_light.png]]
|
|
|
|
And then with the heavy sample.
|
|
#+begin_src matlab :exports code :results silent
|
|
rot_speed = wheavy;
|
|
angle_e = 0;
|
|
m = mheavy;
|
|
|
|
k = kpz;
|
|
c = 1e3;
|
|
Gpz_heavy = linearize(mdl, io, 0.1);
|
|
|
|
k = kvc;
|
|
c = 1e3;
|
|
Gvc_heavy = linearize(mdl, io, 0.1);
|
|
|
|
Gpz_heavy.InputName = {'Fu', 'Fv'};
|
|
Gpz_heavy.OutputName = {'Du', 'Dv'};
|
|
Gvc_heavy.InputName = {'Fu', 'Fv'};
|
|
Gvc_heavy.OutputName = {'Du', 'Dv'};
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none :results silent
|
|
figure;
|
|
bode(Gpz_heavy, Gvc_heavy);
|
|
#+end_src
|
|
|
|
#+NAME: fig:coupling_simscape_heavy
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/coupling_simscape_heavy.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:coupling_simscape_heavy
|
|
#+CAPTION: caption
|
|
#+RESULTS: fig:coupling_simscape_heavy
|
|
[[file:figs/coupling_simscape_heavy.png]]
|
|
|
|
Plot the ratio between the main transfer function and the coupling term:
|
|
#+begin_src matlab :results silent :exports none
|
|
freqs = logspace(-2, 3, 1000);
|
|
|
|
figure;
|
|
hold on;
|
|
plot(freqs, abs(squeeze(freqresp(Gvc_light('Du', 'Fu'), freqs, 'Hz')))./abs(squeeze(freqresp(Gvc_light('Dv', 'Fu'), freqs, 'Hz'))));
|
|
plot(freqs, abs(squeeze(freqresp(Gpz_light('Du', 'Fu'), freqs, 'Hz')))./abs(squeeze(freqresp(Gpz_light('Dv', 'Fu'), freqs, 'Hz'))));
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Coupling ratio');
|
|
legend({'Voice Coil', 'Piezoelectric'})
|
|
#+end_src
|
|
|
|
#+NAME: fig:coupling_ratio_simscape_light
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/coupling_ratio_simscape_light.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:coupling_ratio_simscape_light
|
|
#+CAPTION: caption
|
|
#+RESULTS: fig:coupling_ratio_simscape_light
|
|
[[file:figs/coupling_ratio_simscape_light.png]]
|
|
|
|
#+begin_src matlab :results silent :exports none
|
|
freqs = logspace(-2, 3, 1000);
|
|
|
|
figure;
|
|
hold on;
|
|
plot(freqs, abs(squeeze(freqresp(Gvc_heavy('Du', 'Fu'), freqs, 'Hz')))./abs(squeeze(freqresp(Gvc_heavy('Dv', 'Fu'), freqs, 'Hz'))));
|
|
plot(freqs, abs(squeeze(freqresp(Gpz_heavy('Du', 'Fu'), freqs, 'Hz')))./abs(squeeze(freqresp(Gpz_heavy('Dv', 'Fu'), freqs, 'Hz'))));
|
|
hold off;
|
|
xlim([freqs(1), freqs(end)]);
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Coupling ratio');
|
|
legend({'Voice Coil', 'Piezoelectric'})
|
|
#+end_src
|
|
|
|
#+NAME: fig:coupling_ratio_simscape_heavy
|
|
#+HEADER: :tangle no :exports results :results raw :noweb yes
|
|
#+begin_src matlab :var filepath="figs/coupling_ratio_simscape_heavy.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+LABEL: fig:coupling_ratio_simscape_heavy
|
|
#+CAPTION: caption
|
|
#+RESULTS: fig:coupling_ratio_simscape_heavy
|
|
[[file:figs/coupling_ratio_simscape_heavy.png]]
|
|
|
|
*** Low rotation speed and High rotation speed
|
|
#+begin_src matlab :exports code :results silent
|
|
rot_speed = 2*pi/60; angle_e = 0;
|
|
G_low = linearize(mdl, io, 0.1);
|
|
|
|
rot_speed = 2*pi; angle_e = 0;
|
|
G_high = linearize(mdl, io, 0.1);
|
|
|
|
G_low.InputName = {'Fu', 'Fv'};
|
|
G_low.OutputName = {'Du', 'Dv'};
|
|
G_high.InputName = {'Fu', 'Fv'};
|
|
G_high.OutputName = {'Du', 'Dv'};
|
|
#+end_src
|
|
|
|
#+begin_src matlab :results silent
|
|
figure;
|
|
bode(G_low, G_high);
|
|
#+end_src
|
|
|
|
** Identification in the fixed frame
|
|
Let's define some options as well as the inputs and outputs for linearization.
|
|
#+begin_src matlab :exports code :results silent
|
|
%% Options for Linearized
|
|
options = linearizeOptions;
|
|
options.SampleTime = 0;
|
|
|
|
%% Name of the Simulink File
|
|
mdl = 'rotating_frame';
|
|
|
|
%% Input/Output definition
|
|
io(1) = linio([mdl, '/fx'], 1, 'input');
|
|
io(2) = linio([mdl, '/fy'], 1, 'input');
|
|
|
|
io(3) = linio([mdl, '/dx'], 1, 'output');
|
|
io(4) = linio([mdl, '/dy'], 1, 'output');
|
|
#+end_src
|
|
|
|
We then define the error estimation of the error and the rotational speed.
|
|
#+begin_src matlab :exports code :results silent
|
|
%% Run the linearization
|
|
angle_e = 0;
|
|
rot_speed = 0;
|
|
#+end_src
|
|
|
|
Finally, we run the linearization.
|
|
#+begin_src matlab :exports code :results silent
|
|
G = linearize(mdl, io, 0);
|
|
|
|
%% Input/Output names
|
|
G.InputName = {'Fx', 'Fy'};
|
|
G.OutputName = {'Dx', 'Dy'};
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports code :results silent
|
|
%% Run the linearization
|
|
angle_e = 0;
|
|
rot_speed = 2*pi;
|
|
Gr = linearize(mdl, io, 0);
|
|
|
|
%% Input/Output names
|
|
Gr.InputName = {'Fx', 'Fy'};
|
|
Gr.OutputName = {'Dx', 'Dy'};
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports code :results silent
|
|
%% Run the linearization
|
|
angle_e = 1*2*pi/180;
|
|
rot_speed = 2*pi;
|
|
Ge = linearize(mdl, io, 0);
|
|
|
|
%% Input/Output names
|
|
Ge.InputName = {'Fx', 'Fy'};
|
|
Ge.OutputName = {'Dx', 'Dy'};
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports code :results silent
|
|
figure;
|
|
bode(G);
|
|
% exportFig('G_x_y', 'wide-tall');
|
|
|
|
figure;
|
|
bode(Ge);
|
|
% exportFig('G_x_y_e', 'normal-normal');
|
|
#+end_src
|
|
|
|
** Identification from actuator forces to displacement in the fixed frame
|
|
#+begin_src matlab :exports code :results silent
|
|
%% Options for Linearized
|
|
options = linearizeOptions;
|
|
options.SampleTime = 0;
|
|
|
|
%% Name of the Simulink File
|
|
mdl = 'rotating_frame';
|
|
|
|
%% Input/Output definition
|
|
io(1) = linio([mdl, '/fu'], 1, 'input');
|
|
io(2) = linio([mdl, '/fv'], 1, 'input');
|
|
|
|
io(3) = linio([mdl, '/dx'], 1, 'output');
|
|
io(4) = linio([mdl, '/dy'], 1, 'output');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports code :results silent
|
|
rot_speed = 2*pi;
|
|
angle_e = 0;
|
|
G = linearize(mdl, io, 0.0);
|
|
|
|
G.InputName = {'Fu', 'Fv'};
|
|
G.OutputName = {'Dx', 'Dy'};
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports code :results silent
|
|
rot_speed = 2*pi;
|
|
angle_e = 0;
|
|
G1 = linearize(mdl, io, 0.4);
|
|
|
|
G1.InputName = {'Fu', 'Fv'};
|
|
G1.OutputName = {'Dx', 'Dy'};
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports code :results silent
|
|
rot_speed = 2*pi;
|
|
angle_e = 0;
|
|
G2 = linearize(mdl, io, 0.8);
|
|
|
|
G2.InputName = {'Fu', 'Fv'};
|
|
G2.OutputName = {'Dx', 'Dy'};
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports code :results silent
|
|
figure;
|
|
bode(G, G1, G2);
|
|
exportFig('G_u_v_to_x_y', 'wide-tall');
|
|
#+end_src
|
|
|
|
** Effect of the X-Y stage stiffness
|
|
<<sec:effect_stiffness>>
|
|
|
|
*** TODO At full speed, check how the coupling changes with the stiffness of the actuators
|
|
|
|
* Control Implementation
|
|
<<sec:control>>
|
|
|
|
** Measurement in the fixed reference frame
|
|
|
|
* Bibliography :ignore:
|
|
# #+BIBLIOGRAPHY: /home/tdehaeze/MEGA/These/Ressources/references.bib plain option:-a option:-noabstract option:-nokeywords option:-noheader option:-nofooter option:-nobibsource limit:t
|