diff options
Diffstat (limited to '')
| -rw-r--r-- | midterm.tex | 447 |
1 files changed, 447 insertions, 0 deletions
diff --git a/midterm.tex b/midterm.tex new file mode 100644 index 0000000..78e60fa --- /dev/null +++ b/midterm.tex @@ -0,0 +1,447 @@ +\documentclass[a4paper, 10pt]{IEEEtran} + +\usepackage{array} +\usepackage{booktabs} +\usepackage{siunitx} +\usepackage{float} + +\usepackage{amsmath} +\usepackage{amssymb} +\usepackage{oststud} + +\usepackage{graphicx} +\usepackage{subcaption} +\usepackage{pgfplots} +\pgfplotsset{compat=1.18} +\usepackage{tikz} +\usetikzlibrary{3d,perspective,matrix,calc} + +% Bibliography IEEE style +\usepackage[ + backend=biber, + style=ieee, + citestyle=numeric-comp, + sortcites=true +]{biblatex} + +% Allow breaking long URLs in the bibliography +\setcounter{biburllcpenalty}{7000} +\setcounter{biburlucpenalty}{8000} + +% Add a bibliography resource +\addbibresource{RCCO.bib} + +\title{ + Robust $\mathcal{H}_\infty$ Controller Synthesis for Hovering Ducted-Fan VTOL + Micro-UAV +} +\author{ + \IEEEauthorblockN{Naoki Sean Pross} \IEEEauthorrefmark{1} + % \thanks{ + % } + \thanks{ + Midterm project for course \emph{Robust Control and Convex Optimization} + taught by Prof. Dr. R. S. Smith at ETH Zurich. + \IEEEauthorrefmark{1}Student ID:~\texttt{19-150-903} + Email:~\texttt{npross@student.ethz.ch}. + Thanks to T. Rothlin for providing the electromechanical parameters for the + UAV. + } +} + +\begin{document} + +\maketitle + +\begin{abstract} + This paper models an experimental vertical take-off landing (VTOL) micro + unmanned aerial vehicle (UAV) device developed at the University of Applied + Sciences OST RJ in Rapperswil. The model is then converted to a linear time + invariant (LTI) state space plant to design a robust controller using + $\mathcal{H}_\infty$ synthesis that tracks a piecewise constant spatial + position reference. In conclusion simulated preliminary results on the + controller performance are presented. +\end{abstract} + +\section{Introduction} + +The presented VTOL micro UAV is a device developed by T. Rothlin +\cite{rothlin_optimization_2024} at the University of Applied Sciences OST RJ, +Rapperswil that can be steered only by controlling the propeller velocity and +the angle of attack of four flaps under the ducted fan. Currently there is not +a control algorithm for the UAV and the uncertainties from aerodynamics and +gyroscopic effects caused by the propeller make this an interesting candidate +for robust controller synthesis. + +% \subsection{Notation} +% +% Throughout this text we will make use a notation customary to mechanics, using +% boldface letters for vectors and matrices and a hat to indicate unit vectors. +% Spatial coordinates are grouped into a position vector $\vec{P}$ to avoid +% confusing the $x$ and $y$ coordinates with the state space model's output and +% state vectors $\vec{y}$ and $\vec{x}$. + +% \begin{table} +% \centering +% \caption{ +% Summary of notation for physical quantities. +% \label{tab:sym-phy} +% } +% \begin{tabular}{cl} +% \toprule +% Symbol & Physical Quantity \\ +% \midrule +% $\vec{P}$ & Position vector \\ +% $\vec{F}$ & Total force \\ +% $\vec{M}$ & Total torque \\ +% \bottomrule +% \end{tabular} +% \end{table} + +% \subsection{Organization} +% +% Section \ref{sec:modelling} presents a derivation of a model for the plant, and +% then \S\ref{sec:control} illustrates the controller design. Finally, in +% \S\ref{sec:simulation} the closed loop design is simulated and briefly +% discussed. + +\section{System Modelling} +\label{sec:modelling} + +% \subsection{Reference Frames} + +To model the dynamics of the ducted-fan UAV two reference frames are required: +an inertial frame and a body-frame attached to the center of mass +\cite{stengel_flight_2022}. In the inertial frame we work in the base given by +the unit vectors $\{\uvec{\imath},\uvec{\jmath},\uvec{k}\}$, whereas in the +body-frame we use $\{\uvec{x}, \uvec{y}, \uvec{z}\}$ together with the Euler +angles $\vec{\Theta} = [ \phi ~ \vartheta ~ \psi ]^\mathsf{T}$. The angular +velocity $\vec{\Omega} = [p ~ q ~ r]^\mathsf{T}$ is then related to +$\vec{\Theta}$ by $\dot{\vec{\Theta}} = \mx{U} \vec{\Omega}$. To rotate from +the inertial frame to body frame we use the $SO(3)$ matrix $\mx{R}$. Both +$\mx{U}$ and $\mx{R}$ can be found in \cite{stengel_flight_2022}. + +% \begin{equation} +% \dot{\vec{\Theta}} = \begin{bmatrix} +% 1 & S_\phi T_\vartheta & C_\phi T_\vartheta \\ +% 0 & C_\phi & - S_\phi \\ +% 0 & S_\phi / C_\vartheta & C_\phi / C_\vartheta +% \end{bmatrix} +% \vec{\Omega} +% = \mx{u} \vec{\omega}, +% \end{equation} +% wherein we use $S_\alpha$, $C_\alpha$ and $T_\alpha$ as shorthand for $\sin +% \alpha$, $\cos \alpha$ and $\tan \alpha$ respectively. To rotate from the +% inertial frame to body frame we use the $SO(3)$ matrix +% \begin{equation*} +% \mx{R} = \begin{bmatrix} +% C_\vartheta C_\psi & C_\vartheta S_\psi & - S_\vartheta \\ +% S_\phi S_\vartheta C_\psi - C_\phi S_\psi +% & S_\phi S_\vartheta C_\psi + C_\phi C_\psi +% & S_\phi C_\vartheta \\ +% C_\phi S_\theta C_\psi + S_\phi S_\psi +% & C_\phi S_\theta S_\psi - S_\phi C_\psi +% & C_\phi C_\vartheta +% \end{bmatrix}. +% \end{equation*} + +\subsection{Equations of Motion} + +Consider a simplified rigid-body physical model sketched in +Fig.~\ref{fig:phy-sketch}. With respect to the inertial frame the equations of +motion from the Newton-Euler formalism are +\begin{subequations} + \begin{align} + m \ddot{\vec{P}} &= \mt{\mx{R}} \vec{F} \\ + \vec{J} \dot{\vec{\Omega}} &= + - \vec{\Omega} \crossp \vec{J} \vec{\Omega} + \vec{\tau}, + \end{align} +\end{subequations} +where $m$ is the total mass, $\vec{P}$ is the position, $\vec{F}$ and +$\vec{\tau}$ the total force and torque in the body frame respectively, and +finally $\mx{J}$ the moment of inertia, which is assumed to be a diagonal +matrix. + +\begin{figure} + \centering + % \includegraphics[width=.7\linewidth, trim={0 10 0 0}, clip]{fig/sketch_hand} + \input{fig/sketch} + \caption{ + Sketch of the rigid body physical model for the UAV. + % \texttt{TODO: replace with computer-drawn sketch.} + \label{fig:phy-sketch} + } +\end{figure} + +\begin{figure*}[t] + \centering + \begin{tikzpicture}[font = \footnotesize] + \begin{axis}[ + width = 42mm, height = 28mm, + scale only axis, + grid = major, + grid style = {densely dotted, gray!50}, + ylabel = {Flap Angle $\alpha$ / degrees}, + xlabel = {Time $t$ / seconds}, + legend entries = { + $\alpha_1(t)$, + $\alpha_2(t)$, + $\alpha_3(t)$, + $\alpha_4(t)$ + }, + legend columns = 2, + title = {\scshape Step Response of Flap Angles}, + /tikz/line join=bevel, + ] + \addplot+ [mark=none] table [x index = 0, y index = 7] {fig/stepsim.dat}; + \addplot+ [mark=none] table [x index = 0, y index = 8] {fig/stepsim.dat}; + \addplot+ [mark=none] table [x index = 0, y index = 9] {fig/stepsim.dat}; + \addplot+ [mark=none] table [x index = 0, y index = 10] {fig/stepsim.dat}; + \end{axis} + \end{tikzpicture} + \hfill + \begin{tikzpicture}[font = \footnotesize, ] + \begin{axis}[ + width = 42mm, height = 28mm, + scale only axis, + grid = major, + grid style = {densely dotted, gray!50}, + ylabel = {Euler Angle / degrees}, + xlabel = {Time $t$ / seconds}, + legend entries = { + $\phi(t)$, + $\vartheta(t)$, + $\psi(t)$, + }, + legend columns = 2, + legend style = {anchor = south east, at = {(.95,.05)}}, + title = {\scshape Step Response of Euler Angles}, + ] + \addplot+ [thick, mark=none] table [x index = 0, y index = 11] {fig/stepsim.dat}; + \addplot+ [thick, mark=none] table [x index = 0, y index = 12] {fig/stepsim.dat}; + \addplot+ [thick, mark=none] table [x index = 0, y index = 13] {fig/stepsim.dat}; + \end{axis} + \end{tikzpicture} + \hfill + \begin{tikzpicture}[font = \footnotesize] + \pgfplotsset{set layers} + \begin{axis}[ + width = 42mm, height = 28mm, + scale only axis, + grid = major, + grid style = {densely dotted, gray!50}, + axis y line* = left, + ylabel = {Coordinate / meters}, + xlabel = {Time $t$ / seconds}, + legend entries = { + $x(t)$, + $y(t)$ + }, + title = {\scshape Step Response of Position}, + ] + \addplot+ [thick, mark=none] table [x index = 0, y index = 1] {fig/stepsim.dat}; + \addplot+ [thick, mark=none] table [x index = 0, y index = 2] {fig/stepsim.dat}; + \end{axis} + \begin{axis}[ + width = 42mm, height = 28mm, + scale only axis, + % grid = major, + % grid style = {densely dotted, gray!50}, + axis y line* = right, + axis x line = none, + ylabel = {Velocity / \unit{\meter\per\second}}, + ] + \addplot+ [thick, mark=none, densely dotted, opacity=.8] + table [x index = 0, y index = 4] {fig/stepsim.dat}; + \addplot+ [thick, mark=none, densely dotted, opacity=.8] + table [x index = 0, y index = 5] {fig/stepsim.dat}; + \end{axis} + \end{tikzpicture} + \caption{ + Simulated step responses along $x$ of the closed loop $\mathcal{H}_\infty$ + controller design. In the third plot from the left, continuous lines show + the position error while dotted lines with the right axis are the velocity + in the corresponding directions. + \label{fig:step} + } +\end{figure*} + +In the body frame we model the total force $\vec{F}$ acting on the UAV +following \cite{naldi_design_2010} by considering a thrust force from the +ducted-fan $\vec{F}_T = - k_T \omega^2 \uvec{z}$ as a function of the +propeller's angular velocity $\omega$. Because of the geometry we approximate +the generated air velocity field in the duct as being constant and collinear +to $\uvec{z}$ \cite{pereira_hover_2008} +\begin{equation} + \vec{\nu} = \sqrt{\frac{F_T}{2a\varrho \pi^2}} \uvec{z} + = \frac{\omega}{\pi} \sqrt{\frac{k_T}{2a\varrho}} \uvec{z}. +\end{equation} +The \emph{drag} and \emph{lift} forces generated by the interaction of each +flap with the air velocity field are then $\vec{F}_d = \frac{1}{2} \varrho S +C_d\nu^2 \uvec{z}$ and $\vec{F}_\ell = \frac{1}{2} \varrho S C_\ell \nu^2 +\uvec{n}$ respectively ($\uvec{n}$ being given by the orientation of the flap). +The density of air $\varrho$ is assumed to be constant, while the drag +coefficients $S$, $C_d$ and $C_\ell$ depend on the angle of attack $\alpha$ of +the flap making them quite difficult to determine. Hereinafter, under the +assumption that $\alpha$ is small they will be approximated with $C_d = c_d +\alpha^2 + c_{0}$, $C_\ell = c_\ell \alpha$ \cite{stengel_flight_2022, +pereira_hover_2008} and $S$ is considered a constant. Finally, gravity adds a +term $\vec{F}_g = mg \mt{\mx{R}} \uvec{k}$. + +For the total torque it is derived from the geometry that each flap induces +$\vec{\tau}_{f} = (\frac{1}{3}a \uvec{t} + d\uvec{z}) \crossp (\vec{F}_d + +\vec{F}_\ell)$, with $\uvec{t} = \uvec{n} \crossp \uvec{z}$. In addition there +is also a torque induced by the gyroscopic procession $\vec{\tau}_g = +\mx{R}^\mathsf{T} \omega J_r \uvec{k} \crossp \vec{\Omega}$ with $J_r$ being +the inertia of the propeller with respect to its spin axis +\cite{naldi_design_2010}. The resulting total quantities are then +\begin{subequations} + \begin{align} + \vec{F} &= mg \mt{\mx{R}} \uvec{k} + -k_T\omega^2 \uvec{z} \nonumber \\ + &\quad + \frac{\varrho S \nu^2}{2} \sum_i + \left(c_d \alpha^2_i + c_0 \right) \uvec{z} + + c_\ell \alpha_i \uvec{n}_i, \\ + \vec{\tau} &= \omega J_r \mx{R}^\mathsf{T} (\uvec{k} \crossp \vec{\Omega}) + + \frac{\varrho c_\ell S \nu^2 d}{2} + \uvec{z} \crossp \sum_i \alpha_i \uvec{n}_i. + \end{align} +\end{subequations} + +% \subsection{Wind Gusts Model} +% +% \texttt{TODO: add disturbance from wind modelled as simple, slowly changing +% stochastic model following empirical approach as done in +% \cite{chirarattananon_dynamics_2017}. Use 3D state with azimuth $\gamma$, +% zenith $\beta$ and magnitude $W$ to generate force and torque. Maybe +% make cyclostationary?} + + +\subsection{Linearized Dynamics and Plant Model} + +To synthesize a controller we simplify the dynamics to a state-space LTI model +$\dot{\vec{x}} = \mx{A}_x \vec{x} + \mx{B}_{xu} \vec{u}$, $\vec{y} = +\mx{C}_{yx} \vec{x} + \mx{D}_{yu}\vec{u}$ with state and inputs given by +\begin{equation} + \vec{x} = \mt{\begin{bmatrix} + \mt{\vec{P}} & + \mt{\dot{\vec{P}}} & + \mt{\vec{\Theta}} & + \mt{\vec{\Omega}} + \end{bmatrix}}, + \quad + \vec{u} = \mt{\begin{bmatrix} + \mt{\vec{\alpha}} & \omega + \end{bmatrix}} +\end{equation} +respectively. The linearization is performed in a stationary hovering state at +height $h$ above the ground and yaw angle of \ang{45} or $\vec{P}_0 = +-h\uvec{k}$, $\vec{\Theta}_0 = \frac{\pi}{4} \uvec{\psi}$, thus $x_0 = +\mt{[\mt{\vec{P}}_0 ~ \mt{\vec{0}} ~ \mt{\vec{\Theta}}_0 ~ \mt{\vec{0}}]}$ and +$u_0 = \mt{[ \mt{\vec{0}} ~ \omega_0 ]}$, where $\omega_0 \approx \sqrt{mg / +k_T}$ is the propeller's angular velocity to make the UAV hover. + +% \subsection{Plant Model} + +The state is assumed to be known, as there is an inertial measurement unit with +a dedicated sensor fusion chip onboard \cite{rothlin_optimization_2024}. Thus, +we only need to model a measurement delay $T_m$ and the actuators. The +linearized dynamics are extended with 2\textsuperscript{nd} order Padé +approximant $G_{m}(s) \approx e^{-sT_m}$ for the output delay and two transfer +functions $G_\alpha(s)$, $G_\omega (s)$ for the flaps and thruster +respectively. Specifically, $G_\alpha(s)$ is a 2\textsuperscript{nd} order low +pass filter (LPF) with a slight overshoot, while $G_\omega(s)$ is a first order +LPF. All inputs and outputs are then normalized and a source of additive +Gaussian noise $\vec{w}_\alpha \sim \mathcal{N}(\vec{0}, \mx{I})$ is introduced +to model the wiggling in the servo motor of each flap. The resulting +plant has 4 unobservable and 4 uncontrollable modes, however it remains both +detectable and stabilizable. + +\section{Controller Design} +\label{sec:control} + +% \subsection{Control Problem} + +The control objective is to track a piecewise constant position reference +$\vec{r}(t) \in \mathbb{R}^3$ in the inertial frame. The controller should move +the UAV to $\vec{r}(t)$ with control actions that do not exceed the actuator +limits $\alpha \in (\underline{\alpha}, \overline{\alpha})$, $\omega \in [0, +\overline{\omega})$, do not bring the angular components of the state too far +from their set point, i.e. $\vec{\Theta} \in (-\varphi,\varphi)^3$, and do not +exceed the velocity limits $\|\dot{\vec{P}}\|_1 \leq \overline{v}$. + +% \subsection{$\mathcal{H}_\infty$ Synthesis} + +To perform $\mathcal{H}_\infty$ synthesis the system is first rewritten in the +form presented in \cite{skogestad_multivariable_2010}, in terms of control +inputs $\vec{u} \in \mathbb{R}^5$, exogenous inputs $\vec{w} = +\mt{[\mt{\vec{w}}_\alpha ~ \mt{\vec{r}}]} \in \mathbb{R}^7$, control outputs +$\vec{y}_m = \mt{[\mt{(\vec{r} - \vec{y})} ~ \mt{\dot{\vec{P}}} ~ +\mt{\vec{\Theta}} ~ \mt{\vec{\Omega}}]}$ and errors $\vec{e} = +\mt{[\mt{\vec{\alpha}}~\omega~\mt{(\vec{r} - \vec{y})} ~ \mt{\dot{\vec{P}}} ~ +\mt{\vec{\Theta}}]}$ and partitioned accordingly: +\begin{equation} + \begin{bmatrix} \vec{e} \\ \vec{y}_m \end{bmatrix} + = \begin{bmatrix} + \mx{A}_{ew} & \mx{B}_{eu} \\ + \mx{C}_{yw} & \mx{D}_{yu} + \end{bmatrix} + \begin{bmatrix} \vec{w} \\ \vec{u} \end{bmatrix} + = \mx{G}_\text{nom} + \begin{bmatrix} \vec{w} \\ \vec{u} \end{bmatrix}. +\end{equation} + +Then for each component of the error a performance weight function is defined +resulting in 5 transfer functions $W_{P,\alpha}(s)$, $W_{P,\omega}(s)$, +$W_{P,\vec{P}}(s)$, $W_{P,\dot{\vec{P}}}(s)$, $W_{P,\vec{\Theta}}(s)$. The +$\mathcal{H}_\infty$ synthesis is performed using the Riccati method yielding a +stable controller $\mx{K}_\infty$ with $\gamma = +\|\mathcal{L}(\mx{G}_\text{nom}, \mx{K}_\infty)\|_\infty \approx 0.9091$. Here +we used the notation $\mathcal{L}$ for the linear fractional transformation +interconnection. + +\section{Preliminary Results} \label{sec:simulation} + +The resulting closed loop plant $\mx{G}_\text{cl} = +\mathcal{L}(\mx{G}_\text{nom}, \mx{K}_\infty)$ is simulated with a step of +\qty{50}{\cm} along the $x$ spatial coordinate and shown in +Fig.~\ref{fig:step}. Currently, it takes the UAV about 30 seconds to track a +half meter step, while respecting the constraints on the actuators. Although +stable and with good noise rejection characteristics the current design does +not take into account large inaccuracies caused by the highly non linear +equations of motion. For higher speeds necessitate higher pitch and roll +angles, a more performant controller will need to cope with higher modelling +errors. Hence the next controller will be designed with $\mu$-synthesis which +can take into account the aforementioned modelling uncertainties. + +\printbibliography + +\iffalse +\clearpage +\appendix + +\subsection{Linearized Model Details} + +\begin{equation*} + A = \begin{bmatrix} + 0 & I & 0 & 0 \\ + 0 + & \frac{1}{m} \mt{\mx{R}} \frac{\partial \vec{F}}{\partial \dot{\vec{P}}} + & \frac{1}{m} \frac{\partial (\mt{\mx{R}} \vec{F})}{\partial \vec{\Theta}} + & \frac{1}{m} \mt{\mx{R}} \frac{\partial \vec{F}}{\partial \vec{\Omega}} \\ + 0 & 0 + & \frac{\partial \mx{U}}{\partial \vec{\Theta}} \vec{\Omega} + & \mx{U} \\ + 0 + & \minv{\mx{J}} \frac{\partial \vec{\tau}}{\partial \dot{\vec{P}}} + & \minv{\mx{J}} \frac{\partial \vec{\tau}}{\partial \vec{\Theta}} + & \minv{\mx{J}} \left[ + \frac{\partial \vec{\tau}}{\partial \vec{\Omega}} + - \frac{\partial (\vec{\Omega} \crossp \mx{J} \vec{\Omega})}% + {\partial \vec{\Omega}} \right] + \end{bmatrix} +\end{equation*} +\fi + +\end{document} +% vim:ts=2 sw=2 et spell: |