# A low computation saturated controller for a quadcopter position tracking with stability and convoluted input constraints guarantees

Huu-Thinh Do, Franco Blanchini, Ionela Prodan

► To cite this version:

Huu-Thinh Do, Franco Blanchini, Ionela Prodan. A low computation saturated controller for a quadcopter position tracking with stability and convoluted input constraints guarantees. European Journal of Control, 2025, pp.101236. 10.1016/j.ejcon.2025.101236. hal-05074409

HAL Id: hal-05074409

https://hal.science/hal-05074409v1

Submitted on 20 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A low computation saturated controller for a
quadcopter position tracking with stability and
convoluted input constraints guarantees

Huu-Thinh Do<sup>a,*</sup>, Franco Blanchini<sup>b</sup>, Ionela Prodan<sup>a</sup>

<sup>a</sup>Univ. Grenoble Alpes, Grenoble INP<sup>1</sup>, LCIS, Valence 26000, France

<sup>b</sup>Dipartimento di Matematica e Informatica, Università di Udine, Udine 33100, Italy

Abstract

In the class of differentially flat systems, nonlinear dynamics can be trans-
formed into an equivalent linear representation via a coordinate change and
an input mapping. Despite this beneficial linearization, one significant hin-
drance is that the system constraints have a convoluted form in the new
space (the flat output space), which is usually disregarded in the literature
or bypassed via conservative approximation of the feasible domain. In this
paper, based on a pre-stabilization design and a novel input saturation com-
ponent, the nonlinear input constraints are particularly handled for the po-
sition tracking problem of quadcopters in the flat output space. While the
system's stability is shown to hold by Lyapunov arguments, the practical
viability of the method is validated both in simulation and experiments over
a nano-drone platform. Moreover, with the explicit solution introduced for
the proposed saturation function, this flatness-based saturated controller not

Email addresses: huu-thinh.do@lcis.grenoble-inp.fr (Huu-Thinh Do),
franco.blanchini@uniud.it (Franco Blanchini),
ionela.prodan@lcis.grenoble-inp.fr (Ionela Prodan)

¹Institute of Engineering and Management Univ. Grenoble Alpes.
This work has benefited from a French government grant managed by the Agence Nationale
de la Recherche under the France 2030 program, reference ANR-23-IACL-0006 and Inter-
national Research Booster Univ. Grenoble Alpes 2024, project "InteractingCPS" grant
no. C900IDX24V30. Franco Blanchini has been supported by European Union, NextGen-
eration EU (Grant Uniud-DM737). Franco Blanchini has been partially supported by
the Italian Ministry for Research, Project PRIDE, PRIN 2022 by the European Union -
NextGenerationEU grant no. 2022LP77J4.

<footnote>*Corresponding author</footnote>
only ensures stability and constraints satisfaction, but also requires low com-
putational effort, allowing possible embedded implementations.

*Keywords:* saturated control, constraint-handling controller, ellipsoidal invariant set, differential flatness, quadcopters.

# 1. Introduction

*Motivation:* Of all the unmanned aerial vehicles, quadcopters and other classes of multirotors have drawn remarkable attention in both research and practical use thanks to their ability to accomplish vertical take-off and landing as well as static hovering, which are applied in various logistic, surveillance or law-enforcing applications [1, 2, 3]. With the natural decoupling of the dynamics, this class of rotorcrafts is usually governed in a hierarchical architecture. The cascade design includes the position and attitude control (outer loop and inner loop) to handle the translational and the rotational dynamics of the vehicle, respectively, with specific requirements on their different bandwidths or settling times. Thanks to this separation, where the inner loop is required to be stable in higher frequency than the outer loop [4, 5], the work load for control design is significantly relieved by considering the two fully-actuated layers separately. Several validated results provided for both loops can be found in [6, 7] and the references therein. With this architecture, different quadcopter autopilot brands and their specifications are provided in [8]. Although the deep investigation in both layers of control is undeniably important in research, it is commonly required when it comes to agile flights or aggressively high-speed piloting [9, 10]. Meanwhile, in civil and domestic applications (e.g., delivery, irrigation), with commercial drone platforms where high frequency built-in attitude controllers are provided, the navigation problem can be reduced to the customized construction of the position control [11, 8]. However, even when only the outer loop is considered, the guidance and control of these thrust-propelled vehicles remain problematic due to the model's nonlinearity and physical constraints.

*Literature review and open issues:* Typically, to deal with the complexity and input limits, the dynamics of the quadcopter are linearly approximated via the Taylor expansion near its equilibrium points and then governed with classic or modern control approaches (e.g. fuzzy logic, LQR or learning-based approaches) [12, 13, 14, 15]. Although this strategy is proven to be generic and applicable, one potential downside is its sensitivity to uncertainty, which
is brought on by the approximation-based inexact dynamics. Moreover, the closed-loop stability proof is often given only for the approximated system, in lieu of the investigated nonlinear model, hence neglecting the entire left-over remainder of the approximation. On the contrary, as an intrinsic property of the quadcopter, the system's dynamics is known to be *differentially flat* [16, 17, 18]. Namely, all the system's states and inputs can be algebraically expressed in terms of a special output, called *flat output*, and a finite number of its derivatives. Consequently, thanks to this input-output relationship, the system's motion can be exactly described by finite chains of trivial integrators after a diffeomorphism and a dynamic feedback linearization law. In other words, the system can be represented by a linear controllable system in new coordinates, called the *flat output space*. The control problem is now solved by closing the loop for such linear dynamics, converting back the control input to the original representation and applying it to the system. This approximation-free approach undeniably simplifies the system description, yet its corresponding constraints are typically convoluted (in general, becoming nonlinear) and challenging to manage. Hence, control designers end up having to choose whether to deal with nonlinearity in the original model or in the constraints with the new coordinate system. With a view to pursue the latter, several approaches have been employed. In [10, 9], the constraints are taken into account offline (prior to the flight). Specifically, a dynamically-constrained feedforward reference was given to the vehicle in addition to the standard feedback control to improve accuracy and to follow aggressive trajectories. Besides, other solutions include stabilizing the linear model by conservatively sketching the new constraints through some linear approximation and subsequently handling them online. More specifically, by proposing an inner polyhedron inscribed in the feasible domain, the system was governed under the Model Predictive Control (MPC) framework where a constrained optimization problem is solved in-flight to optimally regulate the system [19, 20, 21]. Furthermore, in [22, 23], although the state prediction is computed via the quadcopter's nonlinear model, its linearized representation was employed to ensure the recursive feasibility and the stability of the nonlinear MPC approach via the existence of a feasible local controller. An computationally improved design without terminal ingredients was later studied in [24]. Yet, as a price of solving a nonlinear optimization problem online, these methods have a high computational demand, thus leading to long sampling time. Inspired by the classical approach in [25], the idea of controlling the exactly linearized model was also developed in a nested con-
trol fashion for the quadcopter’s position control [26, 18], which guarantees
the constraint satisfaction by appropriately saturating the amplitude of the
new input. This design was validated to be effective and computationally ef-
ficient, since the online calculation requires only algebraic computation and
the standard saturation operators. However, the nonlinear constraints were,
again, sidestepped by a conservative inner approximation. Illustrations and
comparisons of the approximated sets employed can be found in [27].

*Contributions:* In this work, focusing on the position controller which stabilizes the quadcopter’s translational dynamics, we present an investigation of the system in its flat output space together with the distorted constraint description. Then, based on this new characterization, our goal is to obtain a computationally low-cost control design (lower than 5 to 50 milliseconds as commonly achieved with MPC-based implementations [19, 23, 28, 24]), guaranteeing both constraint satisfaction and stability. To this aim, a Flatness-Based Saturated Control design (FBSC) is proposed to rigorously take into account the non-standard saturation constraints, which, to the best of the authors’ knowledge, is novel to the literature. More specifically, unlike the standard amplitude bound independently imposed on each control after a feedback linearization [26, 18, 22, 21], or the polytopic approximation in [27], we deal with a less conservative constraint set for this control layer at a price of its non-linearity and the convoluted coupling between the new inputs. Additionally, the explicit solution for such a particular saturation function is also introduced with a view to reducing the online computational cost. The efficiency of the solution is shown via the formulation’s simplicity and comparison with the existing MPC solution. Finally, from an experimental standpoint, the gain selection for our position controller, which is directly restrained by the inner loop (attitude control), will also be discussed. Briefly, the contributions are summarized as follows:

* design an optimization-based saturated controller which ensures both stability and constraint satisfaction, based on the linearized dynamics in the new coordinate and its associated nonlinear constraints;

* propose a computational low-cost explicit procedure to implement the controller;

* validate the benefits of the method through simulation and experimental tests by using the Crazyflie 2.1 nano-drone platform. The video for the experiment is available at https://youtu.be/0YcEM280gts.
**Organization:** The remainder of the paper is structured as the following. Section 2 presents the quadcopter model used for outer loop design, its flat characterization and operational constraints. Furthermore, invariance and Lyapunov-based tools are briefly introduced, which will prove instrumental for the proposed controller stability analysis. Section 3 presents the saturated control design together with theoretical proofs for constraint satisfaction and stability. Section 4 shows the validation of the proposed controller through simulations and experiments under different trajectory scenarios. Comparisons and discussions highlight the strengths of our novel flatness-based saturated controller design. Finally, Section 5 summarizes and provides insights on future work.

