Davide Invernizzi
email: davide.invernizzi@polimi.it

Marco Lovera
email: marco.lovera@polimi.it

Luca Zaccarian
email: zaccarian@laas.fr

Integral ISS-based cascade stabilization for vectored-thrust UAVs

Keywords:

We address stabilization of vectored-thrust Unmanned Aerial Vehicles (UAVs): a challenging task due to the peculiar nonlinear underactuated dynamics. According to a wellestablished approach based on the selection of suitable error variables, the error dynamics is described as a pseudo cascade connection where the attitude subsystem indirectly stabilizes the position dynamics. Unlike existing works, we address stability of this cascade using integral Input to State Stability (iISS). With this extension it is possible to employ a quasi time-optimal control law to achieve global asymptotic stabilization of any desired position, thus improving the transient performance with respect to existing control designs. A simulation example shows the performance improvement in comparison with a nested saturations stabilizer.

I. INTRODUCTION

V ECTORED-THRUST Unmanned Aerial Vehicles (UAVs) are platforms endowed with a propulsive system that can deliver a torque in any direction of the aircraft frame but a force only along a fixed axis: the thrust vector. Due to this underactuation, the position dynamics can be stabilized only through suitable attitude motions. There is a significant body of literature, well surveyed in [1], dealing with the stabilization problem for vectored-thrust UAVs as formulated in this paper, which is challenging also because it evolves on a nonlinear manifold corresponding to the set of rigid body motions. Early works on the subject date back to [START_REF] Frazzoli | Trajectory tracking control design for autonomous helicopters using a backstepping algorithm[END_REF] and [START_REF] Hamel | Dynamic modelling and configuration stabilization for an x4-flyer[END_REF], while recent approaches are presented in [START_REF] Roza | A class of position controllers for underactuated VTOL vehicles[END_REF], [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using innerouter loop control paradigms[END_REF], [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF].

The goal of this paper is to propose a control design which overcomes the limitations, in terms of transient performance, of modern strategies developed to globally stabilize the position of underactuated UAVs. In particular, the control designs proposed in [START_REF] Roza | A class of position controllers for underactuated VTOL vehicles[END_REF], [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using innerouter loop control paradigms[END_REF] as well as in our previous works [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF], [START_REF]Dynamic attitude planning for trajectory tracking in underactuated UAVs[END_REF] ensure stability of the intrinsic cascade via a reduction approach [START_REF] Maggiore | Reduction theorems for hybrid dynamical systems[END_REF] assessing the necessary property of "boundedness of solutions" by requiring the position stabilizer to possess Input to State Stabiliy (ISS) with respect to small inputs [START_REF] Chaillet | Combining iISS and ISS with respect to small inputs: the strong iISS property[END_REF]. This small-signal ISS property is typically guaranteed by employing nested saturations-based stabilizers [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using innerouter loop control paradigms[END_REF]: although very robust, this solution is known to yield poor transient performance [START_REF] Marchand | Global stabilization of multiple integrators with bounded controls[END_REF].

In this paper, we propose a hierarchical control design that extends the class of allowable position stabilizers to solutions having less stringent robustness requirements, thereby allowing for a more aggressive behavior. In particular, after introducing the UAV dynamics and the control problem in Section II, we borrow tools from stability analysis of integral ISS cascades [START_REF] Sontag | Comments on integral variants of ISS[END_REF], [START_REF] Arcak | A unifying integral ISS framework for stability of nonlinear cascades[END_REF], and we propose in Section III a control architecture inducing an autonomous pseudo-cascade structure of the closed-loop error dynamics while giving freedom in the selection of the attitude and position stabilizers and in the shape of the interconnection term, as long as they satisfy certain properties always ensuring boundedness of solutions. Then in Section IV, we select these three basic components using high performance strategies. In particular, we suggest using an exponential attitude stabilizer [START_REF] Lee | Exponential stability of an attitude tracking control system on so (3) for large-angle rotational maneuvers[END_REF] and a Quasi Time-Optimal (QTO) stabilizer [START_REF] Forni | A family of global stabilizers for quasi-optimal control of planar linear saturated systems[END_REF] and we prove that this choice makes the UAV position dynamics integral ISS. The QTO stabilizer is the prototypical example of aggressive solutions, achieving excellent transient performance, which cannot be included in existing designs, notably [START_REF] Roza | A class of position controllers for underactuated VTOL vehicles[END_REF], [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using innerouter loop control paradigms[END_REF], [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF], since proving (the previously required) small-signal ISS property is hard if not impossible. A simulation example discussed in Section V highlights the advantage of the proposed solution with respect to a nested saturations approach in a scenario characterized by a significant initial configuration error.

Notation. R (R >0 , R ≥0 ) denotes the set of (positive, nonnegative) real numbers, R n denotes the n-dimensional Euclidean space and R m×n the set of m × n real matrices. The canonical basis in R n is denoted as e i := (0, . . . , 1, . . . , 0) for i ∈ {1, ..., n} and the identity matrix in R n×n is denoted as

I n := [e 1 • • • e i • • • e n ]. Given x = (x 1 , . . . , x n ) ∈ R n , x := x 2 1 + . . . + x 2
n is the Euclidean norm while for a matrix A ∈ R n×n , A F := tr(A T A) is the Frobenius norm and skew(A) := A-A T 2 is the skew-symmetric part of A. The set SO(3) := {R ∈ R 3×3 : R T R = I 3 , det(R) = 1} denotes the thirdorder Special Orthogonal group while S n := {q ∈ R n+1 : q = 1} denotes the n-dimensional unit sphere. where the normalized distance with respect to I 3 , induced by the Frobenius norm, is denoted as 3 and × represents the cross product in R 3 . The inverse of the hat map is the vee map, denoted as (•) ∨ : so(3) → R 3 . We use standard comparison functions from [START_REF] Khalil | Nonlinear systems, 3rd Edition[END_REF]: function α : R ≥0 → R ≥0 is of class K if it is zero at zero, strictly increasing, and continuous.

R SO(3) := 1 8 R -I 3 F = 1 4 tr(I 3 -R) ∈ [0, 1]. Given ω ∈ R 3 , the hat map • : R 3 → so(3) := {Ω ∈ R 3×3 : Ω = -Ω T } is such that ωy = ω × y, ∀y ∈ R
It is of class K ∞ if it is also unbounded. β : R ≥0 × R ≥0 → R ≥0 is of class KL if it is of class K in
the first argument and nonincreasing and converging to zero as its second argument tends to +∞.

II. DYNAMICAL MODEL AND CONTROL PROBLEM

The configuration of a UAV is globally and uniquely described by the pair (R, x) ∈ SO(3) × R 3 , where R and x are, respectively, the rotation matrix and the position vector of a body-fixed frame

F B = (O B , {b 1 , b 2 , b 3 }) with respect to an inertial frame F I = (O I , {i 1 , i 2 , i 3 }).
A reasonable model for control design is described by the following set of differential equations [1]:

ẋ = v m v = -mge 3 + R f c (1) Ṙ = R ω J ω = -ω × Jω = -ωJω + τ c , (2)
where ω ∈ R 3 is the body angular velocity, v ∈ R 3 is the translational velocity of the center of mass, resolved in F I , m ∈ R >0 and J = J T ∈ R 3×3 >0 are the mass and the body inertia matrix of the UAV, respectively, and g = 9.81m/s 2 represents the gravitational acceleration. Finally, f c ∈ R 3 is the control force and τ c ∈ R 3 the control torque, both resolved in F B . The propulsive system of vectored-thrust UAVs allows to produce a (bounded) control force directed only along the positive direction of b 3 , i.e., the components of f c in the body frame must satisfy

f c = T c e 3 , 0 < T c ≤ T M ∈ R >0 , (3)
where T c is the overall thrust that can be assigned by properly modulating the thrust delivered by each propeller. It is assumed that the control torque τ c spans R 3 , i.e., that the rotational dynamics ( 2) is fully In this paper we exploit a cascade-based design for the attitude dynamics and consider ω c := ω as a control input, under the assumption that the attitude dynamics in equation ( 2) can be made sufficiently fast in tracking ω c by a suitable selection of τ c (for a discussion about this approximation, see [1]). In this case, the model for control design reduces to:

ẋ = v (4) m v = -mge 3 + T c Re 3 (5) Ṙ = R ωc (6)
where the attitude motion (6) is independent of the translational one (4), [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using innerouter loop control paradigms[END_REF]. This structure is exploited in the hierarchical control strategy presented in the next section. This work focuses on the stabilization of the underactuated model ( 4)- [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF]. According to this model, it is either possible to stabilize a desired attitude and altitude or a desired position and rotation about the gravity axis e 3 . This paper addresses the latter problem under the following standard assumption.

Assumption 1: 1) the desired trajectory is a constant reference (x d , 0, R d ) ∈ R 3 × R 3 × SO(3) such that R d e 3 = e 3 ; 2) the maximum available thrust T M is larger than the UAV weight, i.e., T M > mg; 3) the state of the system (x, v, R) is available for feedback. While the first part of the assumption is necessary to account for the platform underactuation when constraint (3) is included in the dynamical model, the second part is necessary control authority requirement to have a solvable control problem. In view of Assumption 1, the control problem can be formally stated as follows.

Problem 1: Consider the dynamical model described by equations ( 4)- [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF] and assume that the control thrust T c is bounded as 0 < T c ≤ T M for a strictly positive scalar T M . Find a state feedback control law for u :

= (T c , ω c ) ∈ (0, T M ] × R 3 such that any constant set point (x d , 0, R d ) ∈ R 3 × R 3 × SO(3)
satisfying Assumption 1, is asymptotically stable.

III. CONTROL DESIGN AND STABILITY ANALYSIS

This section is devoted to presenting the control design and the corresponding stability analysis. In [START_REF] Roza | A class of position controllers for underactuated VTOL vehicles[END_REF], [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using innerouter loop control paradigms[END_REF], [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF] a smallsignal ISS property was the basic requirement for the position stabilizer to ensure stability of the closed-loop system. Here we relax this requirement and enlarge the class of allowable position stabilizers towards solutions capable of superior performance. To this end, we invoke recent results on stabilization of nonlinear cascades, based on integral ISS arguments (see equation ( 24) below, and [START_REF] Sontag | Comments on integral variants of ISS[END_REF] for additional details about iISS). The formulation that we propose leaves the door open for different selections of stabilizers of the nominal position and attitude dynamics, therefore parametrizing a family of solutions to Problem 1.

A. Control law and closed-loop dynamics

We propose an attitude planner paradigm, similar to the one presented in our previous work [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF], to tackle the system underactuation. Let us first introduce the stabilization error

e x := x -x d , e v := v, R e := R T p (e x , e v , R d )R, (7)
where R p : R 3 × R 3 × SO(3) → SO(3) is the reference attitude provided by the attitude planner (see Figure 1) defined as:

is guaranteed to never vanish by design. Based on ( 8)-( 9) we select the inputs of dynamics ( 4)-( 6) as

T c = c(R e )e T 3 R T p (e x , e v , R d ) f d (e x , e v ) (10)
ω c = γ R (R e ) + R T e ω p (e x , e v , R d ), (11)
where γ p : R 3 × R 3 → R 3 and γ R : SO(3) → R 3 are static state feedback stabilizers, c(•) is an interconnection term, to be selected, and ω p (e x , e v , R d ) := (R T p Ṙp ) ∨ is the angular velocity reference computed by the attitude planner (see [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF], [START_REF]Dynamic attitude planning for trajectory tracking in underactuated UAVs[END_REF] for details).

Proposition 1: Consider dynamics ( 4)-( 6) with the stabilization errors in [START_REF]Dynamic attitude planning for trajectory tracking in underactuated UAVs[END_REF] and controller ( 10)- [START_REF] Sontag | Comments on integral variants of ISS[END_REF]. If f d (e x , e v ) = 0 ∀(e x , e v ) ∈ R 6 , then the closed-loop dynamics reads:

ėx = e v (12) m ėv = γ p (e x , e v ) + ∆R(R e , e x , e v , R d ) f d (e x , e v ) (13) Ṙe = R e γR (R e ) (14)
where

∆R(R e , e x , e v , R d ) := c(R e )R p (e x , e v , R d )R e R T p (e x , e v , R d ) -I 3 . ( 15
)
Figure 1. Cascade structure of the error dynamics (equations ( 12)-( 14)).

Proof. We start by inspecting the position dynamics in (4), ( 5) and refer to the position and velocity errors defined in [START_REF]Dynamic attitude planning for trajectory tracking in underactuated UAVs[END_REF]. By taking their time derivative along ( 4)-( 5), we get

ėx = e v (16) m ėv = m ( v -vd ) = -mge 3 + T c Re 3 . ( 17
)
By noting that T c e 3 = c(R e )R T p f d , equation [START_REF] Mayhew | Quaternion-based hybrid control for robust global attitude tracking[END_REF] becomes

m ėv = -mge 3 + c(R e )RR T p f d . (18)
The velocity error dynamics can be further rewritten as

m ėv = -mge 3 + c(R e )R p R T p RR T p f d = -mge 3 + c(R e )R p R e R T p f d . (19)
Finally, adding and subtracting f d from equation ( 9), we get

m ėv = γ p (e x , e v ) + c(R e )R p R e R T p -I 3 f d , (20)
corresponding to ∆R(R e , e x , e v , R d ) in [START_REF] Khalil | Nonlinear systems, 3rd Edition[END_REF]. For the attitude dynamics, using the definition of R e in (7), we get:

Ṙe = ṘT p R + R p Ṙ = -ωp R T p R + R T p R ω (21) = R T p R( ωc -(R p R) T ωp R T p R) = R e ( ωc -R T e ωp R e ). (22
) Finally, using controller [START_REF] Sontag | Comments on integral variants of ISS[END_REF] and using property (R e ω p ) ∧ = R T e ωp R e , the closed-loop ( 12)-( 14) is obtained.

Remark 1: The position error ( 12)-( 13) is affected by the attitude error through the term ∆R f d , which weighs the mismatch between the desired force f d in [START_REF] Chaillet | Combining iISS and ISS with respect to small inputs: the strong iISS property[END_REF] and the control force resolved in F I , i.e., R f c . As shown in [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF], the interconnection function c(•) in [START_REF] Khalil | Nonlinear systems, 3rd Edition[END_REF] helps improving the transient performance by shaping the mismatch term.

B. Cascade stability analysis and stabilizers

Control law ( 10)-( 11) induces a pseudo cascade structure of the closed-loop dynamics (Figure 1) while leaving freedom in the selection of the attitude and position stabilizers γ R (•) and γ p (•, •), as well as in the selection of the connection term ∆R, which depends upon the scaling function c(•). Clearly, their selection must guarantee asymptotic stabilization of the equilibrium (R e , e x , e v ) = (I 3 , 0, 0). The following definition [START_REF] Arcak | A unifying integral ISS framework for stability of nonlinear cascades[END_REF] will be used to characterize the admissible iISS gains when the speed of convergence of the perturbing subsystem is exponential.

Definition 1:

A function η(•) is said to be of class-H I if it is of class-K and satisfies 1 0 η(s) s ds < ∞.
The following result is a slightly modified version of [START_REF] Arcak | A unifying integral ISS framework for stability of nonlinear cascades[END_REF]Corollary 2], adapted to our system evolving on R 6 × SO(3).

Lemma 1: Let x ∈ R 6 and R ∈ SO(3) and consider the cascade ẋ = f (x, R), Ṙ = Q(R) where f : R 6 ×SO(3) → R 6 and Q : SO(3) → T SO(3) are smooth vector fields with f (0, I 3 ) = 0 and Q(I 3 ) = 0. Suppose that the equilibrium R = I 3 is locally exponentially stable (LES) for Ṙ = Q(R), namely there exist c R , λ R ∈ R >0 and ∈ (0, 1) such that R(0

) SO(3) ≤ implies R(t) SO(3) ≤ c R R(0) SO(3) exp(-λ R t), t ≥ 0. ( 23
)
If the dynamics ẋ = f (x, R) is iISS with respect to input R with a class-H I gain, namely, if there exist a class-K ∞ function α(•), a class-KL function β (•), and a class-H I gain η(•) such that, for all t ≥ 0,

α( x(t) ) ≤ β ( x(0) ,t) + t 0 η( R(τ) SO(3) )dτ, (24)
then (x, R) = (0, I 3 ) is locally asymptotically stable with domain of attraction including the set

(x, R) : x ∈ R 6 , R ∈ SO(3), R SO(3) ≤ . Proof.
The proof hinges upon the iISS property guaranteed by inequality (24). In particular, from (23) we have that R(0

) SO(3) ≤ implies ∞ 0 η( R(τ) SO(3) )dτ ≤ ∞ 0 η c R R(0) SO(3) exp(-λ R τ) dτ = 1 λ R c R R(0) SO(3) 0 η(s) s ds (25)
where we used s

:= c R R(0) SO(3) exp(-λ R τ). Then, because η(•) is of class-H I , the function ν( s) := 1 λ R s 0 η(s) s ds (26)
is well-defined for all s ∈ R ≥0 and it is of class-K because ν(0) = 0 and η(s) s > 0 for all s > 0. Thus, from property (24), we can write:

α( x(t) ) ≤ β ( x(0) ,t) + ν(c R R(0) SO(3) ). ( 27
)
This proves that all solutions are bounded. Furthermore, since ∞ 0 η( R(τ) SO(3) )dτ is bounded, (24) implies that x(t) → 0 as t → ∞ following the arguments in [START_REF] Sontag | Comments on integral variants of ISS[END_REF]Prop. 6]. System ( 12)-( 14) can be written as in Lemma 1 by defining e p := (e x , e v ) ∈ R 6 and R := R e ∈ SO(3). Therefore, the stabilization of the cascade ( 12)-( 14) is possible if one can provide a suitable attitude stabilizer γ R (•) in ( 14) guaranteeing (local) exponential stability 1 of the equilibrium R e = I 3 , a position stabilizer γ p (•, •) and a scaling function c(•) for which the position error dynamics ( 12)-( 13) is iISS with respect to input R e in the sense of (24). In this spirit, we define the following properties.

Property 1: The attitude stabilizer R e → γ R (R e ) is continuous and such that the equilibrium point R e = I 3 is locally exponentially stable for [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF] with domain of attraction containing R ∈ SO(3) : R SO(3) ≤ for some ∈ (0, 1). Property 2: The position stabilizer (e x , e v ) → γ p (e x , e v ) is continuously differentiable and the scaling function R e → c(R e ) is continuous and they are such that 1 Due to topological obstructions of SO(3), only a local exponential result can be achieved with continuous static state feedback [START_REF] Koditschek | The Application of Total Energy as a Lyapunov Function for Mechanical Control Systems[END_REF]. Nonetheless, it can be shown that the basin of attraction can be extended up to the set R ∈ SO(3) : R SO(3) < 1 which does not contain only rotations about any axis of πrad away from I 3 . Moreover, the results can be extended to global ones by using hybrid stabilizers [START_REF] Mayhew | Quaternion-based hybrid control for robust global attitude tracking[END_REF].

1) there exist class-K ∞ and KL functions α(•) and β (•) and a class-H I gain η(•) such that the (e x , e v )-subsystem in equations ( 12), ( 13), ( 15) is iISS from R e as in (24); 2) there exist saturation levels M i ∈ R >0 , i ∈ {1, 2, 3}, satisfying

