Assalé Adjé

Pierre-Loïc Garoche

Victor Magron
email: v.magron@imperial.ac.uk

Property-based Polynomial Invariant Generation using Sums-of-Squares Optimization

While abstract interpretation is not theoretically restricted to specific kinds of properties, it is, in practice, mainly developed to compute linear over-approximations of reachable sets, aka. the collecting semantics of the program. The verification of user-provided properties is not easily compatible with the usual forward fixpoint computation using numerical abstract domains. We propose here to rely on sums-of-squares programming to characterize a property-driven polynomial invariant. This invariant generation can be guided by either boundedness, or in contrary, a given zone of the state space to avoid. While the target property is not necessarily inductive with respect to the program semantics, our method identifies a stronger inductive polynomial invariant using numerical optimization. Our method applies to a wide set of programs: a main while loop composed of a disjunction (ifthen-else) of polynomial updates e.g. piecewise polynomial controllers. It has been evaluated on various programs.

Introduction

With the increased need for confidence in software, it becomes more than ever important to provide means to support the verification of specification of software. Among the various formal verification methods to support these analysis, a first line of approaches, such as deductive methods or SMT-based model checking, provide rich languages to support the expression of the specification and then try to discharge the associate proof obligation using automatic solvers. The current state of the art of these solvers is able to manipulate satisfiability problems over linear arithmetics or restricted fragments of non linear arithmetics. Another line of approaches, such as static analysis also known as abstract interpretation, restricts, a priori, the kind of properties considered during the computation: these methods typically perform interval arithmetic analysis or rely on convex polyhedra computations. In practice this second line of work seems more capable of manipulating and generating numerical invariants through the computation of inductive invariants, while the first line of approaches hardly synthesize these required invariants through satisfiability checks.

However, when it comes to more than linear properties, the state of the art is not well developed. In the early 2000s, ellipsoid analyses [START_REF] Feret | Static analysis of digital filters[END_REF], similar to restricted cases of Lyapunov functions, were designed to support the study of a family of Airbus controllers. This exciting result was used to provide the analysis of absence of runtime errors but could hardly be adapted to handle more general user provided specifications for polynomial programs.

However proving polynomial inequalities is NP-hard and boils down to show that the infimum of a given polynomial is nonnegative. Still, one can obtain lower bounds of such infima by decomposing certain nonnegative polynomials into sums-of-squares (SOS). This actually leads to solve hierarchies of semidefinite relaxations, introduced by Lasserre in [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF]. Recent advances in semidefinite programming allowed to extensively apply these relaxations to various fields, including parametric polynomial optimization, optimal control, combinatorial optimization, etc. (see e.g. [START_REF] Parrilo | Semidefinite programming relaxations for semialgebraic problems[END_REF][START_REF] Laurent | Sums of squares, moment matrices and optimization over polynomials[END_REF] for more details).

While these approaches were mentioned a decade ago in [START_REF] Cousot | Proving Program Invariance and Termination by Parametric Abstraction, Lagrangian Relaxation and Semidefinite Programming[END_REF] and mainly applied to termination analysis, they hardly made their way through the software verification community to address more general properties.

Contributions. Our contribution allows to analyze high level properties defined as a sublevel set of polynomials functions, i.e. basic semialgebraic sets. This class of properties is rather large: it ranges from boundedness properties to the definition of a bad region of the state space to avoid. While these properties, when they hold, are meant to be invariant, i.e. they hold in each reachable state, they are not necessarily inductive. Our approach rely on the computation of a stronger inductive property using SOS programming. This stronger property is proved inductive on the complete system and, by construction, implies the target property specified by the user. We develop our analysis on discrete-time piecewise polynomial systems, capturing a wide class of critical programs, as typically found in current embedded systems such as aircrafts.

Organization of the paper. The paper is organized as follows. In Section 2, we present the programs that we want to analyze and their representation as piecewise polynomial discrete-time systems. Next, we recall in Section 3 the collecting semantics that we use and introduce the polynomial optimization problem providing inductive invariants based on target polynomial properties. Section 4 contains the main contribution of the paper, namely how to compute effectively such invariants with SOS programming. Practical computation examples are provided in Section 5. Finally, we explain in Section 6 how to derive template bases from generated invariants.

Polynomial programs and piecewise polynomial discrete-time systems

In this section, we describe the programs which are considered in this paper and we explain how to analyze them through their representation as piecewise polynomial discrete-time dynamical systems. We focus on programs composed of a single loop with a possibly complicated switch-case type loop body. Moreover we suppose without loss of generality that the analyzed programs are written in Static Single Assignment (SSA) form, that is each variable is initialized at most once.

Definitions. We recall that a function f from R d to R is a polynomial if and only if there exists k ∈ N, a family

{c α | α = (α 1 , . . . , α d ) ∈ N d , |α| = α 1 + . . . + α d ≤ k} such that for all x ∈ R d , f (x) = |α|≤k c α x α1 1 . . . x α d d .
By extension a function f : R d → R d is a polynomial if and only if all its coordinate functions are polynomials. Let R[x] stands for the set of d-variate polynomials.

In this paper, we consider assignments of variables using only parallel polynomial assignments (x 1 , . . . , x d ) = T (x 1 , . . . , x d ) where (x 1 , . . . , x d ) is the vector of the program variables. Tests are either weak polynomial inequalities r(x 1 , . . . , x d ) ≤ 0 or strict polynomial inequalities r(x 1 , . . . , x d ) < 0. We assume that assignments are polynomials from R d to R d and test functions are polynomials from R d to R. In the program syntax, the notation will be either <= or <. The form of the analyzed program is described in Figure 1.

x ∈ X in ; w h i l e ( r 0 1 ( x ) 0 and . . . and r 0 n 0 ( x ) 0 ) { c a s e ( r 1 1 ( x ) 0 and . . . and r 1 n 1 ( x ) 0 ) : x = T 1 ( x ) ; c a s e . . . c a s e ( r 1 i ( x ) 0 and . . . and r 1 n i ( x ) 0 ) : x = T i ( x ) ; } Fig. 1. One-loop programs with nested conditional branches

A set C ⊆ R d is said to be a basic semialgebraic set if there exist g 1 , . . . , g m ∈ R[x] such that C = {x ∈ R d | g j (x)
0, ∀ j = 1, . . . , m}, where is used to encode either a strict or a weak inequality.

As depicted in Figure 1, an update T i : R d → R d of the i-th condition branch is executed if and only if the conjunction of tests r i j (x) 0 holds. In other words, the variable x is updated by T i (x) if the current value of x belongs to the basic semialgebraic set

X i := {x ∈ R d | ∀j = 1, . . . , n i , r i j (x) 0} . ( 1
)
Piecewise Polynomial Systems. Consequently, we interpret programs as constrained piecewise polynomial discrete-time dynamical systems (PPS for short).

The term piecewise means that there exists a partition {X i , i ∈ I} of R d such that for all i ∈ I, the dynamics of the system is represented by the following relation, for k ∈ N:

if x k ∈ X i ∩ X 0 , x k+1 = T i (x k ) . ( 2
)
We assume that I is finite and that the initial condition x 0 belongs to some compact basic semialgebraic set X in . For the program, X in is the set where the variables are supposed to be initialized in. Since the test entry for the loop condition can be nontrivial, we add the term constrained and X 0 denotes the set representing the conjunctions of tests for the loop condition. The iterates of the PPS are constrained to live in X 0 : if for some step k ∈ N, x k / ∈ X 0 then the PPS is stopped at this iterate with the terminal value x k .

We define a partition as a family of nonempty sets such that:

i∈I X i = R d , ∀ i, j ∈ I, i = j, X i ∩ X j = ∅ . ( 3
)
From Equation (3), for all k ∈ N * there exists a unique i ∈ I such that x k ∈ X i . A set X i can contain both strict and weak polynomial inequalities and characterizes the set of the n i conjunctions of tests polynomials r i j . Let r i = (r i 1 , . . . , r i ni ) stands for the vector of tests functions associated to the set X i . We suppose that the basic semialgebraic sets X in and X 0 also admits the representation given by Equation (1) and we denote by r 0 the vector of tests polynomials (r 0 1 , . . . , r 0 n0 ) and by r in the vector of test polynomials (r in 1 , . . . , r in nin ). To sum up, we give a formal definition of PPS.

Definition 1 (PPS).

A constrained polynomial piecewise discrete-time dynamical system (PPS) is the quadruple (X in , X 0 , X , L) with:

-X in ⊆ R d is the compact basic semialgebraic set of the possible initial conditions; -X 0 ⊆ R d is the basic semialgebraic set where the state variable lives; -X := {X i , i ∈ I} is a partition as defined in Equation (3); -L := {T i , i ∈ I} is the family of the polynomials from R d to R d , w.r.t. the partition X satisfying Equation (2).

From now on, we associate a PPS representation to each program of the form described at Figure 1. Since a program admits several PPS representations, we choose one of them, but this arbitrary choice does not change the results provided in this paper. In the sequel, we will often refer to the running example described in Example 1. Example 1 (Running example). The program below involves four variables and contains an infinite loop with a conditional branch in the loop body. The update of each branch is polynomial. The parameters c ij (resp. d ij ) are given parameters. During the analysis, we only keep the variables x 1 and x 2 since oldx 1 and oldx 2 are just memories.