**Notations:** Denote bold upper-case letters as matrices with suitable dimension. Next, $I_n$ and $0_n$ denote an $n \times n$ identity and zero square matrix, respectively. Similarly, $0_{a \times b}$, represents a matrix of dimension $a \times b$ all of whose components are zero. diag($\cdot$) generates a diagonal matrix created by the employed components. Next, vectors are represented by bold letters (e.g, $\mathbf{v}, \mathbf{x}$). $\|\mathbf{v}\|_Q$ and $\|\mathbf{v}\|_2$ denote the weighted norm $\sqrt{\mathbf{v}^\top Q \mathbf{v}}$ and the Euclidean norm $\sqrt{\mathbf{v}^\top \mathbf{v}}$, respectively. Then, $\mathcal{B}(r), \mathcal{B}_M(\varepsilon) \subset \mathbb{R}^n$ defines the balls characterized by $\mathcal{B}(r) = \{\mathbf{x} \in \mathbb{R}^n : \|\mathbf{x}\|_2^2 \le r\}$ and $\mathcal{B}_M(\varepsilon) = \{\mathbf{x} \in \mathbb{R}^n : \|\mathbf{x}\|_M^2 \le \varepsilon\}$, respectively. Next, with $\mathbf{x} = [x_1, ..., x_n]^\top \in \mathbb{R}^n$, the expression $\mathbf{x} \le 0$ implies the set of component-wise inequalities $x_i \le 0, i \in \{1, ..., n\}$. $\mathbf{M} \succ 0$ and $\mathbf{M} \succeq 0$ imply that $\mathbf{M}$ is positive definite and semi-definite, respectively. Similarly, $\mathbf{M} \prec 0$ and $\mathbf{M} \succeq 0$ imply $-\mathbf{M} \succ 0$ and $-\mathbf{M} \succeq 0$, respectively. For a set $\mathcal{X}$, $\partial \mathcal{X}$ denotes its boundary. The subscript $m$ implies the signal is known by measurement. Finally, the superscript “ref” denotes the desired reference signal for the system to track (e.g., $\mathbf{u}^{\text{ref}}$).

## 2. Model representation and prerequisites

An overview of the general flatness-based control scheme is provided in Figure 1. One of the key contributions resides in the fact that the control synthesis (the blue block in Figure 1) is done in the flat output space (the red block) where the dynamics is linear. The reference for the linearized dynamics together with the feedback signals are given to a novel saturated (the blue block) to compute the stabilizing control action, $\mathbf{v}$. Then, such input, $\mathbf{v}$, in the flat output space, is transformed back to the real input, $\mathbf{u}$, via a flatness-based transformation (the green block). This input $\mathbf{u}$ is, in
turn, applied to the quadcopter (the white block in Figure 1), closing the loop for the entire scheme.

<figure><img src="image_2.png" /><figcaption>Figure 1: Flatness-based control scheme for the quadcopter position control.</figcaption></figure>

In this section, we refer to the white and green blocks of Figure 1 (with a particular attention to the constraints’ representation in the flat output space), while in Section 3 we will handle the blue block. Hence, in the following, the quadcopter model characterization in the flat output space is recapitulated. Then, some invariance and Lyapunov-based tools are also provided.

## 2.1. Quadcopter’s translational model and its characterization in the flat output space

Consider the position control problem of a quadcopter [29, 30]:

$$ \begin{bmatrix} \ddot{x} \\ \ddot{y} \\ \ddot{z} \end{bmatrix} = \begin{bmatrix} T(\cos \phi \sin \theta \cos \psi_m + \sin \phi \sin \psi_m) \\ T(\cos \phi \sin \theta \sin \psi_m - \sin \phi \cos \psi_m) \\ T \cos \phi \cos \theta - g \end{bmatrix} \triangleq h_{\psi}(\mathbf{u}), \quad (1) $$

where $x, y, z$ represent the positions of the quadcopter in the global frame’s three axes. $\psi_m$ is the measured yaw angle, which is considered as a known parameter for this system, $g$ is the gravitational acceleration while $\mathbf{u} = [T, \phi, \theta]^T \in \mathbb{R}^3$ collects the three inputs of the system, including, respectively, the thrust provided by the propellers, the control roll and pitch angles needed to steer the system towards its reference. Denote also $\mathcal{U} \subset \mathbb{R}^3$ the constraint set for $\mathbf{u} \in \mathbb{R}^3$, which is described as [23]:

$$ \mathcal{U} = \{ \mathbf{u} : 0 \le T \le T_{\text{max}}, |\phi| \le \phi_{\text{max}}, |\theta| \le \theta_{\text{max}} \} \quad (2) $$

where $T_{\text{max}} > 0$, $(\phi_{\text{max}}, \theta_{\text{max}}) \in (0; \pi/2)^2$ are constant bounds of the inputs.
**Assumption 1.** In the hierarchical scheme, the attitude controller is assumed to be capable of tracking the control roll and pitch angles $(\phi, \theta)$ given by the position control loop, i.e.:

$$ \phi_m = \phi \text{ and } \theta_m = \theta, \qquad (3) $$

where $\phi_m, \theta_m$ denote the actual value of the roll and pitch angles, respectively.

With Assumption 1, the control variables now become the thrust $T$, the roll $\phi$ and the pitch angles $\theta$. Meanwhile, the actual angles $\phi_m$ and $\theta_m$ will be controlled by the built-in attitude controller to track $\phi$ and $\theta$. The description of such an inner-loop employed for the experiments is also provided in Section 4.

The model (1) for thrust-propelled systems is known to be differentially flat [31, 9]. Namely, all state and input variables can be differentially parameterized by a fictitious output, called the flat output, and a finite number of its derivatives. Indeed, by choosing the flat output $\boldsymbol{\sigma}$ as:

$$ \boldsymbol{\sigma} = [\sigma_1, \sigma_2, \sigma_3]^\top \triangleq [x, y, z]^\top \in \mathbb{R}^3, \qquad (4) $$

we achieve the following flat representation:

$$ x = \sigma_1, y = \sigma_2, z = \sigma_3, \qquad (5a) $$

$$ T = \sqrt{\ddot{\sigma}_1^2 + \ddot{\sigma}_2^2 + (\ddot{\sigma}_3 + g)^2}, \qquad (5b) $$

$$ \phi = \arcsin((\ddot{\sigma}_1 \sin \psi_m - \ddot{\sigma}_2 \cos \psi_m)/T), \qquad (5c) $$

$$ \theta = \arctan((\ddot{\sigma}_1 \cos \psi_m + \ddot{\sigma}_2 \sin \psi_m)/(\ddot{\sigma}_3 + g)). \qquad (5d) $$

Then by exploiting (5b)-(5d), the linearizing feedback can be deduced as:

$$ T = \sqrt{v_1^2 + v_2^2 + (v_3 + g)^2}, \qquad (6a) $$

$$ \phi = \arcsin((v_1 \sin \psi_m - v_2 \cos \psi_m)/T), \qquad (6b) $$

$$ \theta = \arctan((v_1 \cos \psi_m + v_2 \sin \psi_m)/(v_3 + g)), \qquad (6c) $$

or in a more compact form:

$$ \boldsymbol{u} = \boldsymbol{\beta}_{\psi}(\boldsymbol{v}), \qquad (7) $$

where $\boldsymbol{v} = [v_1, v_2, v_3]^\top \in \mathbb{R}^3$ is the new input of the transformed system. Then, under the condition $v_3 \ge -g$ and the mapping (7), the system (1) yields:

$$ \dot{\boldsymbol{\xi}} = \boldsymbol{A}\boldsymbol{\xi} + \boldsymbol{B}\boldsymbol{v}, \qquad (8) $$
with $\xi = [\sigma^\top, \dot{\sigma}^\top]^\top = [x, y, z, \dot{x}, \dot{y}, \dot{z}]^\top \in \mathbb{R}^6$ and the matrices $A = \begin{bmatrix} 0_3 & I_3 \\ 0_3 & 0_3 \end{bmatrix}$,
$B = \begin{bmatrix} 0_3 \\ I_3 \end{bmatrix}$. With the dynamics (8) and the feedback law (6), in the next subsection, the constraints on $v$ will be analyzed to provide a context for the controller design developed later.

## 2.2. Input constraints description in the flat output space

Denote the constraint set for the input $v$ in (8) as:

$$ \mathcal{V} = \{v \in \mathbb{R}^3 : \beta_{\psi}(v) \in \mathcal{U} \text{ as in (2)}\}. \qquad (9) $$

Next, we show by a counterexample that $\mathcal{V}$ as in (9) is a non-convex set. More specifically, $\mathcal{V}$ contains the two points $v_\pm = h_\psi([T_{max}, \pm\phi_{max}, \pm\theta_{max}]^\top)$ but it does not contain their midpoint $0.5(v_- + v_+)$. Thus, the non-convexity hinders us from embedding the control (8) into any convex optimization based framework. Second, the description of $\mathcal{V}$ requires a constant update of the yaw angle $\psi_m$, which is usually time-varying in practical applications. This requirement is evidently demanding in terms of computation. For those reasons, $\mathcal{V}$ is rarely exploited in the literature. Indeed, as an alternative that is independent of the yaw angle, a convex subset of $\mathcal{V}$, denoted as $\mathcal{V}_c$, is frequently employed, as [32, 18]:

$$ \mathcal{V}_c = \left\{ v \in \mathbb{R}^3 : \begin{bmatrix} v_1^2 + v_2^2 + (v_3 + g)^2 - T_{max}^2 \\ v_1^2 + v_2^2 - (v_3^2 + g)^2 \tan^2 \epsilon_{max} \end{bmatrix} \le 0, \right. \\ \left. \epsilon_{max} \triangleq \min(\theta_{max}, \phi_{max}) \text{ and } v_3 \ge -g \right\}. \qquad (10) $$

In detail, the containment of $\mathcal{V}_c$ in (10), inside $\mathcal{V}$ as in (9) can be shown by considering the following inequalities (see Figure 2). Firstly, for the constrained angles, by applying Cauchy-Schwarz inequality in (6b) and (6c), we always have:

$$ |\sin \phi| \le \sin \epsilon(v) \text{ and } |\tan \theta| \le \tan \epsilon(v), \qquad (11) $$