M 3 < mg, 3 ∑ i=1 M i ≤ T M -mg (28)
and such that the components of γ p are bounded as

|γ p i (e x , e v )| ≤ M i ∀(e x , e v , i) ∈ R 6 × {1, 2, 3}. Remark 2:
The second item in Property 2 and the continuous differentiability of γ p (•, •) are required in order to obtain a planner reference R p (e x , e v , R d ) as defined in equations ( 8)-( 9) that is well-defined and [6, Section IV-C] for additional details). Based on the above properties and Lemma 1, the main result of the paper is then given by the following theorem.

C 1 ∀(e x , e v , R d ) ∈ R 3 × R 3 × SO(3) (see
Theorem 1: Consider the closed-loop system described by ( 4)-( 6) controlled by ( 8)- [START_REF] Sontag | Comments on integral variants of ISS[END_REF]. If γ R (•), γ p (•, •) and c(•) are selected according to Properties 1 and 2, then for any desired constant trajectory (x d , 0, R d ) ∈ R 3 × R 3 × SO(3) satisfying Assumption 1, the control law ( 10)-( 11) solves Problem 1, in particular the point (e x , e v , R e ) = (0, 0, I 3 ) is asymptotically stable with domain of attraction containing

R 3 × R 3 × R ∈ SO(3) : R SO(3) ≤ .
Proof. The proof follows the same steps as the proof of Lemma 1. Indeed, since the planner reference is C 1 and well-defined (see Remark 2), the control input ω c in ( 11) is well-defined as well, and Proposition 1 ensures that the steps in Lemma 1 can be completed. Moreover, thanks to (28) and Assumption 1, the control force is bounded by T M .

IV. SAMPLE STABILIZERS DESIGN

In this section we present a selection of the basic components of the control law, namely γ p (•, •), γ R (•), c(•) and show that they satisfy Properties 1 and 2, thereby guaranteeing the applicability of Theorem 1. Specifically, we adopt a quasi time-optimal strategy for position stabilization, adapted from [START_REF] Forni | A family of global stabilizers for quasi-optimal control of planar linear saturated systems[END_REF], which consists of a locally Lipschitz state feedback that behaves linearly in a neighborhood of the origin, but coincides with the time-optimal bang-bang feedback in the presence of large velocity errors. The QTO stabilizer is a representative candidate of aggressive solutions with excellent transient performance, but limited robustness properties: while we are not able to prove that the QTO enjoys a small-signal ISS property, which makes it unfit for existing strategies, we show that it enjoys integral ISS as in (24).

A. Attitude stabilization

The next proposition gives an example of an attitude stabilizer for which Property 1 is satisfied.

Proposition 2: Given γ R (R e ) := - k R √ 1+tr(R e ) skew(R e ) ∨ ,
where k R ∈ R >0 is a scalar gain, Property 1 is satisfied for any ∈ (0, 1). In particular, the (closed-loop) trajectories of ( 14) converge exponentially to R e = I 3 for all initial conditions starting in set S R := R ∈ SO(3) : R SO(3) < 1 .

Proof. The proof follows from Lyapunov arguments, by considering the candidate Lyapunov function [13, equation 12.4]), and continuously differentiable in S R . Indeed, by taking its time derivative along the dynamics [START_REF] Forni | A family of global stabilizers for quasi-optimal control of planar linear saturated systems[END_REF] and by exploiting the trace operator property tr(A ŷ) = tr(skew(A) ŷ) = -2(skew(A) ∨ ) T y, A ∈ R 3×3 , y ∈ R 3 , we get:

V R (R e ) := k R (2 -1 + tr(R e )) which is quadratic in the error, namely, k R R e 2 SO(3) ≤ V R (R e ) ≤ 2k R R e 2 SO(3) (see
VR (R e ) = - k R 2 √ 1+tr(R e ) tr( Ṙe )) = (14) - k R 2 √ 1+tr(R e ) tr(R e γR (R e )) = - k R 2 √ 1+tr(R e ) tr (skew(R e ) γR (R e )) = k R (skew(R e ) ∨ ) T √ 1+tr(R e ) γ R (R e ) = -γ R (R e ) 2 , (29)
which is continuous and negative definite in 3) ∀R e ∈ S R , equation ( 29) can be written as