x 1 , x 2 ∈ [a 1 , a 2 ] × [b 1 , b 2 ] ; oldx 1 = x 1 ; oldx 2 = x 2 ; w h i l e (-1 <= 0 ) { oldx 1 = x 1 ; oldx 2 = x 2 ; c a s e : oldx 1 ^2 + oldx 2 ^2 <= 1 : x 1 = c 11 * oldx 1 ^2 + c 11 * oldx 2 ^3 ; x 2 = c 21 * oldx 1 ^3 + c 22 * oldx 2 ^2 ; c a s e : -oldx 1 ^2 -oldx 2 ^2 < -1 x 1 = d 11 * oldx 1 ^3 + d 12 * oldx 2 ^2 ; x 2 = d 21 * oldx 1 ^2 + d 22 * oldx 2 ^2 ; } } The associated PPS corresponds to the quadruple (X in , X 0 , {X 1 , X 2 }, {T 1 , T 2 }),
where the set of initial conditions is:

X in = [a 1 , a 2 ] × [b 1 , b 2 ] ,
the system is not globally constrained, i.e. the set X 0 in which the variable x = (x 1 , x 2 ) lies is:

X 0 = R d ,
the partition verifying Equation (3) is:

X 1 = {x ∈ R 2 | x 2 1 + x 2 2 ≤ 1}, X 2 = {x ∈ R 2 | -x 2 1 -x 2 2 < -1}
, and the polynomials relative to the partition {X 1 , X 2 } are:

T 1 (x) = c 11 x 2 1 + c 12 x 3 2 c 21 x 3 1 + c 22 x 2 2 and T 2 (x) = d 11 x 3 1 + d 12 x 2 2 d 21 x 2 1 + d 22 x 2 2 .

Program invariants as sublevel sets

The main goal of the paper is to decide automatically if a given property holds for the analyzed program, i.e. for all its reachable states. We are interested in numerical properties and more precisely in properties on the values taken by the d-uplet of the variables of the program. Hence, in our point-of-view, a property is just the membership of some set P ⊆ R d . In particular, we study properties which are valid after an arbitrary number of loop iterates. Such properties are called loop invariants of the program. Formally, we use the PPS representation of a given program and we say that P is a loop invariant of this program if:

∀ k ∈ N, x k ∈ P ,
where x k is defined at Equation (2) as the state variable at step k ∈ N of the PPS representation of the program. Our approach addresses any property expressible as a polynomial level set property. This section defines formally these notions and develop our approach: synthesize a property-driven inductive invariant.

Collecting Semantics as postfixpoint characterization

Now, let us consider a program of the form described in Figure 1 and let us denote by S the PPS representation of this program. The set R of reachable values is the set of all possible values taken by the state variable along the running of S. We define R as follows:

R = k∈N T k | X 0 (X in ) (4)
where T | X 0 is the restriction of T on X 0 and T | X 0 is not defined outside X 0 . To prove that a set P is a loop invariant of the program is equivalent to prove that R ⊆ P . We can rewrite R inductively:

R = X in ∪ i∈I T i R ∩ X i ∩ X 0 . ( 5
)
Let us denote by ℘(R d ) the set of subsets of R d and introduce the map F :

℘(R d ) → ℘(R d ) defined by: F (C) = X in ∪ i∈I T i C ∩ X i ∩ X 0 (6)
We equip ℘(R d ) with the partial order of inclusion. The infimum is understood in this sense i.e. as the greatest lower bound with respect to this order. The smallest fixed point problem is:

inf C ∈ ℘(R d ) | C = F (C) .
It is well-known from Tarski's theorem that the solution of this problem exists, is unique and in this case, it corresponds to R. Tarski's theorem also states that R is the smallest solution of the following Problem:

inf C ∈ ℘(R d ) | F (C) ⊆ C .
Note also that the map F corresponds to a standard transfer function (or collecting semantics functional) applied to the PPS representation of a program. We refer the reader to [START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF] for a seminal presentation of this approach.

To prove that a subset P is a loop invariant, it suffices to show that P satisfies F (P ) ⊆ P . In this case, such P is called inductive invariant.

Considered properties: sublevel properties P κ,α

In this paper, we consider special properties: those that are encoded with sublevel sets of a given polynomial function.

Definition 2 (Sublevel property). Given a polynomial function κ ∈ R[x] and α ∈ R ∪ {+∞}, we define the sublevel property P κ,α as follows:

P κ,α := {x ∈ R d | κ(x) α} .
where denotes ≤ when α ∈ R and denotes < for +∞. The expression κ(x) < +∞ expresses the boundedness of κ(x) without providing a specific bound α.

Example 2 (Sublevel property examples).

Boundedness. When one wants to bound the reachable values of a system, we can try to bound the l 2 -norm of the system:

P • 2 2 ,∞ with κ(x) = x 2 2 .
The use of α = ∞ does not impose any bound on κ(x). Safe set. Similarly, it is possible to check whether a specific bound is matched. Either globally using the l 2 -norm and a specific α: P • 2 2 ,α , or bounding the reachable values of each variable: P κi,αi with κ i : x → x i and α i ∈ R. Avoiding bad regions. If the bad region can be encoded as a sublevel property k(x) ≤ 0 then its negation -k(x) ≤ 0 characterize the avoidance of that bad zone. Eg. if one wants to prove that the square norm of the program variables is always greater than 1, then we can consider the property P κ,α with κ(x) = 1 -x 2 2 and α = 0.

A sublevel property is called sublevel invariant when this property is a loop invariant. This turns out to be difficult to prove loop invariant properties while considering directly R, thus we propose to find a more tractable overapproximation of R for which such properties hold.

Approach: compute a P κ,α -driven inductive invariant P

In this subsection, we explain how to compute a d-variate polynomial p ∈ R[x] and a bound w ∈ R, such that the polynomial sublevel sets P := {x ∈ R d | p(x) ≤ 0} and P κ,w satisfy: R ⊆ P ⊆ P κ,w ⊆ P κ,α .

(7)

The first (from the left) inclusion forces P to be valid for the whole reachable values set. The second inclusion constraints all elements of P to satisfy the given sublevel property for a certain bound w. The last inclusion requires that the bound w is smaller than the desired level α. When α = ∞, any bound w ensures the sublevel property. Now, we derive sufficient conditions on p and w to satisfy Equation (7). We decompose the problem in two parts. To satisfy the first inclusion, i.e. ensure that P is a loop invariant, it suffices to guarantee that F (P ) ⊆ P , namely that P is an inductive invariant. Using Equation (5), P is an inductive invariant if and only if:

X in ∪ i∈I T i P ∩ X i ∩ X 0 ⊆ P ,
or equivalently:

X in ⊆ P , ∀ i ∈ I, T i P ∩ X i ∩ X 0 ⊆ P . ( 8
)
Thus, we obtain:

p(x) ≤ 0 , ∀x ∈ X in , ∀ i ∈ I , p (T i (x)) ≤ 0 , ∀x ∈ P ∩ X i ∩ X 0 . ( 9
)
Now, we are interested in the second and third inclusions at Equation (7) that is the sublevel property satisfaction. The condition P ⊆ P κ,w ⊆ P κ,α can be formulated as follows:

κ(x) ≤ w ≤ α , ∀x ∈ P . ( 10
)
We recall that we have supposed that P is written as

{x ∈ R d | p(x) ≤ 0} where p ∈ R[x]
. Finally, we provide sufficient conditions to satisfy both (9) and (10).

Consider the following optimization problem:

       inf p∈R[x],w∈R w , s.t. p(x) ≤ 0 , ∀x ∈ X in , ∀ i ∈ I , p (T i (x)) ≤ p(x) , ∀x ∈ X i ∩ X 0 , κ(x) ≤ w + p(x) , ∀x ∈ R d . ( 11
)
We remark that α is not present in Problem (11). Indeed, since we minimize w, either there exists a feasible w such that w ≤ α and we can exploit this solution or such w is not available and we cannot conclude. However, from Problem (11), we can extract (p, w) and in the case where the optimal bound w is greater than α, we could use this solution with another method such as policy iteration [START_REF] Adjé | A Sums-of-Squares Extension of Policy Iterations[END_REF]. In practice, we rely on sum-of-squares programming to solve a relaxed version of Problem (11).

Sums-of-Squares Programming for Invariant Generation

We first recall some basic background about sums-of-squares certificates for polynomial optimization. Let R[x] 2m stands for the set of polynomials of degree at most 2m and Σ[x] ⊂ R[x] be the cone of sums-of-squares (SOS) polynomials, that is Σ[x] := { i q 2 i , with q i ∈ R[x] }. Our work will use the simple fact that for all p ∈ Σ[x], then p(x) ≥ 0 for all x ∈ R d i.e. Σ[x] is a restriction of the set of the nonnegative polynomials. For q ∈ R[x] 2m , finding a SOS decomposition q = i q 2 i valid over R d is equivalent to solve the following matrix linear feasibility problem:

q(x) = b m (x) T Q b m (x) , ∀x ∈ R d , ( 12
)
where b m (x) := (1, x 1 , . . . , x d , x 2 1 , x 1 x 2 , . . . , x m d ) (the vector of all monomials in x up to degree m) and Q being a semidefinite positive matrix (i.e. all the eigenvalues of Q are nonnegative). The size of Q (as well as the length of b m ) is

d+m d . Example 3. consider the bi-variate polynomial q(x) := 1 + x 2 1 -2x 1 x 2 + x 2 2 . With b 1 (x) = (1, x 1 , x 2 ),
one looks for a semidefinite positive matrix Q such that the polynomial equality q(x) = b 1 (x) T Q b 1 (x) holds for all x ∈ R 2 . The matrix