where $\epsilon(v) \triangleq \arctan(\sqrt{(v_1^2 + v_2^2)/(v_3 + g)^2})$ interprets as a $v$-dependent angle upper bounding the roll $\phi$ and the pitch angle $\theta$. Then, with the assumption $|\theta|, |\phi| \le \pi/2$ as in (2) and by imposing $\epsilon(v) \le \epsilon_{max} \triangleq \min(\theta_{max}, \phi_{max})$, we can ensure $[|\theta| - \theta_{max}, |\phi| - \phi_{max}]^\top \le 0$ thanks to the local monotonicity of $\tan(\cdot)$ and $\sin(\cdot)$ functions. Secondly, for the constrained thrust $T$, it is straightforward that if $v_1^2 + v_2^2 + (v_3 + g)^2 - T_{max}^2 \le 0$ holds, we have $0 \le T \le T_{max}$. Therefore, $\mathcal{V}_c$ is a subset of $\mathcal{V}$. Moreover, the convexity of $\mathcal{V}_c$
<figure><img src="image_3.png" /><figcaption>Figure 2: Quadcopter’s input constraints in the flat output space with $\psi_m = 0$ (rad), $g = T_{max}/1.45 = 9.81$ (m/s²), $\phi_{max} = \theta_{max} = \pi/18$ (rad)).</figcaption></figure>

can be proven by showing that $V_c$ is the intersection of two convex sets: a ball described by $v_1^2 + v_2^2 + (v_3 + g)^2 - T_{max}^2 \le 0$ and a convex cone formed by two inequalities: $v_3 \ge -g$ and $v_1^2 + v_2^2 - (v_3 + g)^2 \tan^2 \epsilon_{max}$. An illustration of $V$ and $V_c$ is given in Figure 2.

**Remark 1.** Let us underline that, the new input constraint set $V_c$ does not appear in the standard form of the component-wise saturation on the signal magnitude [33, 34, 35], but in a convoluted nonlinear representation. This concern is the motivation to adapt the conventional saturation function towards a specified version particularized for the set $V_c$, which will be introduced in the next section.

The problem has now been reduced to the control of the linear dynamics (8) subject to a non-linear (but convex) constraint $V_c$ as in (10) in the flat output space (the blue block in Figure 1). Then, the real input $u$ is computed via the mapping $u = \beta_{\psi}(v)$ as in (7) and applied to the system.

## 2.3. Invariance and Lyapunov-based tools

In this subsection, we briefly recall some theoretical notions that will later be adopted in the control design.

**Proposition 1.** (Nagumo’s invariance principle [36]) Consider a dynamical system described as:

$$ \dot{\xi} = f(\xi, v), \qquad (12) $$
with $\xi \in \mathbb{R}^n$, $v \in \mathbb{R}^m$ representing the state and the input vector, respectively. Then, the ellipsoid set $\mathcal{E} = \{\xi \in \mathbb{R}^n : \xi^\top P \xi \le \varepsilon\}$ with $P \succ 0, \varepsilon > 0$ is an invariant set associated with the dynamics (12) if and only if:

$$ \xi_b^\top P f(\xi_b, v) \le 0, \forall \xi_b \in \partial \mathcal{E}. \qquad (13) $$

where $\xi_b \in \partial \mathcal{E}$ denotes the boundary points of the set $\mathcal{E}$.

Consider also the gradient-based control which can stabilize a linear system with an associated Lyapunov function.

**Proposition 2.** (Gradient-based control with a quadratic control Lyapunov function [36].)

Consider the following linear time-invariant (LTI) system

$$ \dot{\xi} = A\xi + Bv \qquad (14) $$

with $\xi \in \mathbb{R}^n$, $v \in \mathbb{R}^m$. Then, for any $\alpha > 0$ and $\gamma \ge 1$, by imposing a gradient-based control:

$$ v = -\gamma B^\top P \xi, \qquad (15) $$

satisfying the condition of $P \succ 0$ and:

$$ (A - BB^\top P)^\top P + P(A - BB^\top P) \le -\alpha P, \qquad (16) $$

then the quadratic positive definite function:

$$ V = \xi^\top P \xi, \qquad (17) $$

is a control Lyapunov function for the dynamics (14) associated with the control (15). As a result, this control setup ensures the exponential convergence of the system to the equilibrium point $\xi = \mathbf{0}_{n \times 1}, v = \mathbf{0}_{m \times 1}$.

Finally, recall the S-procedure transforming quadratic inequalities to linear matrix inequalities (LMI).

**Proposition 3.** (S-procedure [37]) With $\xi \in \mathbb{R}^n$, consider the two following scalar-valued quadratic functions:

$$ \mathcal{F}_i(\xi) = \xi^\top Q_i \xi + 2\mu_i^\top \xi + w_i, \quad i \in \{1, 2\}. \qquad (18) $$

Then, with $Q_i = Q_i^\top$, the two following statements are equivalent:
1. $\mathcal{F}_1(\xi) \geq 0 \Rightarrow \mathcal{F}_2(\xi) \geq 0$.

2. $\exists \tau \geq 0 \in \mathbb{R}$ such that:

$$ \begin{bmatrix} \mathbf{Q}_1 & \boldsymbol{\mu}_1 \\ \boldsymbol{\mu}_1^\top & w_1 \end{bmatrix} - \tau \begin{bmatrix} \mathbf{Q}_2 & \boldsymbol{\mu}_2 \\ \boldsymbol{\mu}_2^\top & w_2 \end{bmatrix} \succeq 0. \qquad (19) $$

With the above tools, in the subsequent section, a specialized saturated controller for the constraint set $\mathcal{V}_c$ is introduced, along with the stability and constraint satisfaction guarantees.

### 3. Explicit saturated control design

In this section, we propose a novel stabilizing constrained controller that provides noteworthy advantages for real-time embedded implementation.

#### 3.1. Controller characterization

With the goal to stabilize the system (8) under constraints (10), a saturated controller is introduced as follows.

**Proposition 4.** Consider a ball

$$ \mathcal{B}(\rho) = \{ \mathbf{v} \in \mathbb{R}^3 : \|\mathbf{v}\|_2^2 \leq \rho \} \text{ s.t } \mathcal{B}(\rho) \subseteq \mathcal{V}_c \text{ as in (10)}, \qquad (20) $$

some scalars $\alpha > 0, \gamma \geq 1$, a matrix $P \succ 0$ satisfying the condition (16), and the saturated controller defined as:

$$ \mathbf{v} = \text{sat}_{\lambda}(-\gamma \mathbf{B}^{\top} \mathbf{P} \boldsymbol{\xi}), \qquad (21) $$

with sat$_{λ}$(v) given by:

$$ \text{sat}_{\lambda}(\mathbf{v}) = \begin{cases} \mathbf{v} & \text{if } \mathbf{v} \in \mathcal{V}_c, \\ \lambda^*(\mathbf{v})\mathbf{v} & \text{if } \mathbf{v} \notin \mathcal{V}_c \end{cases} \qquad (22a) $$

$$ \text{with } \lambda^*(\mathbf{v}) = \underset{\lambda \mathbf{v} \in \mathcal{V}_c}{\arg \max} \lambda. \qquad (22b) $$

The following results hold for the system (8) under the controller (21):

1. The control action v as in (21) is contained in $\mathcal{V}_c$, i.e.,

$$ \mathbf{v} = \text{sat}_{\lambda}(-\gamma \mathbf{B}^{\top} \mathbf{P} \boldsymbol{\xi}) \in \mathcal{V}_c \forall \boldsymbol{\xi} \in \mathbb{R}^6, \gamma \geq 1. \qquad (23) $$
2. The following ellipsoidal set is rendered invariant:

$$ \mathcal{B}_P(\varepsilon) = \{ \boldsymbol{\xi} \in \mathbb{R}^6 : \|\boldsymbol{\xi}\|_P^2 \leq \varepsilon \}, \qquad (24) $$

where the largest $\varepsilon > 0$ is found by solving:

$$ \varepsilon^* = \arg\max_{\varepsilon} (\varepsilon) \qquad (25a) $$

$$ \text{s.t.} \quad \begin{bmatrix} \mathbf{P} & \mathbf{0}_{6 \times 1} \\ \mathbf{0}_{1 \times 6} & -\varepsilon \end{bmatrix} - \tau \begin{bmatrix} \mathbf{P} \mathbf{B} \mathbf{B}^{\top} \mathbf{P} & \mathbf{0}_{6 \times 1} \\ \mathbf{0}_{1 \times 6} & -\rho \end{bmatrix} \geq 0, \tau \geq 0. \qquad (25b) $$

Moreover, the system (8) achieves exponential stability with domain of attraction $\mathcal{B}_P(\varepsilon)$ as in (24).

PROOF. The constraint satisfaction follows immediately from condition (22). Next, to prove the stability, we show that the quadratic function associated with $\mathbf{P}$ is decreasing (its derivative is negative) inside $\mathcal{B}_P(\varepsilon)$. Then, the set $\mathcal{B}_P(\varepsilon)$ is indeed positively invariant. Hence, overall, once the system is inside the set, its stability and constraints are respected.

Firstly, let us proceed to prove the stability of the system inside $\mathcal{B}_P(\varepsilon)$ by explaining the choice of $\varepsilon$ as in (25). More specifically, by solving (25b), it can be shown that the containment of the state $\boldsymbol{\xi}$ inside the ellipsoid $\mathcal{B}_P(\varepsilon)$ implies that $-\mathbf{B}^{\top} \mathbf{P} \boldsymbol{\xi} \in \mathcal{B}(\rho)$ as in (20), i,e:

$$ \boldsymbol{\xi} \in \mathcal{B}_P(\varepsilon) \Rightarrow -\mathbf{B}^{\top} \mathbf{P} \boldsymbol{\xi} \in \mathcal{B}(\rho) \text{ as in (20).} \qquad (26) $$

Indeed, rewriting (26) explicitly yields:

$$ \|\boldsymbol{\xi}\|_P^2 \leq \varepsilon \Rightarrow \|\mathbf{B}^{\top} \mathbf{P} \boldsymbol{\xi}\|_2^2 \leq \rho \qquad (27) $$