S R (note that γ R (R e ) is well-defined in S R since tr(R) = -1 ⇔ R SO(3) = 1). Furthermore, since it can be shown that γ R (R e ) = k R R e SO(
VR (R e ) = -k 2 R R e 2 SO(3) ≤ -k R 2 V R (R e ) ∀R e ∈ S R . ( 30
)
If we take R e (0)

∈ S V R := {R ∈ SO(3) : V R (R) < 2k R }, tr(R e (0) = -1 since tr(R e ) = -1 ⇔ V R (R e ) = 2k R and V R (R e
) is forced to decrease by continuity by virtue of (30) and γ R (R e ), in turn, is well-defined. Accordingly, for all initial conditions in S V R , the solutions of system (14

) satisfy V R (R e (t)) = V R (R e (0)) exp -k R 2 t ∀t ≥ 0.
It is worth remarking that as subsets, S R and S V R are equal, and therefore S R is positively invariant. Finally, since V R is quadratic, it can be easily shown that R e (t) SO(3) ≤ √ 2 R e (0) SO(3) exp -k R 4 t ∀t ≥ 0, for all initial conditions starting in S R . Finally, Property 1 is satisfied since for any ∈ (0, 1) we have that R ∈ SO(3) : R SO(3) ≤ ⊆ S R .

Remark 3: The stabilizer used in Proposition 2 is only one among several possible alternatives. By exploiting the angle axis parametrization

S 2 × (-π, π) (n, θ ) → R(n, θ ) ∈ SO(3), it is readily seen that γ R (R e (n, θ )) = k R 1-cos(θ ) 2
≤ k R , which reveals that the proposed solution is bounded for any R e ∈ S R and that lim θ →±π γ R = k R . For large attitude errors, this is more aggressive than the standard selection

γ 0 R (R e ) := k R skew(R e ) ∨ which has a vanishing magnitude for θ → ±π, since γ 0 R (R e (n, θ )) = k R | sin(θ )|.
For further details about this comparison see [START_REF] Lee | Exponential stability of an attitude tracking control system on so (3) for large-angle rotational maneuvers[END_REF].

B. Position stabilization

A smooth version of the quasi time-optimal stabilizer for saturated double integrators, originally proposed in [START_REF] Forni | A family of global stabilizers for quasi-optimal control of planar linear saturated systems[END_REF], and also used in [START_REF] Andreetto | Quasi time-optimal hybrid trajectory tracking of an n-dimensional saturated double integrator[END_REF], is employed as the position stabilizer in a decoupled form γ p (e x , e v ) = (γ p 1 (e x 1 , e v 1 ), γ p 2 (e x 2 , e v 2 ), γ p 3 (e x 3 , e v 3 )):

γ p i (e x i , e v i ) := -σ M i k x i e x i + e v i µ |e v i | 2M i , k v i k x i , (31)
for i ∈ {1, 2, 3}, where k x i , k v i ∈ R >0 are scalar PD-like gains assigning the small-signal linear behavior, σ M i (•) denotes a continuously differentiable version of the saturation function satisfying σ M i (s)s > 0 if s = 0, whose components σ M i are bounded away from zero for large values of s and are globally bounded by

|σ M i (s)| ≤ M i with M i ∈ R >0 arbitrarily selected. Finally function µ(a, b) := n √ a n + b n , a, b ∈ R ≤0
, for some n ≥ 1 is a smooth approximation of the maximum of two nonnegative scalars. The next proposition establishes iISS from the attitude error input of the arising position closed-loop ( 12), ( 13), (31).

Proposition 3: Consider the selection c(R e ) := 1 and the QTO position stabilizer (31). Then, Property 2 is satisfied. Proof. The proof is based on [12, Lemma 1] establishing iISS if there exists a C 1 , positive definite and radially unbounded function V (e x , e v ) and a class-H I function η(•) such that V (e x , e v ) ≤ η( R e SO(3) ). To this end, consider the following Lyapunov candidate for the unperturbed system, proposed in [START_REF] Andreetto | Quasi time-optimal hybrid trajectory tracking of an n-dimensional saturated double integrator[END_REF]:

V p (e x , e v ) := 1 2 m 3 ∑ i=1 k x i e 2 v i + 3 ∑ i=1 u i 0 σ M i (s)ds, (32)
which is positive definite and radially unbounded from the stated sector properties of σ M i . Applying [START_REF] Khalil | Nonlinear systems, 3rd Edition[END_REF]Lem. 4.3] and using the structure of the first term of V p , there exist class-K ∞ functions α m (•), α M (•) and a scalar c v ∈ R >0 such that

c v e v 2 + α m ( e x ) ≤ V p (e x , e v ) ≤ α M ( e p ), (33)
where we denoted e p := (e x , e v ). By straightforward derivations we first obtain

∂V p (e p ) ∂ e x i = -k x i σ M i (u i ) (34) ∂V p (e p ) ∂ e v i = mk x i e v i -k x i σ M i (u i ) µ i + ∂ µ i ∂ e v i e v i ( 35
)
where we used the shortcut notation µ i := µ

|e v i | 2M i , k v i k x i
. Using

µ i ≥ 0 and ∂ µ i ∂ e v i
e v i ≥ 0, the time derivative of V p along the unperturbed dynamics ( 12), ( 13), (31) (namely ∆R f d = 0) satisfies

Vp (e x , e v ) = - 1 m 3 ∑ i=1 k x i µ i + ∂ µ i ∂ e v i e v i σ 2 M i (u i ) ≤ 0. ( 36
)
The function µ i is the n-norm of the vector

|e v i | 2M i k v i k x i T
, for which one has the bound

µ i ≤ |e v i | 2M i + k v i k x i
and, by exploiting homogeneity, the same bound can be shown to hold also for

∂ µ i ∂ e v i e v i .
As a consequence, using the bounds and gradients in (33)-(35), for any r > 0, we may establish the upper bound:

e p ≥ r ⇒ ∇ e p V p (e p ) V p (e p ) ≤ c 1 + c 2 e v c v e v 2 + α m (r) ≤ c 3 , (37)
for some scalars c 1 , c 2 ∈ R >0 , and c 3 ∈ R >0 possibly growing unbounded as r approaches zero. Inspired by (37), let us now define a C 1 , positive definite and radially unbounded function

V (e x ,
where we may use ∆R f d as a shortcut for the right hand side interconnection term in [START_REF] Lee | Exponential stability of an attitude tracking control system on so (3) for large-angle rotational maneuvers[END_REF]. Consider now the chain of inequalities ∆R

f d ≤ ∆R f d ≤ ∆R F f d . Then, since ∆R(R e ) F = R T p R e R p -I 3 F = 6 -2tr(R e ) ≤ 2 √
2 R e SO(3) (by exploiting the definition of the Frobenius norm, the trace operator properties and the inequality -1 ≤ tr(A) ≤ 3 for any A ∈ SO(3) [7, eqn. (86)]) and since

f d = γ p (e x , e v ) + mge 3 ≤ ∑ 3 i=1 M i + mg ≤ T M (by selecting ∑ 3
i=1 M i ≤ T Mmg as per (28)), one gets:

∆R f d ≤ 2 √ 2T M R e SO(3) . (40)
Finally, combining bounds (39), (40), we get

V (e p ) ≤ 2 √ 2max (c 3 , c 4 ) T M R e SO(3) := η( R e SO( 3
) ), where function η(•) is clearly of class-H I , and we can apply [12, Lemma 1] to conclude that the (e x , e v )-subsystem is iISS with respect to input R e and that the QTO stabilizer (31) satisfies Property 2.

Remark 4:

The selection c(R e ) := 1 corresponds to the choice adopted by most hierarchical strategies [START_REF] Roza | A class of position controllers for underactuated VTOL vehicles[END_REF], [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using innerouter loop control paradigms[END_REF]. Different solutions can be envisaged to improve performance, e.g., c(R e ) := c-(1-e T 3 R e e 3 )) c

, c > 2, which guarantees reduced position overshoot for large initial attitude errors (see Remark 1 and [6, Section VI-C]). With this choice, it can be shown that the proof of Proposition 3 still works.

V. SIMULATION RESULTS

We present a simulation example to assess the performance of the proposed control law (10)- [START_REF] Sontag | Comments on integral variants of ISS[END_REF] when the attitude and position stabilizers are selected as in Propositions 2 and 3. We also compare the performance obtained with a nested saturations-based stabilizer, which is one of the most common solutions in the recent literature [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using innerouter loop control paradigms[END_REF], [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF], [START_REF] Roza | A class of position controllers for underactuated VTOL vehicles[END_REF] to globally stabilize the position dynamics of UAVs. To better highlight the differences, we refer to ideal conditions in which the UAV is described by the control model ( 4)-( 6). The considered UAV has mass m = 1kg and can deliver a maximum thrust of T M = 40N. The gain of the attitude stabilizer γ R (•) used in Proposition 2 is selected as k R = 40. For the QTO stabilizer (31), the gains can be freely chosen to have a desired behavior in unsaturated conditions and are set to

k v i = 9, k x i = k 2 v i

4

(critical damping). The saturation levels are set to M i = 9 for i ∈ {1, 2, 3} so as to ensure the validity of (28). The nested saturations stabilizer used for comparison is the one employed in [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using innerouter loop control paradigms[END_REF]:

γ p (e x , e v ) := -λ 2 σ 1 k 2 λ 2 e v + λ 1 σ 1 k 1 λ 1 e x (41)
where the gains are set to k 1 = 0.06, k 2 = 9 and the saturation levels to λ 1 = 5, λ 2 = 9. While we have used the same outer saturation level λ 2 ≡ M i and damping gain k 2 ≡ k v i of the QTO stabilizer, we had to follow the guidelines of [5, Prop. 1] in selecting k 1 and λ 1 , to guarantee the small-signal ISS property. Different selections of the control parameters could be considered but the constraints among them [5, eq. ( 26)] ultimately pose an intrinsic limitation to the achievable performance possibly caused by the strong ISS robustness property. As already mentioned, this is a common issue shared by nested saturations-based stabilizers [START_REF] Marchand | Global stabilization of multiple integrators with bounded controls[END_REF].

The initial state of the UAV is x(0) = (4, 4, 4)m, v(0) = (10, 0, 0)m/s and R(0) = I 3 + sin(179.9 • ) ê1 + (1cos(179.9 • ) ê2 1 ), which corresponds to an upside-down configuration with a significantly misplaced position with respect to the desired hovering position x d (0) = (0, 0, 1)m, R d = I 3 . The attitude tracking performance is illustrated in Figure 2 (top) for both controllers. Since the same control law is used for attitude stabilization, the same attitude transient is achieved. On the contrary, the position errors are characterized by a quite different behavior. The aggressive nature of the QTO stabilizer results in an impressively faster response (see Figure 2 (bottom)). This desirable aggressive behavior is confirmed by inspecting Figure 3, where the thrust and the magnitude of the commanded angular velocity are depicted.

VI. CONCLUSIONS

The stabilization problem for vectored-thrust UAVs has been addressed using an inner-outer loop paradigm and integral ISS tools to analyze the arising pseudo cascade. This allows one to employ stabilizers with superior performance with respect to existing strategies while guaranteeing the same basin of attraction. Our stabilization results can be extended to trajectory tracking following the approach in [START_REF] Invernizzi | Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF], which should allow tracking in cluttered environments via suitable planners.

b p 3 ×R d e 1 b p 3 ×R d e 1 × b p 3 b p 3 ×R d e 1 b p 3 ×R d e 1 b p 3 ,

 33133313 R p (e x , e v , R d ) := (8) where b p 3 := f d (e x ,e v ) f d (e x ,e v ) and f d (e x , e v ) := γ p (e x , e v ) + mge 3

Figure 2 .

 2 Figure 2. Stabilization errors: attitude R e SO(3) (top) and position e x .

  [1] M.-D. Hua, T.Hamel, P. Morin, and C. Samson, "Introduction to feedback control of underactuated VTOL vehicles: A review of basic control design ideas and principles," IEEE Control Systems, vol. 33, no. 1, pp. 61-75, 2013.

Figure 3 .

 3 Figure 3. Control inputs: thrust T c (top) and angular velocity ω c .

  e v ) := ln(V p (e x , e v )) if V p (e x , e v ) > e 1 e V p (e x , e v ) if V p (e x , e v ) ≤ e (38) where e denotes Euler's number. Using (37) and the fact that the gradient of V p is bounded in any compact set e p ≤ r, we obtain ∇V (e p ) ≤ c 4 for some c 4 ∈ R >0 . Combining this bound with the complementary one in (37), inspecting the perturbed dynamics (12), (13), (31) and exploiting the unperturbed bound (36), we may write V (e p ) ≤ max (c 3 , c 4 ) ∆R(R e , e x , e v , R d ) f d (e x , e v ) ,

Research supported in part by ANR via grant HANDY, number ANR-18-CE40-0010.


Davide Invernizzi

Marco Lovera

Luca Zaccarian

Dynamic attitude planning for trajectory tracking in thrust-vectoring UAVs

Keywords: UAVs, thrust vectoring, trajectory tracking

This paper addresses the trajectory tracking control problem for underactuated UAVs, with specific focus on vehicles with thrustvectoring capabilities. According to the different actuation mechanisms, the most common UAV platforms can achieve only a partial decoupling of attitude and position tasks. Since position tracking is of utmost importance for applications involving aerial vehicles, we propose a control scheme in which position tracking is the primary objective. To this end, this work exploits the concept of attitude planner, a dynamical system through which the desired attitude reference is processed to guarantee the satisfaction of the primary objective: the attitude tracking task is considered as a secondary objective which can be realized as long as the desired trajectory satisfies specific trackability conditions.

I. INTRODUCTION

The increasing demand of complex and challenging applications involving Vertical-Take-Off and Landing Unmanned Aerial Vehicles (VTOL UAVs) has led to the design of novel configurations to overcome the maneuverability limitation of standard platforms. In particular, multirotor UAVs with coplanar propellers cannot fulfill at the same time attitude and position tasks. We refer to these platforms as vectored-thrust UAVs because their propulsive system can deliver a control force only along a fixed direction within the airframe.

By designing an actuation mechanism that can change the thrust direction within the airframe, more complex maneuvering can be performed. Such platforms, which we call thrust-vectoring UAVs, can be realized by mounting the propellers in a non-coplanar fashion [START_REF] Crowther | Kinematic analysis and control design for a nonplanar multirotor vehicle[END_REF], [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally bounded input force[END_REF] or by employing servo-actuators to dynamically adjust their orientation [START_REF] Hua | Control of VTOL vehicles with thrust-tilting augmentation[END_REF], [START_REF] Invernizzi | Full pose tracking for a tilt-arm quadrotor UAV[END_REF], [START_REF] Ryll | A novel overactuated quadrotor unmanned aerial vehicle: Modeling, control, and experimental validation[END_REF]. When the thrust vectoring capabilities are limited, the tracking control problem is challenging since only a partial decoupling is possible between the attitude and position objectives. In particular, our work starts by showing that it is impossible to track an arbitrary full pose trajectory (independent attitude and position) also for thrust-vectoring configurations, when the vectoring capability is limited. In this case a viable approach is to prioritize control objectives, for instance by subduing the realization of attitude tracking to the achievement of position tracking [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally bounded input force[END_REF], [START_REF] Hua | Control of VTOL vehicles with thrust-tilting augmentation[END_REF]. While much literature is devoted to the trajectory tracking problem for the vectored-thrust configuration [START_REF] Hamel | Dynamic modelling and configuration stabilization for an X4-Flyer[END_REF], [START_REF]Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF], [START_REF] Lee | Geometric tracking control of a quadrotor UAV on SE(3)[END_REF], [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using inner-outer loop control paradigms[END_REF], [START_REF] Roza | A class of position controllers for underactuated VTOL vehicles[END_REF], [START_REF] Zhao | Quaternion-based trajectory tracking control of VTOL-uavs using command filtered backstepping[END_REF], few works address the thrust-vectoring one [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally bounded input force[END_REF], [START_REF] Hua | Control of VTOL vehicles with thrust-tilting augmentation[END_REF], [START_REF] Invernizzi | Trajectory tracking control of thrustvectoring UAVs[END_REF].

The control design that we propose is developed in two stages. The first step is the design of a control force and torque that ensure robust tracking of any desired trajectory that possesses certain smoothness and boundedness properties. In the second step, we design a dynamic attitude planner to account for the system underactuation and to enforce prioritization of the position tracking objective over the attitude one in the spirit of ideas proposed in [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally bounded input force[END_REF], [START_REF] Hua | Control of VTOL vehicles with thrust-tilting augmentation[END_REF]. In developing our control law we take as well inspiration from [START_REF] Roza | A class of position controllers for underactuated VTOL vehicles[END_REF], in which the emphasis was on defining classes of position and attitude controllers that stabilize the UAV at a given position. In the same way, our control design parametrizes a family of control solutions to the tracking problem for underactuated UAVs provided that they possess some basic properties. Moreover, our contribution extends the results of [START_REF] Roza | A class of position controllers for underactuated VTOL vehicles[END_REF], which only deals with stabilization (constant reference) and vectored-thrust platforms and the results of D. Invernizzi, M. Lovera are with Department of Aerospace Science and Technology, Politecnico di Milano, Via La Masa 34, 20156, Milano, Italy.

L. Zaccarian is with LAAS-CNRS, Université de Toulouse, CNRS, Toulouse, France and with Dipartimento di Ingegneria Industriale, University of Trento, Italy zaccarian@laas.fr [START_REF]Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF], which only deals with a static attitude planning and vectoredthrust platforms. This extension is non-trivially accomplished here by characterizing a non-autonomous closed-loop by way of setvalued maps and it allows considering a large class of attitude control laws that do not necessarily make the attitude error subsystem autonomous nor exponentially stable [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally bounded input force[END_REF], [START_REF] Hua | Control of VTOL vehicles with thrust-tilting augmentation[END_REF]. As in our previous work [START_REF]Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF], the control design leaves freedom also in the choice of the position feedback stabilizer, provided that it guarantees certain properties: a basic requirement is the position dynamics be Input-to-State Stable (ISS) with respect to small inputs. Such property is ensured by employing, e.g., nested saturations-based stabilizers [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using inner-outer loop control paradigms[END_REF]. Then, we show that our control design provides robust tracking with semi-global properties. To achieve this result, we study the stability of the closed-loop interconnection between the attitude and position error systems within the framework of differential inclusions. This approach simplifies the analysis of the corresponding cascaded non-autonomous system: the proof technique relies on casting the control problem as a stability problem for a compact attractor with dynamics satisfying regularity conditions, which ensures robustness of the stability property against a very large class of (sufficiently small) perturbations [START_REF] Goebel | Hybrid Dynamical Systems[END_REF]Chapter 7]. Our proof is based on reduction theorems [START_REF] Maggiore | Reduction theorems for hybrid dynamical systems[END_REF] which have been exploited with a different control strategy in [START_REF] Michieletto | Nonlinear control of multi-rotor aerial vehicles based on the zero-moment direction[END_REF] to address set-point tracking (stabilization). Finally, we show that existing attitude planning strategies developed for vectored-thrust [START_REF] Lee | Geometric tracking control of a quadrotor UAV on SE(3)[END_REF] and thrust-vectoring platforms [START_REF] Invernizzi | Trajectory tracking control of thrustvectoring UAVs[END_REF] fit within the proposed robust control design framework.

Notation. The canonical basis in R n is denoted as ei := (0, . . . , 1, . . . , 0) for i ∈ {1, ..., n} and the identity matrix in R n×n is denoted as

In := [e1 • • • ei • • • en]. Given x = (x1, . . . , xn) ∈ R n , x := x 2 1 + . . . + x 2
n is the Euclidean norm while for a matrix A ∈ R n×n , A F := tr(A T A) is the Frobenius norm and skew(A)

:= A-A T 2 is the skew-symmetric part of A. A ∈ R n×n >0 (R n×n <0
) denotes a positive (negative) definite matrix. The minimum and maximum eigenvalues of A ∈ R n×n are λm(A) and λM (A), respectively, and skew(A)

:= A-A T 2 is the skew-symmetric part of A. Set SO(3) := {R ∈ R 3×3 : R T R = I3, det(R) = 1}
denotes the third-order Special Orthogonal group. The normalized distance of R ∈ SO(3) with respect to I3, induced by the Frobenius norm, is denoted as 3 and × represents the cross product in R 3 . The inverse of the hat map is the vee map, denoted as (•) ∨ : so(3) → R 3 . The set of piecewise-continuous and bounded functions is denoted as L∞. We use standard comparison functions: a function α : R ≥0 → R ≥0 is of class K if it is zero at zero, strictly increasing, and continuous. It is of class

R SO(3) := 1 8 R - I3 F = 1 4 tr(I3 -R) ∈ [0, 1]. Given ω ∈ R 3 , the hat map • : R 3 → so(3) := {Ω ∈ R 3×3 : Ω = -Ω T } is such that ωy = ω × y, ∀y ∈ R
K∞ if it is also unbounded. A function β : R ≥0 × R ≥0 → R ≥0 is of class KL if it is of class K in
the first argument and nonincreasing and converging to zero as its second argument tends to +∞.

II. DYNAMICAL MODELS FOR THRUST-VECTORING UAVS

This section recalls the simplified dynamical model of multirotor UAVs that is usually employed for control design purposes. In order to obtain a control design independent from the specific actuation mechanism, we will assume the control wrench wc := (fc, τc), where fc ∈ R3 and τc ∈ R3 are the control force and torque, respectively, to be the control variable. Then, we will resort to approximate models of the wrench map (from physical inputs to delivered wrench), according to most common actuation mechanisms. The equations of motion for control design are: When the range of the map from the set of physical inputs U to the delivered wrench, namely U u → wc(u), spans R 6 , the system is fully actuated. In the following, the control torque τc is assumed to span R 3 , whilst the actuation mechanism allows to deliver the control force fc only on a compact subset of R 3 . Within this category of UAVs, a further classification distinguishes between vectored-thrust [START_REF]Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF] and thrust-vectoring configurations [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally bounded input force[END_REF]. As per the first class, the control force can be delivered only in the common direction of the thrust axes of the propellers:

ẋ = v (1) Ṙ = Rω (2) m v = -mge3 + Rfc (3) J ω = -ωJω + τc, (4)
fc 3 (t) > 0, fc 1 (t) = fc 2 (t) = 0 ∀t ≥ 0. ( 5
)
This configuration is adopted in most multirotor UAVs thanks to the inherent mechanical simplicity combined with good performance in many flight conditions. The second class comprises UAVs with thrust-vectoring capabilities such as multirotor UAVs both with fixed-tilted [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally bounded input force[END_REF] or dynamically tiltable propellers [START_REF] Ryll | A novel overactuated quadrotor unmanned aerial vehicle: Modeling, control, and experimental validation[END_REF]. While the dynamically tiltable configuration makes the system fully actuated if the tilt angle of the servo-actuators is not limited, this is not true for the fixed-tilted configuration, for which the maximum inclination at which the rotors are mounted is limited by the efficiency loss: the power consumption in hover is large and increases with the inclination of the propellers. For the latter configuration, the control force fc has to satisfy:

0 < cos(θM ) ≤ fc(t) T e3 fc(t) = cos(θc(t)) ∀t ≥ 0, ( 6
)
where θM is the maximum tilt angle and • represents the Euclidean norm . Finally, for all the cases above, we assume that the control force is bounded by a positive scalar TM ∈ R>0, namely fc(t) ≤ TM ∀t ≥ 0. This constraint is required to account for actuators saturations, in particular the limited spinning velocity of propellers.

III. CONTROL PROBLEM STATEMENT

This work deals with the tracking control problem for the system described by ( 1)-( 4). The objective is to track a reference

t → (R d (t), ω d (t), x d (t), v d (t)) ∈ T SO(3) × R 3 × R 3 , where T SO(3) := SO(3) × R 3 is the trivial tangent bundle of SO(3).
When considering the actuation limitations shown in the previous section, the tracking of an arbitrary reference trajectory in T SO(3) × R 3 × R 3 is not possible due to the platform underactuation. However, platforms with thrust-vectoring capabilities can achieve a certain degree of decoupling between attitude and position tracking when the corresponding trajectories are compatible in some sense. To suitably represent and exploit these degrees of freedom, let us first compute the steady state (short, ss) control wrench, which is obtained by inverting the system dynamics:

f ss c (t) := mR T d (t) ( vd (t) + ge3) (7) τ ss c (t) := J ωd (t) + ωd (t)Jω d (t). (8)
For arbitrary position and attitude trajectories, it is likely that f ss c will not be compatible with the actuation limitations ( 5) and ( 6).

On the other hand, because position tracking is mandatory in aerial applications, equation [START_REF] Hamel | Dynamic modelling and configuration stabilization for an X4-Flyer[END_REF] suggests that the desired attitude can be properly modified to be compliant with the actuation constraints. Indeed, according to equation [START_REF] Hamel | Dynamic modelling and configuration stabilization for an X4-Flyer[END_REF], the control force is obtained by rotating the vector m ( vd + ge3) by R T d . The rationale behind the proposed control is to prioritize position over orientation tracking, as already suggested by [START_REF] Hua | Control of VTOL vehicles with thrust-tilting augmentation[END_REF] and [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally bounded input force[END_REF]. We will propose a strategy to track the desired attitude as long as it allows to follow the desired position. Whenever this condition cannot be granted, only the closest feasible attitude will be tracked. We will exploit the fact that when the control force is delivered along the vertical direction of the body frame (b3) (vectored-thrust case), the constraint (6) for the thrust-vectoring case is satisfied as well.

Before going on, the following assumptions (not too restrictive for standard applications) are required to hold for the desired trajectory.

Assumption 1: Smoothness and boundedness of the desired trajectory. The desired trajectory t

→ (R d (t), ω d (t), x d (t), v d (t)) ∈ T SO(3) × R 3 × R 3 satisfies Ṙd (t) = R d (t)ω d (t) and ẋd (t) = v d (t)
, for all t ≥ 0, and is continuously differentiable up to the second order. Furthermore, the desired linear acceleration vd (•) and angular velocity ω d (•) are bounded.

IV. ROBUST STABILIZATION OF THE ERROR DYNAMICS

In this section we show the design of a control law that ensures position tracking for any bounded and sufficiently smooth trajectory with no restriction on the initial position error and some restriction on the initial attitude error. Specifically, we will show that the system dynamics can be represented by a feedback interconnection in which an attitude planner provides the reference to the attitude subsystem which, in turn, affects the position error dynamics (see Figure 1).

The attitude planner design plays a central role in our control strategy. This subsystem takes the desired attitude reference and the position errors as inputs and is in charge of computing a reference attitude and angular velocity that are feasible, in the sense that the actuation constraints are satisfied. Specifically, the output of the planner (which is the actual reference to the attitude subsystem) is denoted as (Rp, ωp) ∈ T SO(3).

Property 1: The attitude planner output (Rp, ωp) ∈ T SO( 3) is feasible, in the following sense

Ṙp(t) = Rp(t)ωp(t), ∀t ≥ 0, ωp(•) ∈ C 1 ∩ L∞. (9)

A. Control force and position error dynamics

Since the goal of this work is to obtain a control law that prioritizes position over orientation tracking, we start by inspecting the position error dynamics. By focusing on the dynamics in (1)-( 4), the natural choice of the tracking errors for position and velocity is:

ex := x -x d , ev := v -v d . (10)
Accordingly, the position and velocity errors are expressed in the inertial frame. The control objective is to stabilize the origin (ex, ev) = (0, 0) ∈ R 3 × R 3 . Consider the equations ( 1), ( 3) and the definition of the errors in [START_REF] Invernizzi | Dynamic attitude planning for trajectory tracking in underactuated UAVs[END_REF]. Then,

ėx = ev (11) m ėv = m ( v -vd ) = -m ( vd + ge3) + Rfc. ( 12
)
The system would be fully actuated if one could arbitrarily assign fc, and it would be possible to enforce Rfc = f d , with

f d (ex, ev, vd ) := β (ex, ev) + m ( vd + ge3) , (13)
where β : R 3 × R 3 → R 3 is a static state feedback stabilizer. The corresponding closed-loop dynamics would be described by [START_REF]Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF] and m ėv = β (ex, ev). Thus, the origin would be globally asymptotically stable. In this work we will adopt the same strategy used in [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using inner-outer loop control paradigms[END_REF], corresponding to a nested saturations-based static nonlinear feedback:

β (ex, ev) := λ2sat k2 λ2 ev + λ1sat k1 λ1 ex , (14)
where sat(•) is a smooth decentralized unit saturation function satisfying the properties mentioned in [Sec. II-A] [17], k1, k2 and λ1, λ2 are suitably chosen gains and saturation levels to be selected according to the guidelines of [Prop. 1] [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using inner-outer loop control paradigms[END_REF]. While we develop our theory for this specific stabilizer, we emphasize that similar generalizations to those reported in [START_REF]Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF], [START_REF] Roza | A class of position controllers for underactuated VTOL vehicles[END_REF] are possible, thus allowing for stabilizing laws inducing improved transients.

To account for the fact that it is not possible to obtain Rfc = f d when the control force cannot span R 3 , it is convenient to rewrite the velocity error dynamics as:

m ėv = -m ( vd + ge3) + RR T p Rpfc (15) = -m ( vd + ge3) + ReRpfc. ( 16
)
where Rp is the reference attitude given by the attitude planner and

Re := RR T p ∈ SO(3) (17)
is the so-called left attitude error [START_REF] Bullo | Tracking for fully actuated mechanical systems: a geometric framework[END_REF]. Introducing the corresponding angular velocity error,

eω := ω -ωp, (18)
a natural choice for the control force fc is

fc := R T p Φ (Re, eω, t) f d , (19)
where Φ(Re, eω, t) : T SO(3)×R ≥0 → R 3×3 is a design function, selected later, ensuring that Φ(Re, eω, t) → I3 for (Re, eω) → (I3, 0) and ∀t ≥ 0. Indeed, by adding and subtracting f d in equation ( 16), the closed-loop velocity error dynamics reads:

m ėv = β (ex, ev) + ∆R (Re, eω, t) f d (ex, ev, vd ) (20)
where ∆R(Re, eω, t) := ReΦ (Re, eω, t) -I3.

Written in this form, the position error is clearly affected by the attitude error through the term ∆Rf d , which is the mismatch between the desired force f d and the actual control force resolved in the inertial frame, i.e., Rfc.

Remark 1: The control force is obtained by scaling the desired force f d by a suitable term Φ (Re, eω, t) dependent on the attitude error, and then by applying a rotation R T p to the resulting vector. The idea behind the proposed control law is twofold. When the attitude error is large, the scaling term Φ (Re, eω, t) can be used to reduce the control force and, correspondingly, the overshoot in the position tracking is limited and the transient behavior can be improved. The second point is related to the fact that the desired attitude may be such that the constraints on fc expressed by equations ( 5)-( 6) are still not respected with the control force defined in equation [START_REF] Ryll | A novel overactuated quadrotor unmanned aerial vehicle: Modeling, control, and experimental validation[END_REF]. However, by prioritizing position over orientation tracking one can consider the planner rotation matrix Rp as an additional degree of freedom to ensure that the scaled control force Φ (Re, eω, t) f d is eventually rotated to be compliant with the actuation constraint.

The rationale behind the proposed control law is that if the attitude error dynamics can be made asymptotically stable, then, for t → ∞, Re(t) → I3, Φ (Re(t), eω(t), t) → I3 and the control force delivered in the inertial frame

R(t)fc(t) = Re(t)Φ(t)f d (t) → f d (t)
, which is the force in the inertial frame required to track the desired position trajectory. Indeed, in this case, the mismatch between the desired and actual control force converges to zero as well, namely ∆R(t)f d (t) → 0. Firstly, the choice of Φ must be such that the following property holds true for the attitude mismatch ∆R.

Property 2: (Vanishing perturbations). Consider ∆R defined in [START_REF] Zhao | Quaternion-based trajectory tracking control of VTOL-uavs using command filtered backstepping[END_REF]. Given Va(Re, eω) := eω 2 + Re 2 SO(3) , there exists a bounded class-K function γ(•), satisfying ∀ (Re, eω, t) ∈ T SO(3) × R ≥0 : ∆R(Re, eω, t) ≤ γ (Va (Re, eω)) .

(22) We will now show a possible selection of Φ such that the above property is verified. A wide range of performance-oriented alternative choices are possible as long as they satisfy Property 2, which is needed to ensure that, at convergence, the magnitude of the delivered control force fc converges to the nominal force f ss c [START_REF] Hamel | Dynamic modelling and configuration stabilization for an X4-Flyer[END_REF]. It is straightforward to employ a scaling transformation, dependent on the attitude error alone, as follows:

Φ(Re, eω, t) := c(Re, t)I3, (23)
where c : SO(3) × R ≥0 → R is a properly selected function. The next proposition, whose proof is given in [START_REF] Invernizzi | Dynamic attitude planning for trajectory tracking in underactuated UAVs[END_REF]Appendix], gives an example of such scaling function, which is naturally written in terms of the angle θe between the desired direction bp 3 := Rpe3 and the vertical body axis b3 := Re3:

θe(Re, Rp) := arccos(b T p 3 b3) = arccos(e T 3 R T p ReRpe3). (24) Proposition 1: Given c(Re, t) := -(1-cos(θe(Re,Rp(t))))
, where θe is defined in (24), then, for Φ(Re, eω, t) := c(Re, t)I3 and > 2, Property 2 is satisfied.

B. Control torque and attitude error dynamics

The attitude controller has to ensure the convergence of the attitude tracking errors according to the fully actuated rotational dynamics in equations ( 2), [START_REF] Forni | A family of global stabilizers for quasi-optimal control of planar linear saturated systems[END_REF]. By using matrix multiplication as the group operation, the attitude error Re = RR T p , is employed as the attitude error measure in SO(3), which was already introduced in equation ( 17). The error kinematics are derived from (2) and ( 9) as follows:

Ṙe = ṘR T p + R ṘT p = ReRpêωR T p . (25)
Consider the system defined in (2), ( 4) and the control law

τc := -R T p eR -Kωeω + J ωp + ωpJω, (26)
where Kω ∈ R 3×3 is positive definite and eR := skew(KRRe) ∨ ∈ R 3 is the left-trivialized derivative of the modified trace function ΨK R (Re) := 1 2 tr(KR(I -Re)). The map (•) ∨ : so(3) → R 3 used in the definition of eR represents the inverse of the • map defined in the notation section. It is assumed that matrix KR ∈ R 3×3 is symmetric and that it satisfies:

tr(KR)I3 -KR ∈ R 3×3 >0 . (27)
Using (25) and the control torque (26) with the rotational equations of motion (2), (4), the dynamics of the errors ( 17) and ( 18) reads:

Ṙe = ReRpêωR T p (28) J ėω = -R T p eR -Kωeω -êωJeω -êωJωp. (29)
The control torque (26), first proposed by [START_REF] Bullo | Tracking for fully actuated mechanical systems: a geometric framework[END_REF], has a simpler expression than the one based on the right group error considered in [START_REF] Lee | Geometric tracking control of a quadrotor UAV on SE(3)[END_REF] and no cancellation of non-harmful nonlinearities occurs. The equilibria for the attitude subsystem are the points where the differential of ΨK R and the angular velocity error vanish, i.e., eR = 0, eω = 0. The set of equilibria contains the desired equilibrium (Re, eω) = (I3, 0) and additional undesired configurations corresponding to the other critical points of ΨK R . Due to this intrinsic structure of SO(3), no time-invariant continuous control law can globally stabilize (I3, 0). Nonetheless, by defining the scalar (positive from (27)):

R := λm (tr(KR)I3 -KR) > 0, (30)
it is well known that in the sublevel set

SR := {Re ∈ SO(3) : ΨK R (Re) < R} , (31)
the point Re = I3 is the unique critical point and minimum of ΨK R . Next, the total error energy function

VR (Re, eω) := 1 2 e T ω Jeω + ΨK R (Re) , (32)
will be used in the stability analysis. It can be shown that this function is a Lyapunov candidate for the attitude error dynamics, i.e., it is positive definite about (Re, eω) = (I3, 0), continuously differentiable and radially unbounded in the direction of eω → ∞. The following theorem, the proof of which is reported in Section VII to avoid breaking the flows of the exposition, establishes desirable properties of the attitude stabilizer. Theorem 1: Consider the system described by ( 28)-( 29) and a reference attitude t → (Rp(t), ωp(t)) ∈ T SO(3) satisfying Property 1. For any symmetric matrix KR ∈ R 3×3 satisfying (27) and any matrix Kω ∈ R 3×3 >0 , the equilibrium point (Re, eω) = (I3, 0) is uniformly asymptotically stable (UAS) with basin of attraction containing the set

Sa := {(Re, eω) ∈ T SO(3) : VR (Re, eω) < R} (33)
where R and VR are defined in (30) and (32). Remark 2: By choosing KR = kRI3, kr ∈ R>0, which satisfies (30) with R = 2kR, the sublevel set of ΨK R in (31) contains all the rotations with an angle strictly less than 180 • , namely almost all the configurations in SO [START_REF] Crowther | Kinematic analysis and control design for a nonplanar multirotor vehicle[END_REF]. By increasing the gain kR, also the set of initial conditions eω(t0) included in (33) can be arbitrarily enlarged achieving a semi-global stability property. Robust global asymptotic stability of the desired equilibrium (I3, 0) can be obtained by using a hybrid controller on SO(3) [START_REF] Mayhew | Quaternion-based hybrid control for robust global attitude tracking[END_REF].

Remark 3: As an alternative choice to (26), the control law τc := ωJω+J ωp -ωpω -R T p (eR + KωRpeω) substituted in (2), (4) provides the following closed-loop:

Ṙe = Re êω (34) ėω = -eR -Kω ẽω, (35)
where ẽω := Rp(ω-ωp). In this way the closed loop is autonomous and independent of ωp. Even though this choice would make the stability analysis of the control scheme simpler and would place less restrictions on the position stabilizer (see Remark 6), it would also make the design less general. For instance, the feedback τc cancels non-harmful nonlinearities making it less desirable from the point of view of actuation effort.

C. Complete dynamics

This section presents the main results of the stability analysis for the complete system. Our proof is based on a compact representation of the closed loop wherein the solutions of the time-varying dynamics is embedded into a time-invariant differential inclusion, in ways that are similar to the strategy in [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using inner-outer loop control paradigms[END_REF], even though the approach adopted here does not require the (somewhat stringent) assumption ωp be bounded. By introducing xa := (Re, eω) and xp := (ex, ev), the solutions to the attitude error dynamics (28)-( 29) and the position error dynamics ( 11), ( 20) can be embedded within the solution funnel generated by the following constrained differential inclusion1 :

(A) ẋa ∈ Fa(xa), xa ∈ T SO(3) (36) (P ) ẋp ∈ Fp (xp, xa) , xp ∈ R 3 × R 3 , (37)
where co(•) denotes the closed convex hull, B3 is the closed unit ball and ωM ∈ R>0 is a constant existing thanks to Assumption 1. Moreover, function γ comes from (22), and scalar fM := √ 3λ2 + f ss M is a bound on f d arising from substituting [START_REF] Hamel | Dynamic modelling and configuration stabilization for an X4-Flyer[END_REF] and ( 40) into [START_REF] Lee | Geometric tracking control of a quadrotor UAV on SE(3)[END_REF].

Based on representation (36), (37), asymptotic tracking for the complete dynamics can be proven, under Assumption 1, as stated by the following Theorem.

Theorem 2: Consider the closed loop system described by ( 11), ( 20) and ( 28), (29) controlled by ( 19), (26) and the planner output given by (Rp, ωp)

= (R d , ω d ) ∈ T SO(3), where the desired trajectory t → (R d (t), ω d (t), x d (t), v d (t)) satisfies Assumption 1.
Then, if Φ(•, •, •) is selected according to Property 2, for any symmetric matrix KR satisfying (27), any Kω ∈ R 3×3 >0 , the point (Re, eω, ex, ev) = (I3, 0, 0, 0) is robustly uniformly asymptotically stable with basin of attraction containing the set Sa ×R 3 ×R 3 , where Sa is defined in (33).

Proof: The cascaded interconnection (36), (37) comprises the upper subsystem (36), whose stability properties (with domain of attraction Sa) is established in Theorem 1, and the lower subsystem (37), which is stabilized by the nested saturation feedback [START_REF] Maggiore | Reduction theorems for hybrid dynamical systems[END_REF] proposed in [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using inner-outer loop control paradigms[END_REF]. (Local) stability of the cascade follows from standard reduction theorems for differential inclusions (see, e.g., [START_REF] Maggiore | Reduction theorems for hybrid dynamical systems[END_REF]) whereas attractivity from Sa × R 3 × R 3 can be established using the small signal ISS properties of stabilizer [START_REF] Maggiore | Reduction theorems for hybrid dynamical systems[END_REF] following the same steps as in [17, Proof of Prop. 4]. Finally, stability and attractivity of the point (Re, eω, ex, ev) = (I3, 0, 0, 0) for the closed loop implies robust uniform asymptotic stability from Sa × R 3 × R 3 , by concatenating the implications in [6, Thm 7.12] and [6, Thm 7.21].

V. ATTITUDE PLANNING FOR CONSTRAINTS COMPLIANCE

Theorem 2 establishes robust asymptotic stability of the origin for the error dynamics, regardless of the reference orientation (Rp, ωp) ∈ T SO(3). Nonetheless, it gives no guarantees about the fact that the force fc requested by the control scheme satisfies the constraints characterized in Section II. In this section we propose to select (Rp, ωp) according to a dynamic attitude planner, as represented in Figure 1, which is in charge of properly changing the desired attitude to prioritize position over orientation tracking, while respecting the input constraints at hand. With respect to the approach of [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally bounded input force[END_REF], in which the reference attitude is obtained as the solution of an optimization problem, our design allows to naturally compute a differentiable reference that satisfies Property 1 so that τc in (26) is well defined. Note also that the control design of [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally bounded input force[END_REF], [START_REF] Hua | Control of VTOL vehicles with thrust-tilting augmentation[END_REF] limits the choice of attitude stabilizers to those making the attitude error subsystem autonomous and the desired equilibrium exponentially stable.

Before proceeding with the planner design, we recall the following additional assumptions on the desired trajectory which are needed to ensure a well-defined planner reference at all times.

Assumption 2: The desired trajectory t → (R d (t), ω d (t), x d (t), v d (t)) ∈ T SO(3) × R 3 × R 3 satisfies Assumption 1 and 1) the desired acceleration vd (•) is such that the nominal force ( 7) is bounded by some strictly positive constants f ss m , f ss m , i.e.,

0 < f ss m ≤ f ss c (t) ≤ f ss M < TM ∀t ≥ 0, (40)
and that inf t≥0 (m|g + vd 3 |) > 0; 2) the desired trajectory x d (•) belongs at least to C 4 ; 3) the desired jerk vd (•) belongs to L∞

A. Attitude planner dynamics

There is a natural way to express the reference attitude Rp by noticing that the actuation constraints for all the different configurations are satisfied if the control force is delivered along the positive direction of the b3 axis. Following standard strategies for underactuated UAVs [START_REF] Lee | Geometric tracking control of a quadrotor UAV on SE(3)[END_REF], we introduce a smooth matrix function Rc(f d , t) : R 3 \ {0} × R ≥0 → SO(3) defined as:

Rc(f d , t) := bc 1 (f d , t) bc 2 (f d , t) bc 3 (f d ) (41)
where bc 3 := f d / f d and bc 1 (f d , t) and bc 2 (f d , t) are any two orthogonal unit vectors such that Rc defines a rotation matrix. A possible selection is:

bc 1 := bc 2 (f d , t) × bc 3 (f d ), bc 2 := bc 3 (f d ) × b d (t) bc 3 (f d ) × b d (t) . (42)
The vector b d defines the desired heading direction of the UAV:

b d (t) := cos(ψ d (t)) sin(ψ d (t)) 0 T , (43)
where t → ψ d (t) ∈ R is the desired yaw angle (which may be extracted from a given R d (t)). We note that Rc(f d , t) is well defined as long as f d = 0 0 0 T . Based on Rc, we select the reference attitude Rp, output of the attitude planner, as:

Rp := Rc(f d , t)Rr, (44)
where Rr ∈ SO(3) is an additional state of the dynamic attitude planner. Being an element of SO(3), the differential equation for the relative attitude Rr can be written as:

Ṙr = ωrRr, (45)
where ωr ∈ R 3 is the relative angular velocity, with coordinates in the frame Fc, that will be used as an auxiliary input to control the evolution of Rr. For instance, ωr can be used to track at best the desired attitude R d by exploiting a Lyapunov-based design, once a suitable potential function of the desired attitude R d is provided. Then, the time evolution of Rr can be properly modified in order to satisfy the actuation constraints (see Section VI). Finally, the attitude planner has to provide a corresponding velocity reference ωp, satisfying Ṙp = Rp ωp, which is computed by first introducing the angular velocity

ωc(Rc, f d ) := (R T c Ṙc) ∨ ( 46
)
of the frame FC := {bc 1 , bc 2 , bc 3 }, and then using ( 44) and ( 45) to obtain Ṙp = ṘcRr + Rc Ṙr = RpR T r (ωc + ωr)Rr = Rp ωp. In particular, the above relation provides ωp = R T r (ωc + ωr)Rr which, by exploiting the vee (•) ∨ map gives:

ωp = R T r (ωc(Rc, f d ) + ωr) (47)
.

As the main goal of the attitude planner is to track at best the desired attitude, it is more convenient to rewrite the dynamics (45) in terms of the planner attitude error, i.e.,

R p e := RpR T d . ( 48
)
Using then Ṙd = R d ωd from Assumption 1 and (47), the overall dynamics of the attitude planner becomes:

Ṙp e = R p e R d (ωp -ωd )R T d (49) ωp = (R T c R p e R d ) T (ωc(Rc, f d ) + ωr). ( 50
)
Proposition 2: By selecting the gain λ2 < inf t≥0 |m( vd 3 (t) + g)| in ( 14), the rotation matrix Rc by ( 41)-( 43), is well defined when the desired trajectory satisfies Assumption 2.

Proof: The proposition can be proven by inspecting the following inequality:

f d (t) ≥ |f d 3 (t)| ≥ m| vd 3 (t) + g| -|β3(ex, ev)| (51) ≥ inf t≥0 |m( vd 3 (t) + g)| -β3 M > 0, (52)
which holds thanks to the equivalence β3 M = λ2 for the definition of ( 14) coming from [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using inner-outer loop control paradigms[END_REF] and the assumption that λ2 < inf t≥0 |m( vd 3 (t) + g)|. Then, bc 1 and bc 2 are unit vectors, orthogonal to each other and with bc 3 and well defined ∀t ≥ 0. Indeed, b d (t) belongs to the horizontal plane by definition (43) and the third component of bc 3 never vanishes from (51). Hence, the cross product bc 3 × b d does not vanish either and so bc 3 × b d = 0 ∀t ≥ 0 in equation ( 42).

The following lemma is a useful link between the output (Rp, ωp) of the attitude planner, and the reference motion satisfying the properties of Assumption 1. Its proof is given in [START_REF] Invernizzi | Dynamic attitude planning for trajectory tracking in underactuated UAVs[END_REF]Appendix].

Lemma 1: (Feasibility conditions of the planner output). If the relative angular velocity is bounded and continuously differentiable, i.e., ωr ∈ (L∞ ∩ C 1 ), and the desired angular velocity ω d satisfies Assumption 1, then the reference attitude motion (Rp, ωp) ∈ T SO(3), obtained according to equations ( 44), (47), is feasible, in the sense that it satisfies Property 1.

Remark 4: The decomposition of the reference attitude in equation (44) allows to effectively account for the actuation constraints of the different configurations. The angular velocity ωr of the relative rotation matrix is an additional degree of freedom that can be exploited when the thrust-vectoring capability is not locked, i.e., the control force can be produced in a region around the vertical axis. This additional input can be exploited to track the desired attitude at best while taking into account the constraints. Indeed, it is always possible to select an initial condition Rr(0) such that the actuation constraints are verified. Then, the evolution of the relative attitude can be properly controlled by modifying the angular velocity input ωr to satisfy the constraints.

Remark 5: The requirement x d (•) ∈ C 4 in Assumption 1 allows us to properly define the time derivative of ωp, which is required to apply the attitude control law, as shown in the next section. Indeed, by direct computation from (47):

ωp = -R T r ωrωc + R T r ( ωc + ωr) , (53)
where ωc = (R T c Rc -ω2 c ) ∨ .

B. Special selections

In this section we illustrate the relevance of the proposed control scheme for addressing the input limitations of thrust-vectoring UAVs (case (b) ). A similar analysis about vectored-thrust configurations (case (a)) can be found in [START_REF] Invernizzi | Dynamic attitude planning for trajectory tracking in underactuated UAVs[END_REF], [START_REF]Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF].

Within the present design, we make use of the relative rotation matrix Rr introduced in (44) to perform attitude maneuvers that are compatible with position tracking. By considering that Rc defined in (41) satisfies Rce3 = f d f d , the control force can be written as:

fc = c(Re, t)R T r R T c f d = c(Re, t)R T r f d e3 = c(Re, t) f d R T r e3.
(54) By substituting (54) into the constraint (6), we get:

cos(θc) = c(Re, t) f d R T r e3 c(Re, t) f d = e T 3 Rre3 = e T 3 br 3 , (55)
which shows that in order to satisfy [START_REF] Goebel | Hybrid Dynamical Systems[END_REF], it is sufficient to guarantee that e T 3 br 3 ≥ cos(θM ). In this section we will show how the solution proposed in [START_REF] Invernizzi | Trajectory tracking control of thrustvectoring UAVs[END_REF] to compute the relative angular velocity ωr of equation (45) to account for the conic region constraint [START_REF] Ryll | A novel overactuated quadrotor unmanned aerial vehicle: Modeling, control, and experimental validation[END_REF] and exploit the thrust-vectoring capabilities of tiltrotor configurations can be applied within the present framework. In particular, we will verify that the planner output obtained by selecting

ωr = br 3 × Proj G ω d r × br 3 + b T r 3 ω d r br 3 , (56)
where Proj G : R3 → R 3 is a geometric projection operator 3 and 1+ε , δ , where ε ∈ (0, 1) is a user-defined parameter. In particular, thanks to the projection operator, e T 3 br 3 (t) ≥ cos(θM ) ∀t ≥ 0 and, by virtue of (55), also:

ω d r = Rrω d -ωc -RrR T d e p R (57)
cos (θc(t)) ≥ cos (θM ) ∀t ≥ 0. ( 58
)
As a consequence, if the relative attitude is initialized such that br 3 is inside a cone defined by θM around bc 3 , it will never leave it. Lemma 2: The relative angular velocity computed according to equation ( 56) is such that the planner output (Rp, ωp) ∈ T SO(3) obtained with (44), (47), satisfies Lemma 1.

Proof: By exploiting a smooth Projection operator, ωr defined in (56) is C 1 and its time derivative can be computed as in (53). Then, we can write the following bound on ωr:

ωr ≤ Proj G (ω d r × br 3 ) + ω d r (59)
Notice that the projection operator simply removes the radial component of ḃd r 3 , when it is working, hence Proj G (ω d r × br 3 ) ≤ ω d r because br 3 is a unit vector. Finally, we can conclude the boundedness of ω d r from the following inequality:

ω d r ≤ ω d + ωc + e p R , (60)
in which ω d is bounded according to Assumption 1, ωc is bounded as shown in [START_REF] Invernizzi | Dynamic attitude planning for trajectory tracking in underactuated UAVs[END_REF]Appendix], and e p R is bounded, being the lefttrivialized derivative of a function defined on a compact manifold. , there exists t ∈ R ≥0 , such that

Definition 1: A trajectory t → (R d (t), ω d (t), x d (t), v d (t)) ∈ T SO(3) × R 3 × R 3 satisfying
e T 3 f ss c (t) f ss c (t) ≥ cos(θ b ) ∀t ≥ t, (61)
where f ss c (t) is defined in equation ( 7). According to Definition 1, the attitude motion is compatible with the conic region constraint [START_REF] Goebel | Hybrid Dynamical Systems[END_REF] if the angle between the third desired axis b d 3 and the nominal force R d f ss c (t) is within the cone in which the Proj G operator is not active. The next corollary combines the results of Theorem 2 and Lemma 2.

Corollary 1: Consider the closed loop system described by ( 11), ( 20), (28), (29), (49) controlled by [START_REF] Ryll | A novel overactuated quadrotor unmanned aerial vehicle: Modeling, control, and experimental validation[END_REF], (26), where the planner output (Rp, ωp) is given by ( 44) and (50), with Rc, ωc and ωr selected as in (41)-( 43), ( 46) and ( 56)-(57), respectively. Assume that the desired trajectory t → (R d (t), ω d (t), x d (t), v d (t)) is trackable according to Definition 1, that Φ(•, •, •) is selected according to (23) and satisfies Property 2, and that the gains k1, k2 and saturation levels λ1, λ2 in ( 14) are selected according to [START_REF] Naldi | Robust global trajectory tracking for underactuated VTOL aerial vehicles using inner-outer loop control paradigms[END_REF]Prop. 1] and Proposition 2. Then, for any symmetric matrix KR satisfying (27), Kω ∈ R 3×3 >0 , ε ∈ (0, 1) and k d ∈ R>0, the control force [START_REF] Ryll | A novel overactuated quadrotor unmanned aerial vehicle: Modeling, control, and experimental validation[END_REF] satisfies the actuation constraint [START_REF] Goebel | Hybrid Dynamical Systems[END_REF] and the point (Re, eω, ex, ev, R p e ) = (I3, 0, 0, 0, I3) is robustly asymptotically stable with basin of attraction containing the set Sa × R 3 × R 3 × Sap, where Sa is given by (33) and Sap :

= R p e ∈ SO(3) : e T 3 R T d (0)(R p e ) T f d (0) f d (0) ≥ cos(θM )
, with f d defined in [START_REF] Lee | Geometric tracking control of a quadrotor UAV on SE(3)[END_REF].

Proof: The proof is based on Theorem 2, Lemma 2 and [12, Thm 3]. Lemma 2 guarantees that the output planner reference (Rp, ωp) ∈ T SO(3) computed according to (44) and (56) satisfies Property 1. Therefore, the assumptions of Theorem 2 are satisfied and the point (Re, eω, ex, ev) = (I3, 0, 0, 0) is UAS. Furthermore, because R p e evolves on SO(3) according to (49), (50) with ωr given by (56) and because the projection operator guarantees that equation (58) holds ∀t ≥ 0, the first part of the theorem is demonstrated. Then, following the proof [START_REF] Invernizzi | Trajectory tracking control of thrustvectoring UAVs[END_REF]Thm 3], combining the trackability condition in Definition 1 and the properties of the projection operator together with the expression of ω d r in (57), the theorem can be proved.

Remark 6: (Classes of position and attitude stabilizers). Any position stabilizer β(•, •) satisfying the conditions of [Property 2] [START_REF]Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF] would work in Corollary 1. At the same time, different selections of the control torque (see [START_REF] Bullo | Tracking for fully actuated mechanical systems: a geometric framework[END_REF] for some alternatives beyond the one in Remark 3) could be considered as long as they guarantee UAS of (Re, eω) = (I3, 0) for any planner reference t → (Rp(t), ωp(t)). Following Remark 3, our selection in equation ( 26) is quite general in terms of characteristics of the attitude error subsystem that one has to consider when addressing the stability analysis of the complete control scheme. While this makes our stability analysis quite general, it comes at the expense of more stringent requirements on the position stabilizer (second item in [Property 2] [START_REF]Geometric trajectory tracking with attitude planner for vectored-thrust VTOL UAVs[END_REF]) to guarantee that the planner velocity is bounded (Property 1).

VI. NUMERICAL RESULTS

The control law designed in the previous sections will be applied to the case of a hexarotor UAV in which each propeller is tilted by an angle α with respect to the local x-axis (see Figure 2). The mapping from u := (ω1, . . . , ω6) to the control wrench wc can be compactly written as wc = M (α)ū [10, eqn. 74], where ωi ∈ [0, ωM ] is the angular rate of the i-th propeller, ū := (ω 2 1 , . . . , ω 2 6 ) and M (α) ∈ R 6×6 is a matrix dependent upon physical quantities of the UAV and the tilt angle α. Note that nonfeasible, i.e., negative, angular rates for the individual propellers can be obtained when inverting M (α) for a given control wrench. This can be understood by considering that the force delivered according to M (α)ū spans only a given region in the space around b3, the size of which is closely related to the value of α. The model of case (b) in Section II is reasonable to approximate the actuation constraint of the considered hexarotor.

The simulation model is a multi-body system with seven bodies (a central body and six propellers groups) written in the Modelica modeling language, which accounts for the actuators dynamics as well as aerodynamic effects. A conservative source of drag has been considered in the model to verify robustness, which is one of the main claims of Theorem 2. Furthermore, the control law has been implemented in a 500Hz sample-data scenario [START_REF] Franchi | Full-pose tracking control for aerial robotic systems with laterally bounded input force[END_REF] with a time delay of two samples, together with typical measurement noise experienced in indoor flight conditions with a motion capture system. The mass and inertial matrix of the UAV are m = 1 kg and J = diag (0.008, 0.008, 0.016) kg • m 2 , respectively. By considering α = 20 • , the conic region [START_REF] Goebel | Hybrid Dynamical Systems[END_REF] with θM = 10 • was found to be a reasonable approximation of the wrench map M (α) to ensure its invertibility in a broad operational range. The controller gains are KR = diag(0.6, 0.6, 1.4), Kω = 0.2I3, = 2.1, k d = 2, λ2 = 9, λ1 = 1, k1 = 0.06 k2 = 9, ε = 0.05. The desired position trajectory is a circle x d (t) = (cos(Ω d (t)t), cos(Ω d (t)t), 0)m, see Figure 2, where Ω d (t) is made by two constant intervals connected by a smooth ramp (top of Figure 3). The desired attitude is R d (t) = I3, i.e., the UAV has to fly with level attitude, a requirement that is never compatible with vectored-thrust configurations. The initial conditions correspond to hover at x(0) = (1, 0, 0) m. By inspecting Figure 3 (bottom), the considered trajectory is not always feasible in the sense of Definition 1. Indeed, the inclination angle of the steady state control force [START_REF] Hamel | Dynamic modelling and configuration stabilization for an X4-Flyer[END_REF] with respect to the vertical body axis, i.e., θn := arccos(e T 3 f ss c (t)/ f ss c (t) ), is greater than the maximum tilt angle θM in the interval t ∈ [17, 26]s. It is worth mentioning that θn represents also the inclination that a vectored-thrust UAV would have to reach to guarantee position tracking. Figure 4 shows that, after a transient phase, the position tracking errors remain bounded: aerodynamic drag prevents their actual convergence to zero. The oscillations shown at steady state are induced by the combination of aerodynamic drag (which works in the direction of -v) and the periodic nature of the circular motion. The attitude tracking performance is shown in Figure 5, where the inclination angle of the body axis b3 with respect to the inertial axis e3, i.e., θv := arccos(b T 3 e3), is plotted (top) together with the yaw angle ψ (bottom). Since the desired attitude is R d (t) = I3, the corresponding desired values are θ d v (t) = 0 • and ψ d (t) = 0 • . Additional plots showing the control variables can be found in the extended report [START_REF] Invernizzi | Dynamic attitude planning for trajectory tracking in underactuated UAVs[END_REF]. When the trajectory is feasible, the UAV is capable of flying at almost level attitude (θv(t) ≈ 0.5 • ). On the contrary, during the initial transient and the acceleration phase, when the desired attitude is not compatible with the thrust-vectoring constraint (6) and position tracking, the attitude tracking objective is only partially achieved and the projection operator is working to modify the attitude reference so that position tracking is guaranteed. Note from Figure 2 and Figure 5 that the vehicle is inclined of an angle θv = 11.62 • at t ≈ 21s whereas the nominal angle is θn( t) ≈ 21 • > θM . Therefore, the proposed solution tries to stay as close as possible to the desired attitude (θ d v = 0) even if the trajectory is not trackable. Furthermore, note that the primary objective is not affected during this phase: the position tracking performance is not deteriorated. It is worth mentioning that the tracking errors converge to zero (see Figure 4 and Figure 5, dashed lines) when the control model (CM) is employed in the simulation, as expected from Corollary 1. Finally, the conic region constraint ( 6) is satisfied at all times, i.e., θc(t) ≤ θM , as shown in Figure 3 (bottom, solid line).

VII. PROOF OF THEOREM 1

To prove Theorem 1, we write the attitude error dynamics as a constrained differential inclusion and then we apply an invariance principle to prove the asymptotic stability of the desired attractor. Equations (28), (29) describe the evolution of the attitude error dynamics, the solutions of which are also solutions of the differential inclusion ẋa ∈ Fa, where Fa, defined in (38), is clearly an outer semicontinuous (its graph is closed) and locally bounded set-valued map, and for each xa ∈ T SO(3), Fa (xa) is nonempty and convex. To prove Theorem 1 it will be convenient to intersect Fa (xa) with the tangent cone to T SO(3), denoted as T T SO(3) (xa) at xa. Indeed the solutions to (36) can only flow along the directions in the tangent cone, which simplifies the Lyapunov analysis (see [START_REF] Seuret | A nonsmooth hybrid invariance principle applied to robust eventtriggered design[END_REF] for details). To this end, the following lemma is useful: To prove asymptotic stability, we employ the following invariance principle, which is a corollary of [START_REF] Seuret | A nonsmooth hybrid invariance principle applied to robust eventtriggered design[END_REF]Thm 1].

Proposition 3: Consider ẋ ∈ F (x), x ∈ C, where C is a closed set, F is a set-valued mapping outer semicontinuous and locally bounded relative to C, F (x) is nonempty and convex ∀x ∈ C. Given a compact set A, if there exists a continuously differentiable function V , positive definite and radially unbounded around A relative to C and such that

V (x) = max f ∈F (x)∩T C (x) ∇V (x), f ≤ 0 ∀x ∈ C \ A, ( 63
)
then A is globally stable, namely Lyapunov and Lagrange stable. Furthermore, 1) if there exists an open neighborhood U ⊃ A from which no complete solution γ exists, satisfying V (γ(t))) = V (γ(0)) = 0, then A is asymptotically stable. 2) Any such neighborhood U of the form