Q =   1 0 0 0 1 -1 0 -1 1  
satisfies this equality and has three nonnegative eigenvalues, which are 0, 1, and 2, respectively associated to the three eigenvectors e 0 := (0, 1, 1) , e 1 := (1, 0, 0) and e 2 := (0, 1, -1) . Defining the matrices L := (e 1 e 2 e 0 ) = 1 0 0 0 1 1 0 -1 1 and D = 1 0 0 0 2 0 0 0 0 , one obtains the decomposition Q = L D L and the equality

q(x) = (L b 1 (x)) T D (L b 1 (x)) = σ(x) = 1 + (x 1 -x 2 ) 2 , for all x ∈ R 2 .
The polynomial σ is called a SOS certificate and guarantees that q is nonnegative.

In practice, one can solve the general problem (12) by using semidefinite programming (SDP) solvers (e.g. Mosek [AA00], SDPA [YFN + 10]). For more details about SDP, we refer the interested reader to [START_REF] Vandenberghe | Semidefinite programming[END_REF].

One way to strengthen the three nonnegativity constraints of Problem ( 11) is to consider the following hierarchy of SOS programs, parametrized by the integer m representing the half of the degree of p: Proof. The feasible solution (p m , w m ) is associated with SOS certificates ensuring that the three equality constraints of Problem (13) hold: {σ 0 , σ j } is associated to the first one, {σ i , µ i j , γ i j } is associated to the second one and ψ is associated to the third one. The first equality constraint, namely

                                                               inf p∈R[x]2m,w∈R w , s.t. -p = σ 0 - nin j=1 σ j r in j , ∀ i ∈ I, p -p • T i = σ i - ni j=1 µ i j r i j - n0 j=1 γ i j r 0 j , w + p -κ = ψ , ∀ j = 1, . . . , n in , σ j ∈ Σ[x] , deg(σ j r in j ) ≤ 2m , σ 0 ∈ Σ[x] , deg(σ 0 ) ≤ 2m , ∀ i ∈ I , σ i ∈ Σ[x] , deg(σ i ) ≤ 2m deg T i , ∀ i ∈ I , ∀ j = 1, . . . , n i , µ i j ∈ Σ[x] , deg(µ i j r i j ) ≤ 2m deg T i , ∀ i ∈ I , ∀ j = 1, . . . , n 0 , γ i ∈ Σ[x] , deg(γ i j r 0 j ) ≤ 2m deg T i , ψ ∈ Σ[x] , deg(ψ) ≤ 2m . ( 13
-p m (x) = σ 0 (x) - nin j=1 σ j (x)r in j (x) , ∀x ∈ R d , implies that ∀x ∈ X in , p m (x) ≤ 0. Similarly, one has ∀i ∈ I, ∀x ∈ X i ∩ X 0 , p m (T i (x)) ≤ p m (x) and ∀x ∈ R d , κ(x) ≤ w m + p m (x). Then (p m , w m )
is a feasible solution of Problem (11). The second statement comes directly from Lemma 1.

Computational considerations. Define t := max{deg T i , i ∈ I}. At step m of this hierarchy, the number of SDP variables is proportional to d+2mt d and the number of SDP constraints is proportional to d+mt d . Thus, one expects tractable approximations when the number d of variables (resp. the degree 2m of the template p) is small. However, one can handle bigger instances of Problem ( 13) by taking into account the system properties. For instance one could exploit sparsity as in [START_REF] Waki | Sums of squares and semidefinite programming relaxations for polynomial optimization problems with structured sparsity[END_REF] by considering the variable sparsity correlation pattern of the polynomials {T i , i ∈ I}, {r i j , i ∈ I, j = 1, . . . , n i }, {r 0 j , j = 1, . . . , n 0 }, {r in j , j = 1, . . . , n in } and κ.

Benchmarks

Here, we perform some numerical experiments while solving Problem (13) (given in Section 4) on several examples. In Section 5.1, we verify that the program of Example 1 satisfies some boundedness property. We also provide examples involving higher dimensional cases. Then, Section 5.2 focuses on other properties, such as checking that the set of variable values avoids an unsafe region. Numerical experiments are performed on an Intel Core i5 CPU (2.40 GHz) with Yalmip being interfaced with the SDP solver Mosek.

Checking boundedness of the set of variables values

Example 4. Following Example 1, we consider the constrained piecewise discrete-

time dynamical system S = (X in , X 0 , {X 1 , X 2 }, {T 1 , T 2 }) with X in = [0.9, 1.1]× [0, 0.2], X 0 = {x ∈ R 2 | r 0 (x) ≤ 0} with r 0 : x → -1, X 1 = {x ∈ R 2 | r 1 (x) ≤ 0} with r 1 : x → x 2 -1, X 2 = {x ∈ R 2 | r 2 (x) < 0} with r 2 = -r 1 and T 1 : (x 1 , x 2 ) → (c 11 x 2 1 + c 12 x 3 2 , c 21 x 3 1 + c 22 x 2 2 ), T 2 : (x 1 , x 2 ) → (d 11 x 3 1 + d 12 x 2 2 , d 21 x 2 1 + d 22 x 2 2
). We are interested in showing that the boundedness property P • 2 2 ,α holds for some positive α. Here we illustrate the method by instantiating the program of Example 1 with the following input: a 1 = 0.9, a 2 = 1.

(i) 1 , x (i)
2 ) (i = 1, . . . , N ) inside the box X in = [0.9, 1.1] × [0, 0.2] (see the corresponding square of dots on Figure 2). The other dots are obtained after successive updates of each point (x

(i) 1 , x (i)
2 ) by the program of Example 1. The sets of dots in Figure 2 are obtained with N = 100 and six successive iterations.

At step m = 3, Program (13) yields a solution (p 3 , w 3 ) ∈ R 6 [x] × R together with SOS certificates, which guarantee the boundedness property, that is

x ∈ R =⇒ x ∈ P 3 := {p 3 (x) ≤ 0} ⊆ P • 2 2 ,w3 =⇒ x 2 2 ≤ w 3 . One has p 3 (x) := -2.510902467-0.0050x 1 -0.0148x 2 +3.0998x 2 1 -0.8037x 3 2 -3.0297x 3 1 +2.5924x 2 2 + 1.5266x 1 x 2 -1.9133x 2 1 x 2 -1.8122x 1 x 2 2 + 1.6042x 4 1 + 0.0512x 3 1 x 2 -4.4430x 2 1 x 2 2 - 1.8926x 1 x 3 2 + 0.5464x 4 2 -0.2084x 5 1 + 0.5866x 4 1 x 2 + 2.2410x 3 1 x 2 2 + 1.5714x 2 1 x 3 2 - 0.0890x 1 x 4 2 -0.9656x 5 2 + 0.0098x 6 1 -0.0320x 5 1 x 2 -0.0232x 4 1 x 2 2 + 0.2660x 3 1 x 3 2 + 0.7746x 2 1
x 4 2 +0.9200x 1 x 5 2 +0.6411x 6 2 (for the sake of conciseness, we do not display p 4 and p 5 ).

Figure 2 displays in light gray outer approximations of the set of possible values X 1 taken by the program of Example 4 as follows: (a) the degree six sublevel set P 3 , (b) the degree eight sublevel set P 4 and (c) the degree ten sublevel set P 5 . The outer approximation P 3 is coarse as it contains the box [-1.5, 1.5] 2 . However, solving Problem (13) at higher steps yields tighter outer approximations of R together with more precise bounds w 4 and w 5 (see the corresponding row in Table 2).

We also succeeded to certify that the same property holds for higher dimensional programs, described in Example 5 (d = 3) and Example 6 (d = 4).

Example 5. Here we consider

X in = [0.9, 1.1] × [0, 0.2] 2 , r 0 : x → -1, r 1 : x → x 2 2 -1, r 2 = -r 1 , T 1 : (x 1 , x 2 , x 3 ) → 1/4(0.8x 2 1 + 1.4x 2 -0.5x 2 3 , 1.3x 1 + 0.5x 2 3 , 1.4x 2 +0.8x 2 3 ), T 2 : (x 1 , x 2 , x 3 ) → 1/4(0.5x 1 +0.4x 2 2 , -0.6x 2 2 +0.3x 2 3 , 0.5x 3 + 0.4x 2 1 ) and κ : x → x 2 2 .
Example 6. Here we consider

X in = [0.9, 1.1] × [0, 0.2] 3 , r 0 : x → -1, r 1 : x → x 2 2 -1, r 2 = -r 1 , T 1 : (x 1 , x 2 , x 3 , x 4 ) → 0.25(0.8x 2 1 + 1.4x 2 -0.5x 2 3 , 1.3x 1 + 0.5, x 2 2 -0.8x 2 4 , 0.8x 2 3 + 1.4x 4 , 1.3x 3 + 0.5x 2 4 ), T 2 : (x 1 , x 2 , x 3 , x 4 ) → 0.25(0.5x 1 + 0.4x 2 2 , -0.6x 2 1 + 0.3x 2 2 , 0.5x 3 + 0.4x 2 4 , -0.6x 3 + 0.3x 2 4 ) and κ : x → x 2 2 .
Table 1 reports several data obtained while solving Problem (13) at step m, (2 ≤ m ≤ 5), either for Example 4, Example 5 or Example 6. Each instance of Problem ( 13) is recast as a SDP program, involving a total number of "Nb. vars" SDP variables, with a SDP matrix of size "Mat. size". We indicate the CPU time required to compute the optimal solution of each SDP program with Mosek.