Then by applying the $\mathcal{S}$-procedure, as in Proposition 3, the condition is (27) equivalently rewritten as:

$$ \exists \tau \geq 0 \text{ s.t.:} \quad \begin{bmatrix} \mathbf{P} & \mathbf{0}_{6 \times 1} \\ \mathbf{0}_{1 \times 6} & -\varepsilon \end{bmatrix} - \tau \begin{bmatrix} \mathbf{P} \mathbf{B} \mathbf{B}^{\top} \mathbf{P} & \mathbf{0}_{6 \times 1} \\ \mathbf{0}_{1 \times 6} & -\rho \end{bmatrix} \geq 0. \qquad (28) $$

Therefore, the conditions (26)-(27) can be assured by solving (25).

Moreover, from the definition (22), with a fixed choice of $\gamma \geq 1$, it can be noted that:

$$ \gamma \lambda^*(\gamma \mathbf{v}) \geq 1, \forall \mathbf{v} \neq \mathbf{0}_{3 \times 1}, \mathbf{v} \in \mathcal{V}_c. \qquad (29) $$
The explanation for (29) is twofold. First, given $v \in \mathcal{V}_c$, one always has:

$$ \gamma^{-1} \in \Delta \triangleq \{\lambda \in [0, 1] : \gamma\lambda v \in \mathcal{V}_c\} \quad (30) $$

Then, by definition (22), $\lambda^*(\gamma v) = \max \Delta$, hence $\lambda^*(\gamma v) \ge \gamma^{-1}$, leading to (29).

This property later will serve as an ingredient to show the non-positivity of the Lyapunov function's derivative inside the set $\mathcal{B}_P(\varepsilon)$ as in (24). More specifically, with the fact that $P$ satisfies the condition (16), consider the time-derivative of the Lyapunov function $V(\xi) \triangleq \xi^\top P \xi$ associated with the controller (21):

$$ \dot{V}(\xi) = 2\xi^\top P (A\xi + B\text{sat}_\lambda(-\gamma B^\top P\xi)). \quad (31) $$

Then, according to the definition of the saturation function in (22), the closed-loop behavior inside $\mathcal{B}_P(\varepsilon)$ is dependent on the two following cases:

**Unsaturated input:** If the control input (21) is not saturated, we have:

$$ v = \text{sat}_\lambda(-\gamma B^\top P \xi) = -\gamma B^\top P \xi. \quad (32) $$

Thus, the stability is straightforward, since the equation (31) now yields:

$$ \dot{V}(\xi) \le -\alpha V(\xi), \quad \forall \xi \in \mathcal{B}_P(\varepsilon). \quad (33) $$

**Saturated input:** If saturated, the control input (21) is computed as:

$$ v = \text{sat}_\lambda(-\gamma B^\top P \xi) = -\lambda^*(-\gamma B^\top P \xi)\gamma B^\top P \xi. \quad (34) $$

Then, taking into account (16), equation (31) becomes:

$$ \begin{aligned} \dot{V}(\xi) &= 2\xi^\top P A \xi + 2\xi^\top P B \text{sat}_\lambda(-\gamma B^\top P \xi), \\ &\le \xi^\top (-\alpha P + 2PBB^\top P)\xi + 2\xi^\top P B \text{sat}_\lambda(-\gamma B^\top P \xi), \\ &\le -\alpha \|\xi\|_P^2 - 2\|B^\top P \xi\|_2^2 (\gamma \lambda^*(-\gamma B^\top P \xi) - 1). \end{aligned} \quad (35) $$

Meanwhile, by replacing $v = -\gamma B^\top P \xi$ into the condition (29), one has$^2$:

$$ \gamma \lambda^*(-\gamma B^\top P \xi) \ge 1 \quad \forall \xi \in \mathcal{B}_P(\varepsilon) \setminus \{\mathbf{0}_{6 \times 1}\}. \quad (36) $$

<footnote>$^2$For $\xi = \mathbf{0}_{6 \times 1}$, we have back the unsaturated input case.</footnote>
As a consequence, (35) yields:

$$ \dot{V}(\xi) \leq -\alpha V(\xi) \quad \forall \xi \in \mathcal{B}_P(\varepsilon). \qquad (37) $$

Therefore, from (33) and (37), $V(\xi)$ is shown to be monotonically decreasing inside $\mathcal{B}_P(\varepsilon)$.

Secondly, the Nagumo's invariance condition for the set $\mathcal{B}_P(\varepsilon)$ can be evidently proven since, particularly in this case, the condition (13) yields:

$$ \xi_b^\top P (A\xi_b + B\text{sat}_\lambda(-\gamma B^\top P\xi_b)) = \dot{V}(\xi_b)/2 \leq 0, \qquad (38) $$

which holds $\forall \xi_b \in \partial\mathcal{B}_P(\varepsilon)$ because $\dot{V}(\xi) < 0 \forall \xi \in \mathcal{B}_P(\varepsilon)$ as proven previously for the stability. Furthermore, condition (36) along the invariance property ensures asymptotic stability with domain of attraction $\mathcal{B}_P(\varepsilon)$.

With the above demonstration, we summarize the synthesis procedure for the controller (21) as the following:

1. Choose $\rho > 0$, such that $\mathcal{B}(\rho) \subset \mathcal{V}_c$ as in (20);

2. Choose $\alpha > 0$ and solve (16) for a stabilizing symmetric matrix $P$;

3. With any choice $\gamma \geq 1$, the ellipsoid $\mathcal{B}(\varepsilon)$ as in (24) is guaranteed to be positively invariant.

**Remark 2.** It is apparent that no eigenvalue of the matrix $A$ in (8) has positive real part, and that the pair $(A, B)$ is stabilizable. Hence, as stated in [38, 39], for any fixed value of $\rho > 0$, there always exists a control law which makes the equilibrium points $(\xi = \mathbf{0}_{6 \times 1}, v = \mathbf{0}_{6 \times 1})$ globally asymptotically stable. Thus, with a proper choice of $\alpha$, one can enlarge the invariant set $\mathcal{B}_P(\varepsilon)$ defined in (24) arbitrarily just with the nominal controller satisfying:

$$ v = -B^\top P \xi \in \mathcal{B}(\rho) \forall \xi \in \mathcal{B}_P(\varepsilon). \qquad (39) $$

However, having a large domain of attraction along with the linear controller (39) produces very poor local performance, namely, close to the steady state. To tackle such drawback, we scaled up the nominal control (39) with $\gamma > 1$ and introduced the saturation function proposed in (22) to maintain continuity, invariance and constraint satisfaction.

**Remark 3.** It is also noteworthy that the saturation function $\text{sat}_\lambda(v)$ as in (22a) is Lipschitz continuous due to the convexity of $\mathcal{V}_c$ given in (10). This
property allows one to practically implement the controller without any chattering, singularities or discontinuous actuation. Moreover, input saturation in the literature is typically considered as a nonlinear component, which complicates the stability analysis and causes the system to behave unexpectedly or undesirably [40, 33]. Yet, in this paper, the novel saturation function for the quadcopter proposed in (22) helps to secure the constraint satisfaction. Even when the input $v$ is saturated on the boundary of the set $V_c$ as in (10), the stability is still maintained thanks to the controller characterization proposed previously.

Up to this point, a control strategy ensuring both exponential stability and constraint satisfaction has been presented. In the next subsection, we illustrate the synthesis procedure and provide some insights concerning the effect of different parameter values.

## 3.2. Invariant set characterization

As stated above, the proposed design process will be carried out to illustrate its feasibility and to analyze the parameters' effect. With the first step, instead of empirically choosing a small ball $\mathcal{B}(\rho)$ as in (20) inscribed in $V_c$, one candidate solution can be found as follows.

**Proposition 5.** The largest ball, of the form (20) and inscribed in $V_c$ as in (10), can be computed by solving the following semi-definite programming problem:

$$ \rho^* = \arg\max_{\rho} \rho \qquad (40) $$