U := {x ∈ C : V (x) < } , > 0, (64)
is contained in the basin of attraction of A. Proof of Theorem 1. Since all solutions to (28), (29) are also solution to (36), we prove Theorem 1 by applying Proposition 3 to (36) with A = {I3, 0} =: {x * a }, with the Lyapunov candidate VR in (32) which is clearly positive definite and radially unbounded. In particular, for V = VR, we first prove (63) and then we prove item 2) of Proposition 3 with U = Sa, according to (31) and (64). Verification of (63). From Lemma 3, we may evaluate VR in ( 63

γ(t) will exit the set G for some small t > t0, since eR = 0 ∀xa ∈ U \ {x * a }. As a consequence, VR(γ(t)) is forced to decrease again. We can conclude that there is no complete solution that keeps VR(γ(t)) constant and different from zero. Therefore, from Proposition 3 the proof is complete.

VIII. CONCLUSIONS

In this paper we proposed a priority-oriented robust control paradigm for a general class of underactuated UAVs, in which position tracking is the primary objective. To this end, we exploited a dynamic attitude planner in charge of providing a modified attitude reference, which guarantees that the control force, required to track the desired position, can always be delivered. A numerical simulation has been performed by considering a hexarotor with tilted propellers to verify the robustness of the control law against unmodeled effects such as unknown mass distribution, aerodynamic drag and disturbance torques. Future work will be oriented to improve the position tracking performance by considering highperformance feedback stabilizers [START_REF] Forni | A family of global stabilizers for quasi-optimal control of planar linear saturated systems[END_REF] which require a non-trivial extension of the robust design proposed here.

where x ∈ R 3

 3 and v ∈ R 3 are position and velocity of the center of mass in the inertial frame, respectively, R = [ b 1 b 2 b 3 ] ∈ SO(3) is the rotation matrix describing the attitude of the body fixed frame FB := {b1, b2, b3} with respect to the inertial frame, ω ∈ R 3 is the body angular velocity, m ∈ R>0 and J = J T ∈ R 3×3 >0 are the UAV mass and inertia matrix with respect to the principal axes.

Fig. 1 .

 1 Fig. 1. Control scheme: cascade interconnection.

1 m

 1 with a slight abuse of notation 2 and where Fa(xa), Fp(xp) are Fa(xa) := (38) Rp ∈ co (SO (3)) ωp ≤ ω M ReRpêωR T p -J -1 R T p eR + Kωeω + êωJeω + êωJωp , Fp(xp, xa) := ev β(ex, ev) + fM γ(Va(Re, eω))B3

  with e p R := k d skew(R p e ) ∨ , satisfies Property 1. The projection operator keeps the planar vector b ⊥ r 3 := [ e T 1 br 3 e T 2 br 3 ] inside a circle of radius δ := sin(θM ) by acting on the vector field ω d r in the region δ √

  Assumption 2 is compatible with the position tracking task for the system of case (b) if, given δ := sin(θM ), ε ∈ (0, 1) and θ b := arcsin δ √ 1+

Fig. 2 .

 2 Fig. 2. Trajectory followed by the hexacopter.

Fig. 3 .

 3 Fig. 3. Ω d (t) (top) -Maximum tilt angle θ M , nominal angle θn of the steady state force with respect b 3 and control force angle θc (bottom).

Fig. 4 .

 4 Fig. 4. Position tracking error ex (dashed lines: control model).

Fig. 5 .

 5 Fig. 5. Attitude tracking -Inclination (top) and yaw angle (bottom) (solid lines: multi-body model -dashed lines: control model).

Lemma 3 :

 3 Given Fa defined in (38) and the closed set T SO(3), we have the following:Fa (xa) ∩ T T SO(3) (xa) = (62) Rp ∈ SO (3) ωp ≤ ω M RpReR T p êω -J -1 R T p eR + Kωeω + êωJeω + êωJωp .Proof of Lemma 3. For a smooth manifold, the tangent cone is equivalent to the tangent space, namely T T SO(3) (xa) = Tx a T SO(3). Since T SO(3) is a trivial bundle, i.e., T SO(3) = SO(3) × R 3 , then we can write T (Re,eω ) T SO(3) TR e SO(3) × R 3 . Hence, equality (62) follows from the observation that ReRpêωR T p ∈ TR e SO(3) ←→ Rp ∈ SO(3).

  ) with F = Fa and Fa (xa) ∩ T T SO(3) (xa) given in (38), as follows (where we use e T ω (êωa) = 0 and ΨK R (Re) = -1 2 tr(KRRe) = e T R Rpeω) VR(xa) = -e T ω Kωeω ≤ 0, ∀xa ∈ T SO(3) (65) Verification of item 2) of Proposition 3. According to (64) and (31), select U = {xa ∈ T SO(3) : VR(xa) < } for some > 0. Thanks to the properties of the potential function ΨK R , the sublevel set the form {xa ∈ T SO(3) : VR(xa) < }, where ≤ R, with R defined in (30), contains only the desired equilibrium point. Indeed, VR(xa) < R implies ΨK R (Re) < R for any xa ∈ T SO(3) and only the desired equilibrium point is contained in this sublevel set. Furthermore, since this set is forward invariant from (65) and the viability condition of [6, Pag. 124]) is satisfied, maximal solutions starting in U are complete. We refer now to a solution t → γ(t) starting in U for which VR(γ(0)) = a = 0. Then, if we consider γ(0) ∈ {xa ∈ U : eω = 0}, the function VR(γ(t)) has to decrease in time by continuity, which implies VR(γ(t)) < 0. Instead, if γ(0) ∈ G := {xa ∈ U : eω = 0} \ {x * a }, then, according to the closed loop dynamics Fa(xa) G = Rp ∈ co (SO (3)) ωp ≤ ω M 0 -J -1 R T p eR ,

Following [6, §6.2.1], the regularity of the data of our constrained differential inclusions simplifies the definition of solutions of ẋ ∈ F (x), x ∈ C, corresponding to any absolutely continuous function φ such that φ(t) ∈ F (φ(t)) and φ(t) ∈ C for almost all t.

To be consistent with the formulation, the differential inclusion should be written by exploiting the vectorization, vec( Ṙe) ∈ vec ReRp êωR T p

The geometric projection operator is an extension of the smooth projection operator for systems evolving on R n originally proposed in[START_REF] Cai | A sufficiently smooth projection operator[END_REF], to the case of systems evolving on S n .