The symbol "-" means that the corresponding SOS program could not be solved within one day of computation. These benchmarks illustrate the computational considerations mentioned in Section 4 as it takes more CPU time to analyze higher dimensional programs. Note that it is not possible to solve Problem (13) at step 5 for Example 6. A possible workaround to limit this computational blow-up would be to exploit the sparsity of the system.

Other properties

Here we consider the program given in Example 7. One is interested in showing that the set X 1 of possible values taken by the variables of this program does not meet the ball B of center (-0.5, -0.5) and radius 0.5.

Example 7. Let consider the PPS S = (X in , X 0 , {X 1 , X 2 }, {T 1 , T 2 }) with X in = [0.5, 0.7] × [0.5, 0.7], X 0 = {x ∈ R 2 | r 0 (x) ≤ 0} with r 0 : x → -1, X 1 = {x ∈ R 2 | r 1 (x) ≤ 0} with r 1 : x → x 2 2 -1, X 2 = {x ∈ R 2 | r 2 (x) ≤ 0} with r 2 = -r 1 and T 1 : (x 1 , x 2 ) → (x 2 1 + x 3 2 , x 3 1 + x 2 2 ), T 2 : (x, y) → (0.5x 3 1 + 0.4x 2 2 , -0.6x 2 1 + 0.3x 2 2 ). With κ : (x 1 , x 2 ) → 0.25 -(x 1 + 0.5) 2 -(x 2 + 0.5) 2 , one has B := {x ∈ R 2 | 0 ≤ κ(x)}.
Here, one shall prove x ∈ R =⇒ κ(x) < 0 while computing some negative α such that R ⊆ P κ,α . Note that κ is not a norm, by contrast with the previous examples.

At step m = 3 (resp.m = 4), Program (13) yields a nonnegative solution w 3 (resp. w 4 ). Hence, it does not allow to certify that R ∩ B is empty. This is illustrated in both Figure 3 Let S be the PPS (X in , X

0 , {X 1 , X 2 }, {T 1 , T 2 }) with X in = [-1, 1] × [-1, 1], X 0 = {x ∈ R 2 | r 0 (x) ≤ 0} with r 0 : x → -1, X 1 = {x ∈ R 2 | r 1 (x) ≤ 0} with r 1 : x → x 2 -x 1 , X 2 = {x ∈ R 2 | r 2 (x) ≤ 0} with r 2 = -r 1 and T 1 : (x 1 , x 2 ) → (0.687x 1 + 0.558x 2 -0.0001 * x 1 x 2 , -0.292x 1 + 0.773x 2 ), T 2 : (x, y) → (0.369x 1 + 0.532x 2 -0.0001x 2 1 , -1.27x 1 +0.12x 2 -0.0001x 1 x 2 ).
We consider the boundedness property κ 1 := • 2 2 as well as κ 2 (x) := T 1 (x) -T 2 (x) 2 2 . The function κ 2 can be viewed as the absolute error made by updating the variable x after a possibly "wrong" branching. Such behaviors could occur while computing wrong values for the conditionals (e.g. r 1 ) using floating-point arithmetics. Table 2 indicates the hierarchy of bounds obtained after solving Problem (13) with m = 3, 4, 5, for both properties. The bound w 5 = 2.84 (for κ 1 ) implies that the set of reachable values may not be included in the initial set X in . A valid upper bound of the error function κ 2 is given by w 5 = 2.78.

We finally present further use of the set P defined at Equation ( 7). This sublevel set can be viewed as a template abstraction, following from the definition in [START_REF] Adjé | Coupling policy iteration with semi-definite relaxation to compute accurate numerical invariants in static analysis[END_REF], with a fixed template basis p and an associated 0 bound. This representation allows to develop a policy iteration algorithm [START_REF] Adjé | A Sums-of-Squares Extension of Policy Iterations[END_REF] to obtain more precise inductive invariants.

We now give some simple method to complete this template basis to improve the precision of the bound w found with Problem (13).

Proposition 2 (Template basis completions). Let (p, w) be a solution of Problem (13) and Q be a finite subset of R

[x] such that for all q ∈ Q, p-q ∈ Σ[x]. Then R ⊆ {x ∈ R d | p(x) ≤ 0, q(x) ≤ 0, ∀ q ∈ Q} ⊆ P κ,w ⊆ P κ,α and {x ∈ R d | p(x) ≤ 0, q(x) ≤ 0, ∀ q ∈ Q} is an inductive invariant. Proof. Let Q be the set {x ∈ R d | p(x) ≤ 0, q(x) ≤ 0, ∀ q ∈ Q}. It is obvious that Q ⊆ P = {x ∈ R d | p(x) ≤ 0}
and hence Q ⊆ P κ,w . Now let us prove that Q is an inductive invariant. We have to prove that Q satisfies Equation (8) that is: (i) For all x ∈ X in , q(x) ≤ 0; (ii) For all i ∈ I, for all x ∈ Q ∩ X i ∩ X 0 , q(T i (x)) ≤ 0. For all q ∈ Q, we denote by ψ q the element of Σ[x] such that p -q = ψ q . Let us show (i) and let x ∈ X in . We have q(x) = p(x) -ψ q (x) and since ψ q ∈ Σ[x], we obtain, q(x) ≤ p(x). Now from Proposition 1 and Lemma 1 and since (p, w) is a solution of Problem (13), we conclude that q(x) ≤ p(x) ≤ 0. Now let us prove (ii) and let i ∈ I and x ∈ Q ∩ X i ∩ X 0 . We get q(T i (x)) = p(T i (x)) -ψ q (T i (x)) and since ψ q ∈ Σ[x], we obtain q(T i (x)) ≤ p(T i (x)). Using the fact that (p, w) is a solution of Problem (13) and using Proposition 1 and Lemma 1, we obtain q(T i (x)) ≤ p(T i (x)) ≤ p(x). Since x ∈ Q ⊆ P = {y ∈ R d | p(y) ≤ 0}, we conclude that q(T i (x)) ≤ 0.

Actually, we can weaken the hypothesis of Proposition 2 to construct an inductive invariant. Indeed, after the computation of p following Problem (13), it suffices to take a polynomial q such that p -q ≥ 0. Nevertheless, we cannot compute easily such a polynomial q. By using the hypothesis p -q ∈ Σ[x], we can compute q by sum-of-squares. Proposition 2 allows to define a simple method to construct a basic semialgebraic inductive invariant set. Then the polynomials describing this basic semialgebraic set defines a new templates basis and this basic semialgebraic set can be used as initialisation of the policy iteration algorithm developed in [START_REF] Adjé | A Sums-of-Squares Extension of Policy Iterations[END_REF]. Note that the link between the templates generation and the initialisation of policy iteration has been addressed in [START_REF] Adjé | Policy iteration in finite templates domain[END_REF].

Example 9. Let us consider the property P • 2 2 ,∞ and let (p, w) be a solution of Problem (13). We have κ(x) = 1≤j≤k x 2 j and w + p -κ = ψ where ψ ∈ Σ[x]. In [START_REF] Roux | A generic ellipsoid abstract domain for linear time invariant systems[END_REF], the templates basis used to compute bounds on the reachable values set consists in the square variables plus a Lyapunov function. Let us prove that, in our setting, Q = {x 2 k -w, k = 1, . . . , d} can complete {p} in the sense of Proposition 2. Let k ∈ {1, . . . , d} and let

x ∈ R d , p(x) -(x 2 k -w) = p(x) -κ(x) + w + j =k x 2 j = ψ(x) + j =k x 2 j ∈ Σ[x].
Roux et al. [START_REF] Roux | A generic ellipsoid abstract domain for linear time invariant systems[END_REF] provide an automatic method to compute floating-point certified Lyapunov functions of perturbed affine loop body updates. They use Lyapunov functions with squares of coordinate functions as quadratic invariants in case of single loop programs written in affine arithmetic. In the context of hybrid systems, certified inductive invariants can be computed by using SOS approximations of parametric polynomial optimization problems [START_REF] Lin | Exact safety verification of hybrid systems using sums-of-squares representation[END_REF]. In [START_REF] Prajna | Safety verification of hybrid systems using barrier certificates[END_REF], the authors develop a SOS-based methodology to certify that the trajectories of hybrid systems avoid an unsafe region.

In the context of static analysis for semialgebraic programs, the approach developed in [START_REF] Cousot | Proving Program Invariance and Termination by Parametric Abstraction, Lagrangian Relaxation and Semidefinite Programming[END_REF] focuses on inferring valid loop/conditional invariants for semialgebraic programs c . This approach relaxes an invariant generation problem into the resolution of nonlinear matrix inequalities, handled with semidefinite programming. Our method bears similarities with this approach but we generate a hierarchy of invariants (of increasing degree) with respect to target polynomial properties and restrict ourselves to linear matrix inequality formulations. In [START_REF] Bagnara | Generation of basic semi-algebraic invariants using convex polyhedra[END_REF], invariants are given by polynomial inequalities (of bounded degree) but the method relies on a reduction to linear inequalities (the polyhedra domain). Template polyhedra domains allow to analyze reachability for polynomial systems: in [START_REF] Amin | Reachability analysis of polynomial systems using linear programming relaxations[END_REF], the authors propose a method that computes linear templates to improve the accuracy of reachable set approximations, whereas the procedure in [START_REF] Dang | Reachability analysis for polynomial dynamical systems using the bernstein expansion[END_REF] relies on Bernstein polynomials and linear programming, with linear templates being fixed in advance. Bernstein polynomials also appear in [START_REF] Roux | A polynomial template abstract domain based on bernstein polynomials[END_REF] as polynomial templates but they are not generated automatically. In [START_REF] Srivastava | Program verification using templates over predicate abstraction[END_REF], the authors use SMT-based techniques to automatically generate templates which are defined as formulas built with arbitrary logical structures and predicate conjunctions. Other reductions to systems of polynomial equalities (by contrast with polynomial inequalities, as we consider here) were studied in [START_REF] Müller-Olm | Computing polynomial program invariants[END_REF][START_REF] Rodríguez-Carbonell | Automatic generation of polynomial invariants of bounded degree using abstract interpretation[END_REF] and more recently in [START_REF] Cachera | Inference of polynomial invariants for imperative programs: A farewell to gröbner bases[END_REF].