$$ \text{s.t} \quad \left\{ \begin{array}{l} \begin{bmatrix} \mathbf{I}_3 & \mathbf{0}_{3 \times 1} \\ \mathbf{0}_{1 \times 3} & -\rho \end{bmatrix} - \tau_1 \begin{bmatrix} \mathbf{I}_3 & \bar{\mathbf{g}} \\ \bar{\mathbf{g}}^\top & \bar{\mathbf{g}}^\top \bar{\mathbf{g}} - T_{\max}^2 \end{bmatrix} \succeq 0 \end{array} \right. \qquad (41a) $$

$$ \begin{bmatrix} \mathbf{I}_3 & \mathbf{0}_{3 \times 1} \\ \mathbf{0}_{1 \times 3} & -\rho \end{bmatrix} - \tau_2 \begin{bmatrix} \mathbf{H} & \mathbf{H}^\top \bar{\mathbf{g}} \\ \bar{\mathbf{g}}^\top \mathbf{H} & \bar{\mathbf{g}}^\top \mathbf{H} \bar{\mathbf{g}} \end{bmatrix} \succeq 0 \qquad (41b) $$

$$ \begin{bmatrix} \mathbf{I}_3 & \mathbf{0}_{3 \times 1} \\ \mathbf{0}_{1 \times 3} & -\rho \end{bmatrix} - \tau_3 \begin{bmatrix} \mathbf{0}_3 & -\boldsymbol{\mu} \\ -\boldsymbol{\mu}^\top & -g \end{bmatrix} \succeq 0, \tau_i \ge 0, i \in \{1, 2, 3\} \qquad (41c) $$

with $\mathbf{H} = \text{diag}(1, 1, -\tan \epsilon_{\max})$, $\bar{\mathbf{g}} = [0, 0, g]^\top$ and $\boldsymbol{\mu} = [0, 0, 1/2]^\top$.

**PROOF.** Rewriting the formulation of $\mathcal{B}(\rho)$ in (20), we need the following condition to hold:

$$ \text{if } v \in \mathcal{B}(\rho) \Rightarrow v \in V_c \qquad (42) $$
or explicitly:

$$ \text{if } \|v\|_2^2 \le \rho \Rightarrow \begin{cases} (v + \bar{g})^\top (v + \bar{g}) \le T_{\max}^2 \\ (v + \bar{g})^\top H (v + \bar{g}) \le 0 \\ 2\mu^\top v + g \ge 0. \end{cases} \quad (43) $$

Then, by applying the S-procedure as in Proposition 3, we rewrite (43) as in (41), hence, completing the proof.

An illustration for the optimization problem (40) under the condition (41) is provided in Figure 3 and Table 1.

<figure><img src="image_4.png" /><figcaption>Figure 3: Numerical result for the solution of (40) - a ball maximally inscribed in $V_c$.</figcaption></figure>

Next, to ensure closed-loop stability, it is required to compute a stabilizing matrix $P$ satisfying condition (16). By setting $Q = P^{-1}$, pre and post multiplying the result with $Q$, the condition (16) yields:

$$ QA^\top + AQ - 2BB^\top \le -\alpha Q; \quad Q \succ 0, \quad (44) $$

which is, again, an LMI. Then finally, with the solution of (44), the invariant set as in (24) is characterized by $P = Q^{-1}$ and $\epsilon$ found in (25), finalizing the synthesis procedure.

For the sake of illustration, the simulations are depicted with different values for $\alpha > 0$ (see Figure 4). The numerical results were computed with
<table><caption>Table 1: Parameters for the invariant set construction</caption><thead><tr><th>Symbols</th><th>Values</th></tr></thead><tbody><tr><td>$T_{max} = 1.45g$</td><td>14.2245 m/s<sup>2</sup></td></tr><tr><td>$\theta_{max}, \phi_{max}, \epsilon_{max}$ as in (10)</td><td>0.1745 rad (10°)</td></tr><tr><td>$\rho$ in (20) computed with (40)-(41)</td><td>2.9019</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 4: The projection of $\mathcal{B}_P(\varepsilon)$ into the subspace $(x, \dot{x})$ with different choices of $\alpha$.</figcaption></figure>
Yalmip Toolbox [41] together with MATLAB 2021b. The specifications for the simulation are provided in Table 1.

As expected, the invariant set $\mathcal{B}_P(\varepsilon)$ can be enlarged by tuning down the value of $\alpha > 0$. The downside for this expansion is the reduction of the convergence speed, which can be interpreted from the differential inequalities (33) and (37). Therefore, a crucial point of concern to develop in the future studies is how to best trade the closed-loop performance for the extent of the region of attraction. Besides, computationally, the costliest step of the control implementation is the evaluation of $\lambda^*(v)$ as in (22b) when the input is saturated. One straightforward solution is to employ a reliable optimization solver during the flight to iteratively solve for the result. However, taking advantage of the representation of $\mathcal{V}_c$, an explicit solution can be derived for the controller implementation as shown in Section 3.3.

### 3.3. Implementation algorithm and simulation study

In view of the explicit computation of the proposed saturation function, recall its definition $\text{sat}_\lambda(v)$ from (22):

$$ \text{sat}_\lambda(v) = \begin{cases} v & \text{if } v \in \mathcal{V}_c \\ \lambda^*(v)v & \text{if } v \notin \mathcal{V}_c. \end{cases} \quad (45a) $$

$$ \text{where } \lambda^*(v) = \underset{\lambda v \in \mathcal{V}_c}{\arg\max} (\lambda). \quad (45b) $$

Since (45b) is a convex optimization problem, the value of $\lambda^*(v)$ can be derived directly from the candidate minimizers contained in $\mathcal{M}(v)$ (found from the Karush–Kuhn–Tucker necessary conditions [42]) as the following:

$$ \mathcal{M}(v) = \left\{ \frac{-g}{v_3}, \frac{a_0}{v_3}, \frac{-b_1 \pm \sqrt{b_1^2 - 4a_1c_1}}{2a_1}, \frac{-b_2 \pm \sqrt{b_2^2 - 4a_2c_2}}{2a_2} \right\} \quad (46) $$

$$ \text{with} \quad \begin{cases} a_0 = T_{\max} \cos \epsilon_{\max} - g, \quad a_1 = v_1^2 + v_2^2 - v_3^2 \tan^2(\epsilon_{\max}), \\ b_1 = -2 \tan^2(\epsilon_{\max}) v_3 g, \quad c_1 = - \tan^2(\epsilon_{\max}) g^2, \\ a_2 = v_1^2 + v_2^2 + v_3^2, \quad b_2 = 2v_3 g, \quad c_2 = g^2 - T_{\max}^2. \end{cases} $$

With this finite set of only 6 components, the computation of $\lambda^*(v)$ is reduced to a simple search for a maximum real number:

$$ \lambda^*(v) = \max \mathcal{M}(v) \cap (0; 1) \quad (47) $$
Note that $\mathcal{M}(\mathbf{v})$ as in (46) may contain complex numbers that have to be disregarded, since both inputs $\mathbf{v}$ and $\mathbf{u}$ are necessarily real. Then, (47) can be implemented by executing a low-cost search for a maximum value inside $\mathcal{M}(\mathbf{v}) \cap (0; 1)$ (e.g, a bubble sort). Briefly, the function (45a) (or (22a)) can be embedded with the following algorithm:

**Algorithm 1** Compute the function $\text{sat}_\lambda(\mathbf{v})$ as in (22).
**Input:** Variable $\mathbf{v}$, and constants $T_{\text{max}}, \epsilon_{\text{max}}, g$.
**Output:** The result of $\text{sat}_\lambda(\mathbf{v})$ as in (22)
**if** $\mathbf{v} \notin \mathcal{V}_c$ **then**
  Enumerate 6 components of $\mathcal{M}(\mathbf{v})$ as in (46);
  $\lambda^*(\mathbf{v}) \leftarrow \max\{\mathcal{M}(\mathbf{v}) \cap (0; 1)\}$;
  $\mathbf{v} \leftarrow \lambda^*(\mathbf{v})\mathbf{v}$;
**end if**
**return** $\mathbf{v}$

<table><caption>Table 2: Simulation result for $\mathcal{B}_P(\varepsilon)$ as in (24)</caption><tr><td>ε as in (24)</td><td>3.8692</td></tr><tr><td><b>P</b> = <b>Q</b><sup>-1</sup> as in (44)</td><td><b>[</b>0.2109<b>I</b><sub>3</sub> 0.2813<b>I</b><sub>3</sub><b>]</b><br/><b>[</b>0.2813<b>I</b><sub>3</sub> 0.7500<b>I</b><sub>3</sub><b>]</b></td></tr></table>

*Simulation study:* i) We first illustrate the invariance of $\mathcal{B}_P(\varepsilon)$ as in (24) regardless of the choice $\gamma \ge 1$ with the simulation of an origin-tracking problem. The trajectories' projection are given in Figure 5. In this simulation, the parameters in Table 1 were employed with $\alpha = 0.75$ to solve (44), resulting in an invariant ellipsoid $\mathcal{B}_P(\varepsilon)$ characterized by $\mathbf{P}$ and $\varepsilon$ given in Table 2; ii) With the initial point $\xi(0) = [-3.77, -0.46, -3.60, 0.94, -0.24, 2.31]^\top$, we validate also our Algorithm 1 with the numerical solution for the saturation function (22a) provided by the IPOPT solver [43] implemented with Yalmip MATLAB Toolbox. Furthermore, as a comparison in the computation cost, the flatness-based MPC (FB-MPC [27]) setup is adopted by replacing the controller (the blue block) in Figure 1 with the following MPC law:

$$ \arg\min_{\mathbf{v}_k, \dots, \mathbf{v}_{k+N_p-1}} \sum_{j=0}^{N_p-1} \|\xi_{k+j}\|_{\mathcal{Q}_{mpc}}^2 + \|\mathbf{v}_{k+j}\|_{\mathcal{R}_{mpc}}^2 \quad (48) $$
$$ \text{s.t: } \xi_{k+j+1} = A_d \xi_{k+j} + B_d v_{k+j}, \text{ and } v_{k+j} \in \tilde{V}_c $$

where the state are predicted with the discretized model of (8) by means of the Tustin discretization: $A_d = \begin{bmatrix} I_3 & t_s I_3 \\ 0_3 & I_3 \end{bmatrix}$, $B_d = \begin{bmatrix} t_s^2 I_3/2 \\ t_s I_3 \end{bmatrix}$ and the discretization time $t_s = 0.1(s)$; $\tilde{V}_c$ is the tight inner approximation of $V_c$ in (10), both provided in [27]; $N_p = 10$ is the prediction horizon and the weighting matrices are chosen as $Q_{mpc} = \text{diag}(50I_3, 5I_3)$, $R_{mpc} = 5I_3$.

<figure><img src="image_6.png" /><figcaption>Figure 5: The quadcopter's trajectory with different initial states (black points) implemented with different values of $\gamma$ (projected into $(x, \dot{x})$ subspace).</figcaption></figure>

Firstly, as can be seen from Figure 5, the control guarantees the invariance of the set $B_P(\varepsilon)$. The same confirmation can be made in all other subspaces during the simulation. Moreover, the validity of Algorithm 1 can also be seen in Figure 6 where the new input $v$ is always secured inside the constraint set $V_c$. Finally, the effect of the choice for $\gamma$ can also be distinguished from both Figure 7 and Figure 6. More specifically, when $\gamma = 1$, the controller is never saturated due to the property (29) as proven before. At the same time, when $\gamma > 1$, the damping effect appears to be increased together with $\gamma$. However, as can be seen from Figure 7, choosing a larger value for $\gamma$ does not necessarily imply a faster convergence, hence, giving us an indication not to abusively augment the value of this scalar while designing the controller.

**Remark 4.** *From a practical viewpoint, the most influential factor in choosing $\gamma$ relies on the well-known condition of the different frequency responses of the quadcopter's two hierarchical control loops [18, 11]. In detail, it is*
<figure><img src="image_7.png" /><figcaption>Figure 6: Input **u** with its constraints in solid black line (left) and **v** (right) simulated with different values of $\gamma$.</figcaption></figure>