In this paper, we give a formal framework to relate the invariant generation problem to the property to prove on analyzed program. We proposed a practical method to compute such invariants in the case of polynomial arithmetic using sums-of-squares programming. This method is able to handle non trivial examples, as illustrated through the numerical experiments. Topics of further investigation include refining the invariant bounds generated for a specific sublevel property, by applying the policy iteration algorithm. Such a refinement would be of particular interest if one can not decide whether the set of variable values avoids an unsafe region when the bound of the corresponding sums-ofsquares program is not accurate enough. For the case of boundedness property, it would allow to decrease the value of the bounds on the variables. Finally, our method could be generalized to a larger class of programs, involving semialgebraic or transcendental assignments, while applying the same polynomial reduction techniques as in [START_REF] Victor Magron | Certification of Real Inequalities -Templates and Sums of Squares[END_REF].

Lemma 1 .

 1 Let (p, w) be any feasible solution of Problem (11) with w ≤ α or w < ∞ in the case of α = ∞. Then (p, w) satisfies both (9) and (10) with P := {x ∈ R d | p(x) ≤ 0}. Finally, P and P κ,w satisfy Equation (7).

) Proposition 1 .

 1 For a given m ∈ N, let (p m , w m ) be any feasible solution of Problem (13). Then (p m , w m ) is also a feasible solution of Problem (11). Moreover, if w m ≤ α then both P m := {x ∈ R d | p m (x) ≤ 0} and P κ,wm satisfy Equation (7).

Fig. 2 .

 2 Fig. 2. A hierarchy of sublevel sets Pm for Example 4

  (a) and Figure 3 (b), where the light grey region does not avoid the ball B. However, solving Program (13) at step m = 5 yields a negative bound w 5 together with a certificate that R avoids the ball B (see Figure 3 (c)). The corresponding values of w m (m = 3, 4, 5) are given in Table2.

Fig. 3 .

 3 Fig. 3. A hierarchy of sublevel sets Pm for Example 7

Table 1 .

 1 Comparison of timing results for Example 4, 5 and 6

	Degree 2m	4	6	8	10
	Example 4	Nb. vars Mat. size	1513 368	5740 802	15705 1404	35212 2174
	(d = 2)	Time	0.82 s	1.35 s	4.00 s	9.86 s
	Example 5	Nb. vars Mat. size	2115 628	11950 1860	46461 4132	141612 7764
	(d = 3)	Time	0.84 s	2.98 s	21.4 s	109 s
	Example 6	Nb. vars Mat. size	7202 1670	65306 6622	18480 373057	--
	(d = 4)	Time	2.85 s	57.3 s	1534 s	-

Table 2 .

 2 Hierarchies of bounds obtained for various properties

	Benchmark	κ	w2	w3	w4	w5
	Example 4	• 2 2	639	17.4	2.44	2.02
	Example 7	x → 0.25 -x + 0.5 2 2	0.25 0.249 0.0993 -0.0777
	Example 8	• 2 2 x → T 1 (x) -T 2 (x) 2 2	10.2 5.66	2.84 2.81	2.84 2.78	2.84 2.78

c This approach also handles semialgebraic program termination

a The author is supported by the RTRA /STAE Project BRIEFCASE and the ANR ASTRID VORACE Project. b The author is supported by EPSRC (EP/I020457/1) Challenging Engineering Grant.


Modeling, Control and Design Optimization for a Fully-actuated Hexarotor Aerial Vehicle with Tilted Propellers

Sujit Rajappa 1 , Markus Ryll 3 , Heinrich H. Bülthoff 1,2 and Antonio Franchi 3,4 Abstract-Mobility of a hexarotor UAV in its standard configuration is limited, since all the propeller force vectors are parallel and they achieve only 4 DoF actuation, similar, e.g., to quadrotors. As a consequence, the hexarotor pose cannot track an arbitrary trajectory over time. In this paper, we consider a different hexarotor architecture where propellers are tilted, without the need of any additional hardware. In this way, the hexarotor possess a 6 DoF actuation which allows to independently reach positions and orientations in free space and to be able to exert forces on the environment to resist any wrench for aerial manipulation tasks. After deriving the dynamical model of the proposed hexarotor, we discuss the controllability and the tilt angle optimization to reduce the control effort. An exact feedback linearization and decoupling control law is proposed based on the input-output mapping, considering the Jacobian and task acceleration, for non-linear trajectory tracking. The robustness of our approach is validated by simulation results.

I. INTRODUCTION

Research in the field and applications related to unmanned aerial vehicles (UAVs) has been a popular research topic in recent times, see e.g., [START_REF] Franchi | Shared control: Balancing autonomy and human assistance with a group of quadrotor UAVs[END_REF] and references therein. The application possibility to use the UAV for various tasks such as search and rescue operation, exploration, surveillance, cooperative swarm tasks or transportation are all increasing and has been the main research subject with growing interest in the last decade with many industrial collaborations. Lately, the mobile manipulation tasks by aerial vehicles and physical interaction with the environment for various applications has been growing ground within the UAV community. The interaction can be done by direct contact [START_REF] Marconi | Modeling and control of a flying robot interacting with the environment[END_REF]- [START_REF]A nonlinear force observer for quadrotors and application to physical interactive tasks[END_REF], by considering simple grasping/manipulation tasks [START_REF] Pounds | Grasping from the air: Hovering capture and load stability[END_REF], [START_REF] Lindsey | Construction of cubic structures with quadrotor teams[END_REF] and has moved forward to multiple collaborative interactive UAVs [START_REF] Gioioso | The Flying Hand: a formation of uavs for cooperative aerial telemanipulation[END_REF].

Among the many challenges faced by typical UAVs, such as little flight time, limited payload capacity, uncertainties in outside environment etc., an important one is the underactuation, i.e., the inability to exert forces in some directions of the body frame. Quadrotors have been used as the main platform for applications as well as research, though they are also underactuated, i.e., they cannot exert any force parallel to the plane perpendicular to their vertical direction in body frame. This is why a quadrotor needs to roll and pitch to move in any direction.

But when it comes to physical interaction, underactuation might become a serious problem for the capabilities and overall stabilization of the aerial vehicle. As the application complexity is going higher, major breakthroughs and advancements in innovative mechanical designs, actuation concepts, micro-electro mechanical systems, sensor technology and power capacity is always envisioned. Several possibilities have been proposed in the past literature spanning different concepts: ducted-fan designs [START_REF] Naldi | Design and experimental validation of a nonlinear control law for a ducted-fan miniature aerial vehicle[END_REF], tilt-wing mechanisms [START_REF] Oner | Dynamic model and control of a new quadrotor unmanned aerial vehicle with tilt-wing mechanism[END_REF], or tilt-rotor actuations [START_REF] Kendoul | Modeling and control of a small autonomous aircraft having two tilting rotors[END_REF], [START_REF] Sanchez | Autonomous hovering of a noncyclic tiltrotor UAV: Modeling, control and implementation[END_REF]. The concept of tilt-rotor architecture has been much explored to increase flight time [START_REF] Flores | Quad-tilting rotor convertible MAV: Modeling and real-time hover flight control[END_REF] but not for the improvement of the underactuation problem. In [START_REF] Salazar | Modeling and real-time stabilization of an aircraft having eight rotors[END_REF] the underactuation was addressed by four additional rotors at the end of each frame in lateral position. But the position of the rotor increased the complexity of controllability because of the air flow between the vertical and the lateral rotors, resulting in non-linear dynamics along with the increase in payload.

Our own in-lab investigation led to the novel quadrotor design [START_REF] Ryll | Modeling and control of a quadrotor UAV with tilting propellers[END_REF] with tilted propellers by 4 additional actuators included for the tilting thereby creating the possibility to regulate independently the 6 DoFs of the platform. Though underactuation problem was solved by this design, the need of tilting the propellers in order to resist to any external wrench makes it tough for the aerial manipulation task, where forces shall be exerted instantaneously to resist to unexpected external wrenches. Additionally, the use of servomotors for tilting the propellers makes the overall model challenging to control in real scenarios involving physical interaction.

In [START_REF] Voyles | Hexrotor UAV platform enabling dextrous interaction with structurespreliminary work[END_REF], a hexarotor with the propellers rotated about one axis was suggested. Our approach is different, as we present a more general tilt design, a control law for 6 DoFs trajectory tracking, a methodology to optimize the fixed tilting angles for each propeller depending on the task in exam and an improved mechanical design where all the propellers lie in the same plane.

Taking inspiration from all the related work, we propose a novel hexarotor with tilted propeller design, where each rotor is fixedly mounted in a configuration that is rotated about two possible axes. The main objective of this work is full controllability of the UAV's position and orientation by means of tilted propellers, thereby making it completely actuated. The full actuation comes with the cost of a more advanced mechanical design.

The focus and structure of this paper is therefore: (i)