required that the vehicle's inner loop (attitude control) responses considerably faster than its outer loop (position control), to decouple the two problems and solve them hierarchically. Therefore, although the effect cannot be seen in the simulation, it is important to choose the gain $\gamma$ in (21) appropriately to keep the convergence speed of the position control loop from exceeding the that of the attitude controller (Figure 11). In the next section, the undesirable behavior caused by the high gain $\gamma$ will also be demonstrated via the experiments.

<table><caption>Table 3: CPU time of the proposed algorithm w.r.t IPOPT solver and FB-MPC approach</caption><thead><tr><th>Control method</th><th>Average CPU time (s)</th><th>Min (s)</th><th>Max (s)</th><th>Optimization problem</th></tr></thead><tbody><tr><td>FBSC (IPOPT)</td><td>0.0162</td><td>0.0147</td><td>0.0851</td><td>QCQP</td></tr><tr><td>FBSC (Alg. 1)</td><td>1.63E-05</td><td>4.30E-06</td><td>0.0059</td><td>LP</td></tr><tr><td>FB-MPC [27]</td><td>0.015</td><td>0.0131</td><td>0.0248</td><td>QP</td></tr></tbody></table>

Besides, it is also noticeable from Figure 8 and Table 3 that Algorithm 1 generates the correct trajectory as the IPOPT solver does, but the consumed computation time for our solution is significantly reduced (–98% in average CPU time) compared to that of the solver. This complexity reduction can
<figure><img src="image_8.png" /><figcaption>Figure 7: Simulated trajectories in time with different values of $\gamma$, projected on the $z$ axis.</figcaption></figure>

be explained by the change from a Quadratically Constrained Quadratic Programming (QCQP) problem in (22) to a Linear Programming (LP) search for the maximum value described in Algorithm 1. This advantage can also be highlighted when it is compared to the Quadratic Programming (QP) problem (48) within the MPC settings. For more details, in Table 3, it can be seen that compared to the MPC, the FBSC solution with our algorithm significantly reduces the average computation time from $1.5 \times 10^{-2}$ (s) to $1.63 \times 10^{-5}$ (s). The advantage is also shown via the structural nature of the online program. With MPC, a processor equipped with an optimization solver is commonly required to solve the quadratic program online while with our explicit solution, the control can be computed via a closed-form representation as in (45) and can be effortlessly embedded in a microcontroller. These observations, hence, confirm not only the validity of our saturation component but also the computational advantage of the proposed solution (46).
<figure><img src="image_9.png" /><figcaption>Figure 8: The quadcopter’s trajectories and the computation time of Algorithm 1 in comparison with the IPOPT solver ($\gamma = 15$) and FB-MPC method [27].</figcaption></figure>

Before validating the proposed procedure, for completeness, let us address problem of bounding the angular rate of the quadcopter in the following part.

### 3.4. Discussion on bounding the angular rates

In practice, to make the connection with the tracking of the inner loop, it is critical to ensure angular velocity of the vehicle bounded. We will show that this problem can be described in the flat output space. A solution is also provided as the following.

In the literature, regarding the constraints on the rotational dynamics, the angular rates are typically expressed in the body frame as:

$$ \begin{bmatrix} \omega_x \\ \omega_y \\ \omega_z \end{bmatrix} = \begin{bmatrix} 1 & 0 & -\sin\theta \\ 0 & \cos\phi & \sin\phi\cos\theta \\ 0 & -\sin\phi & \cos\phi\cos\theta \end{bmatrix} \begin{bmatrix} \dot{\phi} \\ \dot{\theta} \\ \dot{\psi} \end{bmatrix}. \quad (49) $$

Then, the constraints on the angular rates can be taken into account with the following proposition.

**Proposition 6.** Consider the extended system:

$$ \frac{d}{dt} \underbrace{\begin{bmatrix} \xi \\ v \end{bmatrix}}_{\zeta} = \underbrace{\begin{bmatrix} A & B \\ 0_{3 \times 6} & 0_3 \end{bmatrix}}_{\hat{A}} \zeta + \underbrace{\begin{bmatrix} 0_{6 \times 3} \\ I_3 \end{bmatrix}}_{\hat{B}} \hat{v}, \quad (50) $$
where $\xi, v, A, B$ are from (8). Then, by ensuring:

$$ \mathbf{v} = [\mathbf{0}_{3 \times 6} \ \mathbf{I}_3]\zeta \in \mathcal{V}_c \text{ as in (10)}. \qquad (51a) $$

$$ v_3 = [0 \ 0 \ 1][\mathbf{0}_{3 \times 6} \ \mathbf{I}_3]\zeta \geq T_{min} - g, \qquad (51b) $$

$$ \|\hat{v}\| \leq \bar{\omega} T_{min} \qquad (51c) $$

the following constraints will be satisfied:

$$ \mathbf{u} \in \mathcal{U} \text{ as in (2)}, \qquad (52a) $$

$$ |\omega_x| \leq \bar{\omega}, |\omega_y| \leq \bar{\omega}, \qquad (52b) $$

$$ 0 < T_{min} \leq T \leq T_{max}, \qquad (52c) $$

where $0 < T_{min} < g$ is some lower bound imposed on the thrust $T$ in (5b) in addition to $T_{max}$ in (2). The angular rate w.r.t the $x, y$ axes in the body frame are denoted as $\omega_x, \omega_y$ and $\bar{\omega}$ is their magnitude constraint.

**PROOF.** The proof is deduced from the two following properties. First, by construction in (50) and (8), one has:

$$ \mathbf{v} = \ddot{\sigma}, \quad \hat{v} = \dot{v} = \ddot{\sigma}. \qquad (53) $$

Second, particularly for the considered quadcopter's model, the following relation holds [44, 32, 45, 46]:

$$ |\omega_q| \leq \frac{\|\ddot{\sigma}\|}{T} = \frac{\|\ddot{\sigma}\|}{\sqrt{\ddot{\sigma}_1^2 + \ddot{\sigma}_2^2 + (\ddot{\sigma}_3 + g)^2}}, q \in \{x, y\}. \qquad (54) $$

Hence, with the characterization in Section 2.2, (51a) implies (52a). Additionally, (51c) implies (52b) by means of the relation (54) and the bound:

$$ \frac{\|\ddot{\sigma}\|}{\sqrt{\ddot{\sigma}_1^2 + \ddot{\sigma}_2^2 + (\ddot{\sigma}_3 + g)^2}} \leq \frac{\|\ddot{\sigma}\|}{T_{min}}. $$

Finally, the thrust bounds in (52c) is ensured by imposing (51a)-(51b).

With this proposition, the constraints on the angular rate (52b) can be taken into account by stabilizing (50) while keeping the constraints in (51).

By choosing, for example, $\hat{v} = -\hat{K}\zeta$ with $\hat{K} = \hat{B}^\top \hat{P}$ for some stabilizing matrix $\hat{P}$ satisfying (44), a positive invariant ellipsoid $\hat{\mathcal{E}}$ in the state space $\zeta$ can be found by finding the set satisfying:

$$ \hat{\mathcal{E}} = \{\zeta : \|\zeta\|_{\hat{P}}^2 \leq \hat{\varepsilon}\} \subset \mathcal{Z} = \{\zeta : (51\text{a})-(51\text{b}), \|\hat{K}\zeta\| \leq \bar{\omega}T_{\text{min}}\}. \qquad (55) $$
<figure><img src="image_10.png" /><figcaption>Figure 9: Projection of $\hat{\mathcal{E}}$ in (55) on the $\xi$, $v$ axes with different tuning parameters.</figcaption></figure>

Note that $\mathcal{Z}$ is a convex set. A polytopic approximation $\tilde{\mathcal{Z}}$ of $\mathcal{Z}$ can be obtained by using polytopic inner approximations of $\mathcal{V}_c$ (See [27]) and the 3D ball in (51c). In this manner, the set $\hat{\mathcal{E}}$ can be found as the largest ellipsoid inscribed in the polytope $\tilde{\mathcal{Z}}$ [42]. For illustration, we adopt the LMI (44) for the extended system (50) with different choice of $\alpha$. Denote the corresponding solution from the LMI as $\hat{Q}$ and $\hat{P} = \hat{Q}^{-1}$. Then, the invariant sets characterized as in (55) are depicted in Figure 9. To show the constraint satisfaction and invariance properties, we simulate the system from 100 initial positions on the boundary of $\hat{\mathcal{E}}$ in (55) (see Figure 10).

<figure><img src="image_11.png" /><figcaption>Figure 10: Constraint satisfaction and stability for the extended system (50).</figcaption></figure>
# 4. Experimental validation

Next, to illustrate the applicability of our controller and highlight its advantages, experimental tests are carried out and compared with the results of existing methods in the literature. Analysis considering the practical choice of the $\gamma$ gain mentioned in Remark 4 will also be provided.

## 4.1. *Experimental setup and scenarios*

In the experiments, a nano-drone Crazyflie 2.1 platform will be employed. This aerial drone provides the access to its outer loop controller, which requires exactly the input $\mathbf{u} = [T, \phi, \theta]^\top$ as in (1) together with the desired angular rate for the yaw angle (for our scenario, this variable $\dot{\psi}^{\text{ref}}$ is considered as 0). The packet containing $\mathbf{u}$ and $\dot{\psi}^{\text{ref}} = 0$ will be computed via the proposed controller on a ground station computer and sent to the Crazyflie via a PA USB radio dongle. Note that, before being sent to the drone, the normalized thrust $T$ ($m/s^2$) needs to be converted to the thrust unit of the Crazyflie platform, ranging from 0 to 65535. Details on the conversion can be found on the documentation of Bitcraze AB team with the system identification fitting model [47]. The implementation is carried out using the *Commander* class provided in the API reference with Python 3.9. The computed parameters are sent to an onboard attitude controller, which calculates the desired torque necessary for the tracking. Together with such values, the PWM signals are computed accordingly via an appropriate rotor configuration [48, 49] and applied to the motors. The inner loop architecture of the experimental drone is depicted in Figure 11. Therein, for the attitude controller, two cascade control loops are integrated to stabilize the rotational dynamics [49]: i) the reference for the angular rates $\dot{\phi}^{\text{ref}}$, $\dot{\theta}^{\text{ref}}$ of the roll, pitch angles are computed via a PID controller stabilizing the errors between the measured and the desired values received from the position controller (the proposed outer-loop); ii) these rate references combined with $\dot{\psi}^{\text{ref}}$ will be tracked by another PID loop, which results in the desired torques in the direction of the roll, pitch and yaw angles. These torques and the thrust $T$ then are linearly transformed into PWM signals via a so-called *Control Mixer* specified for the “X”-shaped configuration of the nanodrone. Detailed description and operational parameters can be found in the technical report [48]. Other comparisons between different rotor configurations (e.g. “X”-shaped or “+”-shaped configuration) were studied in [50].
<figure><img src="image_12.png" /><figcaption>Figure 11: The quadcopter with a built-in inner loop (the white block in Figure 1).</figcaption></figure>

For the feedback sensors, an indoor optical motion capture system is employed. More specifically, 8 infrared cameras from Qualisys [51] are used to capture the position and orientation of a rigid body identified with a collection of reflective markers which are glued on the quadcopter. Then, based on the positions of the markers captured by the cameras, the motion of the drone are estimated by Qualisys Track Manager software and sent to the ground station computer in real time via TCP/IP. The cameras' capture rate is chosen at 120Hz for this experiment, while the controller's sampling time is chosen to be 0.075(s). An illustration of our experimental setup is given in Figure 12.

Next, to examine the performance of the controller (21), the following tracking scenarios are introduced as follows.

* Ref. 1: In this scenario, with a fixed choice of $P = Q^{-1}$ satisfying (44), the proposed controller will be employed to track a stationary point $\xi^{\text{ref}} = [0.3, 0.3, 0.8, 0, 0, 0]^\top$. We consider various values for $\gamma$ to analyze the impact of such constant on the quadcopter's tracking performance in real implementation.

* Ref. 2: For comparison, we adopt the time varying reference used for validation in [12] where a hybrid PID-fuzzy controller was used to stabilize the linearly approximated system around its equilibrium point.

* Ref. 3: We also consider the smooth trajectory generated via B-spline parameterization [52], and tracked in an MPC setup.

* Ref. 4: Furthermore, for comparison, we select also the circular trajectory applied in [19] where the authors applied the model predictive
<figure><img src="image_13.png" /><figcaption>Figure 12: LCIS lab experimental setup with optical motion capture system and nanodrones.</figcaption></figure>

control for the linear system (8) in the flat output space but with a conservative box-type constraints as opposed to $V_c$ as in (10). For this last scenario, the flat output given in (4) is parameterized as:

$$ \sigma_1^{\text{ref}}(t) = 0.5 \cos \omega t + 0.2; \sigma_2^{\text{ref}}(t) = 0.5 \sin \omega t $$

$$ \sigma_3^{\text{ref}}(t) = 0.3(m); \omega = 0.3\pi. $$

Finally, in all the scenarios, we simply employ the stabilizing matrix $P = Q^{-1} = \begin{bmatrix} 0.9766I_3 & 0.7813I_3 \\ 0.7813I_3 & 1.2500I_3 \end{bmatrix}$ for the controller (21) as a solution of (44) while choosing $\alpha = 1.25$, so that the resulted invariant set $B_P(\varepsilon)$ as in (24) covers the drone's initial position. The video of the experiment is available at https://youtu.be/0YcEM280gts.

## 4.2. Experimental results and discussion

Overall, as can be seen in Figure 13, 16 and Table 4, all four trajectories were tracked properly with the root-mean-square (RMS) of the tracking errors under 10 centimeters. Moreover, the input constraints are always respected in all scenarios and the computational cost is considerably low
<table><caption>Table 4: Experiment specifications and results</caption><thead><tr><th></th><th>Ref. 1</th><th>Ref. 2</th><th>Ref. 3</th><th>Ref. 4</th></tr></thead><tbody><tr><td>Average computation time (s) ×10<sup>-4</sup></td><td>1.164</td><td>1.321</td><td>1.327</td><td>1.323</td></tr><tr><td>RMS of tracking error (cm)</td><td>10.03</td><td>1.99</td><td>2.57</td><td>3.89</td></tr><tr><td>γ as in (21)</td><td>5</td><td>4.5</td><td>4.5</td><td>4.5</td></tr><tr><td>α; <b>P</b> as in (21)</td><td colspan="4">1.25; <b>[0.9766I<sub>3</sub> 0.7813I<sub>3</sub><br>0.7813I<sub>3</sub> 1.2500I<sub>3</sub>]</b></td></tr><tr><td>Sampling time (s)</td><td colspan="4">0.075</td></tr></tbody></table>

compared to the other constraint-handling controllers with stability guarantees found in the literature [19, 23, 28]. Further details on the experimental results and our proposed controller performance are given in the following.

<figure><img src="image_14.png" /><figcaption>Figure 13: Ref. 1 (blue-dashed line) tracked with different values of $\gamma$.</figcaption></figure>

Firstly, the tracking result for Ref. 1 is depicted in Figure 13 while the corresponding input signals in the original space and in the flat output space are given in Figure 14 and 15, respectively. As expected from the simulation in Figure 7, after increasing to a certain value of $\gamma > 1$ in (21), not only does the settling time not decrease, but also an unwanted oscillating effect appears
(Figure 13). While this effect is not apparent in the simulation, it can be explained as the outcome of the high gain controller that awakes the “parasitic” dynamics existing within the vehicles (e.g. the rotational/attitude dynamics in the hierarchical scheme). Therefore, it is practically verified that $\gamma$ does not need to be chosen excessively large to achieve a favorable performance. Besides, from Figure 14 and 15, the effectiveness of the proposed saturation function is again validated via the restraint $v \in V_c$ given in (10), or consequently, the constraint satisfaction $u \in U$ given in (2).

<figure><img src="image_15.png" /><figcaption>Figure 14: Input $u$ with different values of $\gamma$ with the set-point tracking problem and its bounds $U$ as in (2) (black dashed line).</figcaption></figure>

Next, it is also notable that the average computation time rests approximately at 0.1 milliseconds (see Table 4). In terms of computing time, this result surpasses previously mentioned controllers in the literature (around 5 to 50 milliseconds and above), which also guarantee stability and constraint satisfaction [19, 23, 28, 24]. This satisfactory performance comes from the fact that we provide an explicit formula of the controller from computing the feedback input to saturating such signal via Algorithm 1. Although in this experiment, all the control computation is executed remotely in the ground station computer, these results appear promising for embedding the controller directly on the drone's microcontroller thanks to the simplicity of the proposed scheme.
<figure><img src="image_16.png" /><figcaption>Figure 15: Input **v** in the flat output space, with different values of $\gamma$ in Ref. 1 tracking problem and its constraint set $V_c$ as in (10).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16: Position tracking with different trajectories.</figcaption></figure>

Regarding the performance, the tracking errors provided by the proposed method can be considered commensurate to the FB-MPC, the well-known piece-wise affine MPC (where the model is approximated along the trajectory) [27], or the fuzzy-based method proposed in [12].
<figure><img src="image_18.png" /><figcaption>Figure 17: Input signal **u** in the tracking problem for Ref. 2, Ref. 3 and Ref. 4 with its bounds (black solid line).</figcaption></figure>

Finally, regarding the shortcomings of the proposed scheme, it is noted
that, in the model (1) used for control design, the aerodynamics forces are
neglected. Moreover, although validated via the experiments, it is vital to
analyze the robustness of the systems against external disturbances as well
as uncertainty. Next, as mentioned above, the selection of the high gain
$\gamma$ plays also an important role in the overall stability of the system. For
these analyses, we refer the readers towards our ensuing study [53] which
characterizes the robustness and deals with the gain selection problem via
an adaptive scheme.

**5. Conclusion**

In this work, an explicit saturated controller for the quadcopter system
was presented with the guarantee for both constraint satisfaction and stabil-
ity. Its practical viability was also validated in simulations and experiments.
The controller was constructed by taking advantage of the flatness-based lin-
earization of the system and its corresponding distorted input constraints in
the new coordinates. Hence, on one hand, its validated applicability opens
the door to further developments in computationally low-cost controllers for quadcopter vehicles. Future research will focus on the controller's robustification and optimality when it comes to the selection of feedback gains $\gamma$ and $P$ rather than the current feasible choices. On the other hand, with the linear description in the flat output space, issues concerning cooperative or consensus formation of multiple drones are dynamically connected with the well-developed frameworks of linear multi-agent systems. Hence, the constrained coordination and control of thrust-propelled vehicles are promising research directions to delve into.

## References

[1] A. V. Borkar, V. S. Borkar, A. Sinha, Aerial monitoring of slow moving convoys using elliptical orbits, European Journal of Control 46 (2019) 90–102.

[2] A. A. Nasir, H. D. Tuan, T. Q. Duong, H. V. Poor, UAV-enabled communication using NOMA, IEEE Transactions on Communications 67 (7) (2019) 5126–5138. doi:10.1109/TCOMM.2019.2906622.

[3] Y. Liu, Q. Wang, H. Hu, Y. He, A novel real-time moving target tracking and path planning system for a quadrotor UAV in unknown unstructured outdoor scenes, IEEE Transactions on Systems, Man, and Cybernetics: Systems 49 (11) (2018) 2362–2372.

[4] A. Zagaris, H. G. Kaper, T. J. Kaper, Fast and slow dynamics for the computational singular perturbation method, Multiscale Modeling & Simulation 2 (4) (2004) 613–638.