to discuss in detail and derive the dynamic model for the proposed hexarotor in Sec. II, (ii) to devise and develop the closed-loop controller for the hexarotor which is able to asymptotically track an arbitrary desired trajectory for the position and orientation in 3 dimensional free space in Sec. III, (iii) to optimize the propeller tilt angles depending on the application/trajectory to reduce the overall control effort in the Sec. IV, (iv) to start the design of a novel and feasible hexarotor architecture with tilted propellers in Sec. V, and finally (v) to test the hexarotor model and its theoretical concepts in simulation, presented in Sec. VI. Conclusions and future perspectives are given in Sec. VII.

II. DESIGN AND MODELING

A standard hexarotor possesses six propellers that are all rotating about six parallel axes. Even though this choice increases redundancy and payload, such configuration has an underactuated dynamics similar to a standard quadrotor. In fact, the six propellers create an input force that is always parallel to that axis, no matter the values of the six rotational speeds. In this case a change of the direction of the input force can only be obtained by reorienting the whole vehicle. Therefore, the desired output trajectory can only be defined by a 4-dimensional output, namely the center of mass (CoM) 3D position plus the yaw angle, despite the presence of 6 control inputs. In fact in [START_REF] Mistler | Exact linearization and noninteracting control of a 4 rotors helicopter via dynamic feedback[END_REF] it has been proven that such kind of systems are exactly linearizable with a dynamic feedback using as linearizing output, i.e. the CoM position and the yaw angle. Feedback linearizability also implies that the system is differentially flat and that the linearizing output is the flat output [START_REF] De Luca | Trajectory planning and control for planar robots with passive last joint[END_REF]. The remaining two configuration variables, i.e., the roll and pitch angles, cannot be chosen at will, since they are being determined by the desired trajectory of the CoM, the yaw angle, and their derivatives.

On the converse, the goal of the hexarotor modeling approach presented here is to exploit at best the six available inputs, thus resulting in a system that is fully actuated, i.e., linear and angular accelerations can be set independently acting on the six inputs. In order to obtain full actuation, we remove the constraint for the propellers to rotate about six parallel axes, so that a force in any direction can be generated regardless of the vehicle orientation. Thanks to full actuation, this hexarotor can follow desired 6-DoFs trajectories comprising both the CoM position and, independently, the vehicle orientation described, e.g., by roll, pitch, and yaw, or by a rotation matrix.

Even though a reallocation and reorientation of the six propellers allows for more design flexibility it also increases the number of design parameters thus increasing the design complexity. In order to find a good compromise between full actuation and low number of model parameters, we decide to add some constraints on the parameters, i.e.: • the tilted planes of the propellers are coordinated, in the sense that, each pair of propellers ((1)-( 2), ( 3)-( 4), ( 5)-( 6)) are tilted oppositely about the X-and Y-axis. Detailed explanation can be found in Section IV; These three design choices simplify the design complexity while still allowing a full spectrum of actuation capabilities, as it will be shown in the paper. In the following we formally describe our hexarotor model and all the related design choices.

A. Static System Description

We denote the world inertial frame with

F W : {O W - X W Y W Z W } and with F B : {O B -X B Y B Z B }
the body frame attached to the hexarotor frame, where O B coincides with the hexarotor CoM. Let the frame associated with the i-th propeller be defined as F Pi : {O Pi -X Pi Y Pi Z Pi }, where i = 1 . . . 6. The origin O Pi coincides with the center of spinning and the CoM of the i-th propeller, the axes X Pi and Y Pi define the rotation plane of the propeller, and Z Pi is the axis about which the propeller spins and coincides with the generated thrust force. The propeller frame F Pi is rigidly attached to the hexarotor frame, rather than to the propeller, which spins about Z Pi . In fact, only the direction of the force and torque exerted by the propeller are relevant to our problem. The actual spinning angle of each propeller is not important for the motion, as it will be explained in Sec. II-B.

We shall denote simply by p ∈ R 3 the position of O B in F W , and by B p i ∈ R 3 the position of O Pi in F B , with i = 1 . . . 6. We make the following design choice:

B p i = R Z (i -1) π 3   L xi 0 0   , ∀i = 1 . . . 6 (1)
where R Z (•) is the canonical rotation matrix about a Zaxis and

L xi > 0 is the distance between O Pi and O B .
It is important to note that the parameter L xi is chosen depending on the strength of propeller, size of hexarotor, payload needs and doesn't affect the full-actuation of hexarotor. Furthermore, we chose

L x1 = L x4 , L x2 = L x5 and L x3 = L x6 . With these choices of R Z (•)
Let the matrix W R B ∈ SO(3) represent the ori- entation of F B F W and B R Pi ∈ SO(3) the orientation of F Pi w.r.t. F B , for i = 1 . . . 6.
In order to obtain a minimal parameterization of the propeller orientation we decompose each B R Pi in three consecutive rotations

B R Pi R Z (i-1)π 3 R X (α i )R Y (β i ), ∀i = 1 . . . 6 (2)
where the angular parameters α i and β i represent the tilt angles, that uniquely define the rotation plane of the i-th propeller, X Pi Y Pi or, equivalently, the direction of Z Pi in F B . The angles α i and β i have a clear geometrical interpretation, in fact the i-th propeller plane X Pi Y Pi is obtained from X B Y B by first applying a rotation of α i about the line O B O Pi and then a rotation of β i about Y Pi , which lies on X B Y B and is perpendicular to O B O Pi . The α i and β i rotation is pictorially represented in Fig. 2.

For convenience, we group the tilting angles in two 6-tuples: α = (α 1 , α 2 , α 3 , α 4 , α 5 , α 6 ) and β = (β 1 , β 2 , β 3 , β 4 , β 5 , β 6 ). Considering the design constrains explained earlier in Section II, the tilt w.r.t. X Pi is fixed as

α 1 = -α 2 = α 3 = -α 4 = α 5 = -α 6 and w.r.t Y Pi is β 1 = -β 2 = β 3 = -β 4 = β 5 = -β 6 in a co-ordinated way.
The angle of tilt is further detailed in optimization (Section IV).

In this paper we consider the case in which L xi , α i , β i , for i = 1 . . . 6, are constant during flight, but can be changed during a pre-flight setup, in order, e.g., to minimize the sum of the overall control effort for a specific task, as shown in Sec. IV.

B. Equations of Motion

Utilizing the standard Newton-Euler approach for dynamic systems, it is possible to derive the complete dynamic equations of the hexarotor by considering the forces and torques that are generated by each propeller rotation together with the significant gyroscopic and inertial effects. In the following we recap the standard 1 assumptions that we are considering here:

• O B coincides with the CoM of the hexarotor;

• O Pi coincides with the CoM of the i-th propeller;

• the motors actuating the six propellers implement a fast high-gain local controller which is able to impose a desired spinning speed with negligible transient, thus allowing to consider the spinning rates as (virtual) control inputs in place of the motor torques; • gyroscopic and inertial effects due to the propellers and the motors are considered as second-order disturbances to be rejected by the feedback nature of the controller; • the tilted propellers might cause additional turbulences due to the possible intersection of the airflows. These turbulences are considered as negligible as the possible intersection of the airflows happens not close to the propellers. In fact, tilt configurations have been already proven to be feasible in reality [START_REF] Ryll | First flight tests for a quadrotor UAV with tilting propellers[END_REF]. We will test in simulation (see Section VI) the practicability of these assumptions with the proposed controller on a dynamic model which includes the aforementioned unmodeled effects.

For ease of presentation, in the following we shall express the translational dynamics in F W where as the rotational dynamics is expressed in F B .

1) Rotational dynamics: Denote with ω B ∈ R 3 the angular velocity of F B , with respect to F W , expressed in F B . Then the rotational dynamics is

I B ωB = -ω B × I B ω B + τ + τ ext , (3)
where I B is the hexarotor body inertia matrix, τ ext accounts for external disturbances and unmodeled effects, and τ is the input torque, which is decomposed in

τ = τ thrust + τ drag , (4)
where τ thrust is produced by the six propeller thrusts and τ drag is due to the six propeller drags. The two individual components of (4) are discussed in detail below. a) Torque due to thrusts (τ thrust ): The i-th propeller creates a force vector applied at O Pi and directed along Z Pi , which is expressed in F Pi by

T thrusti = 0 0 k f ω2 i T (5)
where k f > 0 is a constant thrust coefficient and ωi is the spinning velocity of the i-th propeller. The thrust torque, expressed in F B is then

τ thrust = 6 i=1 B p i × B R Pi T thrusti . (6)
b) Torque due to drag (τ drag ): The drag moment generated by the i-th propeller acts in the opposite direction of the propeller angular velocity and is expressed in F Pi by k m > 0 is the propeller drag coefficient. The factor (-1) i is used since half of the propellers rotate clockwise and the other half rotates counter-clockwise. This is done in order to have an automatic counterbalance of the drag torques at hovering. The drag torque due to the six propellers expressed in F B is then

T drag i = 0 0 (-1) i k m ω2 i T , (7)
τ drag = 6 i=1 B R Pi T drag i . (8)
Putting together ( 8) and ( 6) in ( 4) we can write

τ = τ (α, β)u, (9)
where τ (α, β) ∈ R 3×6 is the Jacobian matrix that relates the input torque τ to the control input

u = [ω 2 1 ω2 2 ω2 3 ω2 4 ω2 5 ω2 6 ] T ∈ R 6×1 , (10)
i.e., the squares of the rotational speeds of each propeller.