[5] S. Bertrand, N. Guénard, T. Hamel, H. Piet-Lahanier, L. Eck, A hierarchical controller for miniature vtol uavs: Design and stability analysis using singular perturbation theory, Control Engineering Practice 19 (10) (2011) 1099–1108.

[6] T. P. Nascimento, M. Saska, Position and attitude control of multi-rotor aerial vehicles: A survey, Annual Reviews in Control 48 (2019) 129–146.

[7] N. Cao, A. F. Lynch, Inner-outer loop control for quadrotor uavs with input and state constraints, IEEE Transactions on Control Systems Technology 24 (5) (2016) 1797–1804. doi:10.1109/TCST.2015.2505642.
[8] P. Foehn, E. Kaufmann, A. Romero, R. Penicka, S. Sun, L. Bauersfeld, T. Laengle, G. Cioffi, Y. Song, A. Loquercio, et al., Agilicious: Open-source and open-hardware agile quadrotor for vision-based flight, Science robotics 7 (67) (2022) eabl6259.

[9] S. Sun, A. Romero, P. Foehn, E. Kaufmann, D. Scaramuzza, A comparative study of nonlinear MPC and differential-flatness-based control for quadrotor agile flight, IEEE Transactions on Robotics (2022).

[10] E. Tal, S. Karaman, Accurate tracking of aggressive quadrotor trajectories using incremental nonlinear dynamic inversion and differential flatness, IEEE Transactions on Control Systems Technology 29 (3) (2020) 1203-1218.

[11] A. Intwala, Y. Parikh, A review on vertical take off and landing (vtol) vehicles, International Journal of Innovative Research in Advanced Engineering (IJIRAE) 2 (2) (2015) 187-191.

[12] F. Santoso, M. A. Garratt, S. G. Anavatti, Hybrid PD-fuzzy and pd controllers for trajectory tracking of a quadrotor unmanned aerial vehicle: Autopilot designs and real-time flight tests, IEEE Transactions on Systems, Man, and Cybernetics: Systems 51 (3) (2019) 1817-1829.

[13] C. Zhu, J. Chen, M. Iwasaki, H. Zhang, Event-triggered deep learning control of quadrotors for trajectory tracking, IEEE Transactions on Industrial Electronics (2023).

[14] L. Derafa, A. Benallegue, L. Fridman, Super twisting control algorithm for the attitude tracking of a four rotors uav, Journal of the Franklin Institute 349 (2) (2012) 685-699.

[15] N. A. Ofodile, M. C. Turner, Anti-windup design for input-coupled double integrator systems with application to quadrotor uav's, European Journal of Control 38 (2017) 22-31.

[16] J. Levine, Analysis and control of nonlinear systems: A flatness-based approach, Springer Science & Business Media, 2009.

[17] D. Ma, Y. Xia, G. Shen, Z. Jia, T. Li, Flatness-based adaptive sliding mode tracking control for a quadrotor with disturbances, Journal of the Franklin Institute 355 (14) (2018) 6300-6322.
[18] N. T. Nguyen, I. Prodan, L. Lefèvre, Flat trajectory design and tracking with saturation guarantees: a nano-drone application, International Journal of Control 93 (6) (2020) 1266-1279.

[19] M. Greeff, A. P. Schoellig, Flatness-based model predictive control for quadrotor trajectory tracking, in: 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems, IEEE, 2018, pp. 6740–6745.

[20] H. Lu, C. Liu, L. Guo, W.-H. Chen, Constrained anti-disturbance control for a quadrotor based on differential flatness, International Journal of Systems Science 48 (6) (2017) 1182–1193.

[21] M. W. Mueller, R. D'Andrea, A model predictive controller for quadcopter state interception, in: 2013 European Control Conference (ECC), IEEE, 2013, pp. 1383–1389.

[22] N. T. Nguyen, I. Prodan, L. Lefèvre, A stabilizing NMPC design for thrust-propelled vehicles dynamics via feedback linearization, in: 2019 American Control Conference (ACC), IEEE, 2019, pp. 2909–2914.

[23] N. T. Nguyen, I. Prodan, L. Lefèvre, Stability guarantees for translational thrust-propelled vehicles dynamics through NMPC designs, IEEE Transactions on Control Systems Technology 29 (1) (2020) 207–219.

[24] M. A. Gomaa, O. De Silva, G. K. Mann, R. G. Gosine, Computationally efficient stability-based nonlinear model predictive control design for quadrotor aerial vehicles, IEEE Transactions on Control Systems Technology 31 (2) (2022) 615–630.

[25] A. R. Teel, Global stabilization and restricted tracking for multiple integrators with bounded controls, Systems & control letters 18 (3) (1992) 165–171.

[26] H. Xie, D. Veitch, Nested saturation control of multiple vector integrators and its application to motion control of uavs, International Journal of Robust and Nonlinear Control 30 (1) (2020) 246–265.

[27] H. T. Do, I. Prodan, Indoor experimental validation of MPC-based trajectory tracking for a quadcopter via a flat mapping approach, 2023 European Control Conference (ECC) (2023).
[28] N. T. Nguyen, I. Prodan, Stabilizing a multicopter using an nmpc design with a relaxed terminal region, IFAC-PapersOnLine 54 (6) (2021).

[29] S. Wang, A. Polyakov, G. Zheng, Quadrotor stabilization under time and space constraints using implicit pid controller, Journal of the Franklin Institute 359 (4) (2022) 1505-1530.

[30] L. Amezquita-Brooks, E. Liceaga-Castro, M. Gonzalez-Sanchez, O. Garcia-Salazar, D. Martinez-Vazquez, Towards a standard design model for quad-rotors: A review of current models, their accuracy and a novel simplified model, Progress in Aerospace Sciences 95 (2017) 1-23.

[31] N. T. Nguyen, I. Prodan, F. Stoican, L. Lefevre, Reliable nonlinear control for quadcopter trajectory tracking through differential flatness, IFAC-PapersOnLine 50 (1) (2017) 6971-6976.

[32] V. Freire, X. Xu, Flatness-based quadcopter trajectory planning and tracking with continuous-time safety guarantees, IEEE Transactions on Control Systems Technology (2023).

[33] T. Hu, Z. Lin, Control systems with actuator saturation: analysis and design, Springer Science & Business Media, 2001.

[34] F. Blanchini, Feedback control for linear systems with state and control bounds in the presence of disturbance, IEEE Trans. Aut. Control 35 (1990) 1131-1135.

[35] R. Wang, J. Liu, Trajectory tracking control of a 6-dof quadrotor uav with input saturation via backstepping, Journal of the Franklin Institute 355 (7) (2018) 3288-3309.

[36] F. Blanchini, S. Miani, Set-theoretic methods in control, Vol. 78, Springer, 2008.

[37] S. Boyd, L. El Ghaoui, E. Feron, V. Balakrishnan, Linear matrix inequalities in system and control theory, SIAM, 1994.

[38] E. D. Sontag, H. J. Sussmann, Nonlinear output feedback design for linear systems with saturating controls, in: 29th IEEE Conference on Decision and Control, IEEE, 1990, pp. 3414-3416.
[39] A. Saberi, A. A. Stoorvogel, P. Sannuti, Control of linear systems with regulation and input constraints, Springer Science & Business Media, 2012.

[40] T. Hu, Z. Lin, B. M. Chen, Analysis and design for discrete-time linear systems subject to actuator saturation, Systems & control letters 45 (2) (2002) 97–112.

[41] J. Lofberg, Yalmip : a toolbox for modeling and optimization in matlab, in: 2004 IEEE International Conference on Robotics and Automation, 2004, pp. 284–289. doi:10.1109/CACSD.2004.1393890.

[42] S. Boyd, S. P. Boyd, L. Vandenberghe, Convex optimization, Cambridge university press, 2004.

[43] A. Wächter, L. T. Biegler, On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming, Mathematical programming 106 (1) (2006) 25–57.

[44] M. Hehn, R. D'Andrea, Real-time trajectory generation for quadcopters, IEEE Transactions on Robotics 31 (4) (2015) 877–892.

[45] N. T. Nguyen, Reliable hierarchical control for multicopter systems, Ph.D. thesis, Université Grenoble Alpes (2019).

[46] V. Marguet, C. K. Dinh, F. Stoican, I. Prodan, Indoor formation motion planning using b-splines parametrization and evolutionary optimization, Control Engineering Practice 152 (2024) 106066.

[47] J. Förster, System identification of the crazyflie 2.0 nano quadrocopter, B.S. thesis, ETH Zurich (2015).

[48] C. Luis, J. L. Ny, Design of a trajectory tracking controller for a nanoquadcopter, Tech. rep., Mobile Robotics and Autonomous Systems Laboratory, Polytechnique Montreal (2016).

[49] W. Giernacki, M. Skwierczyński, W. Witwicki, P. Wroński, P. Kozier- ski, Crazyflie 2.0 quadrotor as a platform for research and education in robotics and control engineering, in: 2017 22nd International Confer- ence on Methods and Models in Automation and Robotics, IEEE, 2017, pp. 37–42.
[50] M. Idrissi, M. Salami, F. Annaz, A review of quadrotor unmanned aerial vehicles: applications, architectural design and control algorithms, Journal of Intelligent & Robotic Systems 104 (2) (2022) 22.

[51] D. Senior, Qualisys track manager: User manual, Laboratory Memorandum (2004).

[52] I. Prodan, F. Stoican, C. Louembet, Necessary and sufficient LMI conditions for constraints satisfaction within a B-spline framework, in: IEEE 58th Conference on Decision and Control, IEEE, 2019, pp. 8061–8066.

[53] H.-T. Do, F. Blanchini, I. Prodan, Enhancements on a saturated control for stabilizing a quadcopter: Adaptive and robustness analysis in the flat output space, in: 2023 62nd IEEE Conference on Decision and Control (CDC), 2023, pp. 739–744. doi:10.1109/CDC49753.2023.10384084.