2) Translational dynamics: Thanks to the assumption on the location of the hexarotor and propeller centers of mass, we can express the translational dynamics in F W , using the standard Newton-Euler formulation, as [START_REF] Kendoul | Modeling and control of a small autonomous aircraft having two tilting rotors[END_REF] where f ext represents external disturbances and unmodeled effects, and F (α, β) ∈ R 3×6 is the Jacobian matrix that relates u with the total force produced by the propellers (expressed in body frame), i.e.,

mp = m   0 0 -g   + W R B F (α, β)u + f ext
F (α, β)u = 6 i=1 B R Pi T thrusti . (12)
Notice that in a standard hexarotor α i = β i = 0, for all i = 1 . . . 6. This implies that F (α, β) has rank equal to one (the total force is always directed on the Z B axis).

III. CONTROL DESIGN

The control problem considered here is an output tracking problem. In particular, the hexarotor is tasked to track, simultaneously, a desired trajectory p d (t) with the CoM position p and a given orientation R d (t) with the body orientation W R B . The available control inputs are the six spinning rates of the propellers u defined in [START_REF] Oner | Dynamic model and control of a new quadrotor unmanned aerial vehicle with tilt-wing mechanism[END_REF].

Neglecting the external forces and torques (which are handled by the feedback nature of the control) we rewrite here the hexarotor dynamical model, that is used for the control design

p = 0 0 -g T + 1 m W R B F (α, β)u (13) ωB = -I -1 B (ω B × I B ω B ) + I -1 B τ (α, β)u (14)
W ṘB = W R B [ω B ] ∧ (15)
with [•] ∧ being the usual map from R 3 to so(3).

A. Exact Feedback Linearization and Decoupling Control

In order to apply a feedback linearization technique we rewrite ( 13)-( 14) in a matricial form

p ωB = f + J R J (α, β) u = f + J (α, β)u (16)
where f ∈ R 6 is the drift vector due to the gravity and the rotational inertia, J R =

1 m W R B 0 0 I -1 B ∈ R 6×6 , J (α, β) = F (α,β) τ (α,β)
∈ R 6×6 , and the 6 × 6 matrix J (α, β) will be referred to as the input Jacobian2 . If J (α, β) is invertible we choose the control input as

u = J -1 (α, β) (-f + v) ( 17
)
where v is an additional input, thus obtaining

p ωB = v = v p v R , (18)
i.e. the system is exactly linearized via a static feedback. Fig. 3 shows the control scheme architecture.

In order to obtain an exponential convergence to 0 of the position error p-p d = e p one can choose a linear controller

v p = pd -K p1 ėp -K p2 e p -K p3 t t0 e p , (19)
where the diagonal positive definite gain matrixes K p1 , K p2 , K p3 define Hurwitz polynomials. Now considering the orientation tracking, a popular used parameterization is to resort to Euler angles. However it is well known that they are prone to singularity problems. Keeping this in mind, the controller for the rotational configuration is developed directly on so(3) and thereby it avoids

e R = 1 2 [R T d W R B -W R T B R d ] ∨ , (20)
and the tracking error of the angular velocity e ω ∈ R 3 is given by

e ω = ω B -W R T B R d ω d . (21)
In order to obtain an asymptotic convergence to 0 of the rotational error e R one can choose the following controller

v R = ωd -K R1 e ω -K R2 e R -K R3 t t0 e R ( 22
)
where the diagonal positive definite gain matrixes K R1 , K R2 , K R3 define Hurwitz polynomials also in this case.

B. Discussion on the Invertibility of J (α, β)

The previous control design relies on the invertibility of the matrix J (α, β). To have the 6 DoFs, it is necessary to create acceleration in all direction (i.e. ẍ, ÿ, z, ṗ, q, ṙ). This is possible only if the Jacobian matrix J (α, β) in ( 17) is invertible at all time instances. This implies ρ J = rank(J (α, β)) = rank(J R . J (α, β)) ≡ rank( J (α, β)) = 6, ∀t > 0. Here J R is a nonsingular square matrix as seen in ( 16) and therefore doesn't affect rank(J (α, β)). Therefore J (α, β) is the only rank affecting component and is directly dependent on the propeller tilting configuration α i and β i . Necessary conditions for α i and β i are:

• α i = ±π/2, where rank( J (α, β)) < 6 ∀ β i . Due to the high non linearity of J (α, β) sufficient conditions for the invertibility are hard to find. Figure 4 shows det( J (α, β)) of the presented prototype in Sec.V. All combinations lying in the XY-plane are no feasible solutions.

IV. OPTIMIZATION

The angles α i and β i for i = 1 . . . 6 are adjusted during the pre-flight setup as described in Sec. II. This gives the possibility to change the angles depending on the needs of a particular trajectory. In this section, we consider this capability to optimize the angles α i and β i depending on a predefined desired trajectory to reduce the control effort. As a reminder, the main motive is the full controllability in position and orientation. This comes with the cost of a higher control effort. The objective of this section is therefore to reduce this parasitic effect.

The predominant energy consuming parts of the hexarotor are the propeller motors. Minimizing the control effort through the norm of the control output ||u|| by optimizing the particular α i and β i will as well reduce the energy consumption in flight. To reduce the complexity of the optimization, α i and β i shall be changed in a coordinated way as explained before. We decided to use the same α and β respectively for α i and β i ∈ i = 1 . . . 6, but with different signs for the individual joints. An overview of the compared configurations can be found in table II. The coordinated variation of α i and β i offers two additional advantages: (i) no asymmetries in the hexarotor body and (ii) no or a minimum change of the CoM. Considering these constraints, the optimization problem can be defined as:

min α,β t f 0 ||u||dt (23)
Subject to:

0 < α < π 2 (24) 0 < β < π 2 (25) 0 < ωi , for i ∈ 1..6 (26)
Here ( 24) and ( 25) are defining the lower and upper bounds for α i and β i , while (26) ensures a positive rotation speed ωi for all propellers. The presented minimization problem is a multi-dimensional constrained nonlinear optimization problem and can be solved using the in-build optimization capabilities of MATLAB by exploiting the fmincon-function [START_REF] Cheon | A new PID autotuning strategy with operational optimization for MCFC systems[END_REF].

To compare (minimal control effort) the different configurations shown in table II, we used the presented optimization technique to find the optimal values α and β and the associated ||ω|| min . As trajectory, a typical flight regime has been chosen, which is presented in Section VI-A. The minimum value of the objective function could be found in configuration (b). Therefore all further experiments will be performed by using this configuration:

α = α 1 = -α 2 = α 3 = -α 4 = α 5 = -α 6 and β = β 1 = -β 2 = β 3 = -β 4 = β 5 = -β 6 .
The optimal angles α and β are highly dependent on the desired trajectory. To visualize the influence we conducted a trajectory, where the hexarotor hovers in place but performs a sinusoidal rotation around θ and φ at the same time (see figure 5 (a)). The magnitude of the rotation is increased in 6 steps up to 22.5 • . α and β are increasing accordingly from almost zero values to α = 0.49 rad and β = 0.33 rad for the maximum amplitude.

Figure 6 shows the influence of the optimization itself. For the considered sinusoidal trajectory, we calculated the value of the objective function for a wide variety of α and β . The optimal value is marked by a red circle in Fig. 6.

V. TILTED PROPELLER HEXAROTOR DESIGN

In order to independently control the six DoFs of the UAV, e.g., for tracking a desired position p d (t) and orientation R d (t), the UAV needs to have at least six actuators.

Here we have designed a novel hexarotor architecture with the design concepts discussed in Section II. The six propellers are tilted α i w.r.t. X Pi and β i w.r.t. Y Pi as seen in the Fig. 7. Generally the complexity of such a tilt model would alter the point where the force and torque acts (i.e., the origin of each propeller O Pi ), into a different plane for each propeller P i , thereby making the hexarotor model asymmetric and indirectly increasing the complexity of the system dynamics. A similar problem existed in [START_REF] Voyles | Hexrotor UAV platform enabling dextrous interaction with structurespreliminary work[END_REF]. As shown in Fig. 8, here we have designed the tilting mechanism such that the thrust forces generated by each propeller are in the same point with respect to the main body after tilting. It must be noted here that O Pi (point in which thrust force acts) of each propeller frame F Pi are in the same plane X B Y B as the origin O B of hexarotor body. Each propeller is mounted in an arc frame which is free to rotate in X Pi and Y Pi , so that the tilt angle of α i and β i can be fixed as desired. The radius of the arc(R arc ) is designed equal to the length of the motor (with the propeller attached), so that O Pi always stays at the same location in the X B Y B plane with only its direction vector [X Pi Y Pi Z Pi ] T changing according to the α i and β i orientation. The origin O Pi of each propeller frame is equally spaced with 60 • between each other from the center of the body frame O B to have a symmetric configuration in normal hovering position. It must be noted that the alternating propellers are rotating in clockwise and counter-clockwise rotations, through which the torque produced by each motor shares the same direction with the motor force's in-plane components. The arm in which each propeller set-up is suspended is designed to have a curved architecture with the radius of the curvature, more than the propeller radius (R prop ), so that independent from the tilt of α i and β i , the propellers never come in contact with the arm during flight. The individual α i and β i are fixed following the design constrains discussed in detail in Section II. It is important to note that this hexarotor is made to freely position/orient in free space for 6 DoFs and be able to exert force in the environment to resist any wrench for aerial manipulation without any additional hardware (extra payload).

VI. SIMULATIONS AND RESULTS

Here, we intend to present two simulations performed on the novel tilted propeller hexarotor, designed in Sec. II by applying the controller presented in Sec. III. We aim to prove the important contribution of the papers: (i) Hovering and reorienting the hexarotor, while maintaining a fixed position under the influence of an external force / torque disturbance, (ii) And importantly, we intend to prove the 6 DoFs this hexarotor model achieves with the proposed controller showing the robustness in a more complex trajectory tracking problem. Given the chosen α and β, not all trajectories might be feasible since the negative control outputs u i might occur. This needs to be considered during the pre-trajectory planning step. A. Hovering with external disturbance (f ext /τ ext )

In the first simulation, we tested a hovering trajectory in which the hexarotor maintains a fixed position p but reorients itself changing at the same time the roll φ, pitch θ and yaw ψ angles. This involves hexarotor orienting -12 • w.r.t. X B , 12 • w.r.t. Y B axis and 15 • w.r.t. Z B while still hovering in the position p = [0 0 0] T . Clearly orienting w.r.t. the 3 principal body axes {X B , Y B , Z B } while holding the same position is not feasible in a standard hexarotor UAV. The initial conditions were set to p(t 0 ) = 0, ṗ(t 0 ) = 0, W R B (t 0 ) = I 3 and ω B (t 0 ) = 0. The desired trajectory was chosen as p d (t) = 0 and R d (t) = R X (φ(t))R Y (θ(t))R Z (ψ(t)) with φ(t), θ(t), ψ(t) following the smooth profile having the maximum velocity θmax = 5 • /s and the maximum acceleration θmax = 2.5 • /s 2 . The optimized value of α = 13.6 • and β = 10.6 • obtained from Sec. IV has been used. The gains in Equations ( 19) and ( 22) were set to K p1 = K R1 = 10 I 3 , K p2 = K R2 = 29 I 3 and K p3 = K R3 = 30 I 3 .

Figures 9(a-d) show the result of hovering with external force/torque disturbance. As clearly seen in Fig. 9(c) a constant external force disturbance (f ext = [4 2 1] T N) is applied, along the 3 principal axis {X B , Y B , Z B }, from t = 4 to 9s. Fig. 9(a) shows the position (current (solid line) and desired (dashed line)) brought under control while f ext is applied thanks to the integral term in [START_REF] Ryll | First flight tests for a quadrotor UAV with tilting propellers[END_REF]. Similarly in Fig. 9(d) a constant external torque disturbance (τ ext = [0.2 0.175 0.15] T Nm) is applied, about the 3 principal axes {X B , Y B , Z B }, from t = 12 to 18s. Fig. 9(b) shows the orientation that gets disturbed by this external torque and brought under control after a short transient, thanks to the integral term in (22). The in-zoomed Fig. 9(a) shows that the position tracking error is very minimal in powers of 10 -9 . This simulation provides a first confirmation of the validity of the robustness of the controller during hovering with external disturbance and also the reorientation of the hexarotor while maintaining a fixed position, because of the possibility to create 6 DoFs. This point will also be addressed more thoroughly by the next simulation.

B. 6 DoFs trajectory tracking

In this simulation, we have addressed a more complex trajectory following a square path with vertexes {V 1 , V 2 , V 3 , V 4 , V 5 , V 6 , V 7 }. Each vertex was associated with the following desired positions and orientations3

• V 1 : p d = [0 0 0] T , η d = [0 • 0 • 0 • ] T • V 2 : p d = [2 0 0] T , η d = [-18 • 0 • 0 • ] T • V 3 : p d = [2 3 0] T , η d = [-18 • 12 • 0 • ] T • V 4 : p d = [2 3 1] T , η d = [-18 • 12 • 9 • ] T • V 5 : p d = [2 0 1] T , η d = [-18 • 12 • 0 • ] T • V 6 : p d = [2 0 0] T , η d = [-18 • 0 • 0 • ] T • V 7 : p d = [0 0 0] T , η d = [0 • 0 • 0 • ] T
which were traveled along with rest-to-rest motions with maximum linear/angular velocities of 0.3 m/s and 15 • /s, and maximum linear/angular accelerations of 0.2 m/s 2 and 5 • /s 2 . Figures 10(a-d) show the desired trajectory (p d (t), η d (t)), and the tracking errors (e p (t), e R (t)). The same initial condition as in Section VI-A is considered. The optimized value of α = 26.5 • and β = 19 • for this trajectory obtained from Sec. IV has been used. Here it is clearly illustrated that at the vertex V 4 the hexarotor exploits the 6 DoFs which is one of the main objectives of this titled propeller architecture. Note again how the tracking errors are kept to minimum (in power of 10 -5 ) despite the more complex motion involving several reorientations of the propellers. This confirms again the validity of the controller proposed in the previous Section III.

The interested reader can refer to the video attached to the paper for a more exhaustive illustration of the hovering and 6 DOF hexarotor motion capabilities.

VII. CONCLUSIONS AND FUTURE WORKS In this paper, we have overcome the underactuation, modeling and control issues for a standard UAV with a tilted propeller hexarotor architecture where the propellers can be rotated both w.r.t. X-axis and Y-axis. This concept allows to (i) gain full controllability over the 6 DoFs hexarotor pose in free space, and (ii) optimize the propeller tilt angle with respect to a minimum control effort over a desired trajectory. Carefully analyzing the controllability properties of the dynamic model, resulted in the hexarotor design being combined with trajectory tracking controller based on feedback linearization techniques. A clear validation of the controller's robustness is proved by means of extensive simulations.

Our future works are aimed at (i) developing robust control designs based on different non-linear techniques (e.g., sliding mode, backstepping) for robust stability analysis, (ii) realizing the prototype of the proposed tilted propeller hexarotor in order to experimentally validate the ideas discussed in this paper. Further research is mandatory for feasibility of u. In addition, we are also, (iii) proceeding in the direction of optimizing the tilt angle to get a desired force (F x , F y , F z ∈ R 3 ) and desired torque (T x , T y , T z ∈ R 3 ) along the 3 principal body axes {X B , Y B , Z B }, which is the futuristic objective for any physical UAV interaction with the environment.

Fig. 1 :

 1 Fig. 1: Schematic representation of the hexarotor described in this paper.

1Fig. 2 :

 2 Fig. 2: (a): i-th arm presenting body frame FP i and the generated thrust Tthrust i and drag Tdrag i ; (b) and (c): Visualization of the possible reorientation of the propeller around XP i (case (b)) and YP i (case (c)). The angle of reorientation is denoted with α in (b) and β in (c)

Fig. 3 :

 3 Fig. 3: Control scheme architecture

Fig. 4 :

 4 Fig. 4: Determinant value (z-axis) of J (α, β) of the presented prototype + + + --+ -+ + + + +

Fig. 5 :Fig. 6 :

 56 Fig. 5: (a): Desired sinusoidal trajectories for equal θ and φ. Their amplitude is increased in six steps from 0 • (0 rad) to 22.5 • (≈ 0.39 rad); all other values remain constant (= 0). (b): Optimal values for α and β corresponding to the six trajectories presented in (a)

Fig. 7 :Fig. 8 :

 78 Fig. 7: (a) and (b): Visualization of the possible reorientation of the propeller around XP i (case (a)) and YP i (case (b)). The angle of reorientation is denoted with αi in (a) and βi in (b)

Fig. 9 :

 9 Fig.9: Results of the hovering with external force/torque disturbance. 9(a): Desired (dashed line) and current (solid line) position p d in x(red), y(green) and z(blue). 9(b): Desired (dashed line) and current (solid line) orientation η d in roll(red), pitch(green) and yaw(blue). 9(c-d): external force(f ext ) and torque(τ ext) applied to the hexarotor

Fig. 10 :

 10 Fig. 10: Results of the robust 6 DoFs trajectory tracking. 10(a): Desired (dashed line) and current (solid line) position p d in x(red), y(green) and z(blue). 10(b): Desired (dashed line) and current (solid line) orientation η d in roll(red), pitch(green) and yaw(blue). 10(cd): behavior of the position/orientation tracking errors (ep, eη).

  Center of hexarotor or Center of Mass (CoM) O P i Center of the each propeller group p position of O B (the CoM) in F W Lx i distance between O P i and O B propeller tilt angle about O P i O B β i i-th propeller tilt angle about Y P i ωi i-th propeller spinning velocity about Z P i ω B angular velocity of F B w.r.t. F W expressed in F B τ ext external disturbance torque acting on the hexarotor T thrust ii-th propeller thrust along Z P i T drag i drag due to the i-th propeller along Z P i

	Symbols	Definitions
	m	total mass of the hexarotor
		gravity constant
	O B
	F W	inertial world frame
	F B	hexarotor body frame
	F P i W R B	i-th propeller frame rotation matrix from F B to F W
	B R P i	rotation matrix from F P i to F B
	I B	inertia of the hexarotor frame
	k f	propeller thrust coefficient
	km	propeller drag coefficient
	α i	i-th

and L xi , the O Pi of the 6 propellers lie on the X B Y B plane and are evenly angularly distributed, having a symmetric configuration in normal hovering position.

TABLE I :

 I Main symbols used in the paper

• the CoM and the six propeller centers are coplanar, like in a standard hexarotor; • the propeller centers lie on lines that intersect with the CoM and form 6 equal angles of π/3 each;

In standard hexarotor the input Jacobian J(α, β) has rank equal to four similar to a quadrotor.

Here, for the sake of clarity, we represent orientations by means of the classical roll/pitch/yaw Euler set η ∈ R 3 .

This work has been partially supported by the European Community under Contract ICT 287617 ARCAS.