Denis Efimov

Wilfrid Perruquetti

Anton Shiriaev

On existence of oscillations in hybrid systems

Keywords: Hybrid systems, oscillations, stability/instability

The paper extends the notion of oscillations in the sense of Yakubovich to hybrid dynamics. Several sufficient stability and instability conditions for a forward invariant set are presented. The consideration is illustrated by analysis of a model of two-link compass-gait biped robot.

I. INTRODUCTION

Oscillations constitute one of the main operating modes for many systems in nature or in techniques [START_REF] Fransoise | Oscillations En Biologie: Analyse Qualitative Et Modèles[END_REF], [START_REF] King | Vibrations and Waves[END_REF], [START_REF] Kurths | Synchronization in Oscillatory Networks[END_REF], [START_REF] Burkin | Frequency Methods in Oscillation Theory. Mathematics and Its Applications[END_REF], [START_REF] Rogers | Power System Oscillations[END_REF], [START_REF] Tass | Phase Resetting in Medicine and Biology: Stochastic Modelling and Data Analysis[END_REF]. In some cases it is required to maintain the oscillations, in other cases the oscillations have to be suppressed [START_REF]Control of Oscillations and Chaos: 2000 2nd International Conference Proceedings[END_REF], [START_REF] Fradkov | Introduction to oscillations and chaos[END_REF]. In all cases the conditions of existence of sustained oscillations is of great importance since they allow one to analyze/design a system with desired oscillating properties.

There are many stability theories and definitions of oscillations [START_REF] Fradkov | Introduction to oscillations and chaos[END_REF]. Among them in this work we choose one proposed by Prof. Yakubovich almost 40 years ago [START_REF] Yakubovich | Frequency oscillations conditions in nonlinear systems with stationary single nonlinearity[END_REF]. This approach is rather generic, and it covers periodical and chaotic oscillations.

Contrarily to a pure periodical case, when existence conditions are rather sophisticated [START_REF] Mallet-Paret | The poincaré-bendixson theorem for monotone cyclic feedback systems with delay[END_REF], the conditions of oscillations in the sense of Yakubovich (Y-oscillations) are simple. For Lurie systems, they are formulated in the frequency domain [START_REF] Yakubovich | Frequency oscillations conditions in nonlinear systems with stationary single nonlinearity[END_REF], [START_REF] Yakubovich | Oscillations in systems with discontinuous and hysteresis nonlinearities[END_REF], [START_REF] Yakubovich | Conditions for self-induced oscillations in nonlinear systems[END_REF], [START_REF] Burkin | Frequency Methods in Oscillation Theory. Mathematics and Its Applications[END_REF]. For a generic nonlinear system, the conditions of Y-oscillations are given using Lyapunov arguments [START_REF] Efimov | Oscillatority of nonlinear systems with static feedback[END_REF] or applying the homogeneity framework [START_REF] Efimov | Oscillations conditions in homogenous systems[END_REF]. The main goal of the present paper is to extend the notion of Y-oscillations to hybrid dynamical systems. For this purpose several stability/instability Lyapunov conditions are formulated for hybrid dynamics.

There exist many applications where a system has a hybrid dynamics (continuous and discontinuous) and it is oscillating. The most important one comes from robot locomotion [START_REF] Freidovich | Transverse linearization for an underactuated compass-like biped robot and analysis of the closed-loop system[END_REF], [START_REF] Morris | Hybrid invariant manifolds in systems with impulse effects with application to periodic locomotion in bipedal robots[END_REF]. The phenomenon has a hybrid nature due to impacts occurring when a leg is hitting the ground, and the main operating mode is a periodical oscillation. The design/analysis of robot locomotion [START_REF] Freidovich | Transverse linearization for an underactuated compass-like biped robot and analysis of the closed-loop system[END_REF], [START_REF] Morris | Hybrid invariant manifolds in systems with impulse effects with application to periodic locomotion in bipedal robots[END_REF], [START_REF] Shiriaev | Computing a transverse linearization for mechanical systems with two and more underactuated degrees of freedom[END_REF], [START_REF] Shiriaev | Virtual constraints: a tool for orbital stabilization of nonlinear systems theory[END_REF], [START_REF] Shiriaev | Transverse linearization for hybrid controlled mechanical systems with one passive link[END_REF] or a sliding-mode system [START_REF] Aguilar | Generating self-excited oscillations via two-relay controller[END_REF] as a periodically oscillating system is rather sophisticated. However, relaxing the periodicity requirement and considering Y-oscillations, it is possible to develop more constructive conditions for analysis and design of robot motion.

The problem is introduced and illustrated for a two-link compass-gait biped robot model in Section 2. Some preliminaries are given in Section 3. The main result on conditions of Y-oscillations in hybrid systems is presented in Section 4. Several examples of hybrid oscillating systems are given in Section 5.

The first two authors are with the Non-A project at INRIA -LNE, Parc Scientifique de la Haute Borne, 40 avenue Halley, Bât.A Park Plaza, 59650 Villeneuve d'Ascq, France,{denis.efimov; wilfrid.perruquetti}@inria.fr and with LAGIS (UMR-CNRS 8146), Ecole Centrale de Lille, BP 48, Cité Scientifique, 59651 Villeneuve-d'Ascq, France.

The third author is with the Department of Engineering Cybernetics, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway, e-mail: Anton.Shiriaev@itk.ntnu.no.

II. MOTIVATING EXAMPLE

Dynamics of a two-link compass-gait biped robot with a control torque applied at the hip (schematically shown in Fig. 1) can be described by the hybrid system [START_REF] Freidovich | A passive 2-DOF walker: Hunting for gaits using virtual holonomic constraints[END_REF], [START_REF] Freidovich | Transverse linearization for an underactuated compass-like biped robot and analysis of the closed-loop system[END_REF]:

           ẋ1 = x 2 ; p 1 ẋ2 -p 2 (c 13 ẋ4 + s 13 x 2 4 ) -p 4 sin(x 1 ) = u; ẋ3 = x 4 ; p 3 ẋ4 -p 2 (c 13 ẋ2 -s 13 x 2 2 ) + p 5 sin(x 3 ) = -u, x ∈ C;            x + 1 = x - 3 ; x + 3 = x - 1 ; A x + 2 x + 4 = B x - 2 x - 4 , x ∈ D, (1)
where

x = [x 1 x 2 x 3 x 4 ] T ∈ R 4 is the state vector of the robot, u ∈ R is the control torque; s 13 = sin(x 1 -x 3 ), c 13 = cos(x 1 -x 3 ), A = p 1 -p 2 c -p 3 -p 2 c - -p 2 c - p 3 , B = p 7 c --p 6 -p 6 -p 6 0 , c -= cos(x - 1 -x - 3 ); C = R 4 , D = {x ∈ R 4 : cos(x 1 + ψ) = cos(x 3 + ψ)
} define the sets with continuous and discrete dynamics respectively; p i , i = 1, . . . , 7 are the robot parameters related with its physical counterparts presented in Fig. 1 as

p 1 = (m H + m)l 2 + ma 2 , p 2 = mlb, p 3 = mb 2 , p 4 = (m H l + mb + ml)g, p 5 = mbg, p 6 = mab, p 7 = m H l 2 + 2mal; ψ is
the slope of the walking surface. The standard abbreviations

x -= x(t -) = lim ε→0 x(t -|ε|), x + = x(t + ) = lim ε→0 x(t + |ε|)
are used to denote the values before and after jumps of a solution. The trajectories of (1) with a unique switch at each passage of D correspond to physical behavior of a compass-gait biped robot. For ψ = 2.87π/180 and u = 0 the system has two (known) nontrivial periodic solutions [START_REF] Freidovich | Transverse linearization for an underactuated compass-like biped robot and analysis of the closed-loop system[END_REF]. If the value of ψ is changed or if it is a time-varying sequence (that corresponds to a walking on an irregular surface), then the conditions of existence of periodical solutions is an open question. That is more, it is hard to address this problem with existing theoretical approaches for a time-varying surface slope.

The tools for detection and stabilization of periodical oscillating modes for a constant ψ are largely reported in the literature (see [START_REF] Freidovich | A passive 2-DOF walker: Hunting for gaits using virtual holonomic constraints[END_REF], [START_REF] Freidovich | Transverse linearization for an underactuated compass-like biped robot and analysis of the closed-loop system[END_REF], [START_REF] Morris | Hybrid invariant manifolds in systems with impulse effects with application to periodic locomotion in bipedal robots[END_REF], [START_REF] Shiriaev | Computing a transverse linearization for mechanical systems with two and more underactuated degrees of freedom[END_REF], [START_REF] Shiriaev | Virtual constraints: a tool for orbital stabilization of nonlinear systems theory[END_REF] and references therein). These analysis and design methods are complex and, frequently, only local. The problem complexity is originated by hybrid nature of the dynamics. It is hard to search, analyze and stabilize such a type of behavior. However, if we would skip the requirement that the motion has to be periodical, allowing other types of "oscillating behavior" (that is a natural relaxation for a time-varying ψ), then existence conditions could be more simple and global depending of course on an alternatively defined oscillation concept. Below such a reduction is demonstrated applying Y-oscillation concept.

III. PRELIMINARY RESULTS

This section has three parts. The first one deals with the hybrid system formalism following [START_REF] Cai | Smooth lyapunov functions for hybrid systems part i: Existence is equivalent to robustness[END_REF], [START_REF] Cai | Smooth lyapunov functions for hybrid systems part ii: (pre)asymptotically stable compact sets[END_REF], [START_REF] Ye | Stability theory for hybrid dynamical systems[END_REF]. The second part introduces (pre)asymptotic stability definition and its equivalent Lyapunov characterization from [START_REF] Cai | Smooth lyapunov functions for hybrid systems part ii: (pre)asymptotically stable compact sets[END_REF]. The third part is devoted to sufficient conditions of Lyapunov-type instability for hybrid systems.

In this work, R denotes the real numbers, R + corresponds to the nonnegative real numbers, Z and Z + are stated for integers and nonnegative integers respectively. The symbol | • | denotes an absolute value for a real scalar or Euclidean norm for a vector, |x| A = inf y∈A |y -x| is the distance from a point x ∈ R n to a set A ⊂ R n . A continuous function α : R + → R + belongs to class K if α(0) = 0 and it is strictly increasing, it is from class K ∞ if α ∈ K and it is radially unbounded.

A. Hybrid systems

The considered hybrid system has the following form: These assumptions allow one to define solutions of the system (2) on hybrid time domains [START_REF] Cai | Smooth lyapunov functions for hybrid systems part i: Existence is equivalent to robustness[END_REF], [START_REF] Cai | Smooth lyapunov functions for hybrid systems part ii: (pre)asymptotically stable compact sets[END_REF].

ẋ ∈ F(x) x ∈ C; x + ∈ G(x) x ∈ D, (2)
A subset E ⊂ R + × Z + is a compact hybrid time domain if E = J-1
j=0 ([t j , t j+1 ], j) for some finite sequence of times 0 ≤ t 0 ≤ . . . ≤ t J . It is a hybrid time domain if for all (T, J) ∈ E, E ∩ ([0, T ] × {0, 1, . . . , J}) is a compact hybrid time domain. A hybrid arc is a function φ : E ⇒ O defined on a hybrid time domain E, and such that φ(•, j) is locally absolutely continuous for each j. A hybrid arc φ : dom φ ⇒ O is a solution (trajectory) of (2) if φ(0, 0) ∈ C ∪ D and:

(S1) φ(t, j) ∈ F[φ(t, j)] for all j ∈ Z + and almost all t ∈ R + such that (t, j) ∈ dom φ and φ(t, j) ∈ C;

(S2) φ(t, j + 1) ∈ G[φ(t, j)] for all (t, j + 1) ∈ dom φ and φ(t, j) ∈ D.

The conditions C ∪ D = O guarantees existence of solutions for all initial conditions. A solution φ to the hybrid system is called nontrivial if domφ contains at least one point different from (0, 0); it is called maximal if its domain cannot be extended. A solution is called complete if its domain is unbounded (complete solutions are maximal), let S(x) denote the set of all maximal solutions for (2) starting from a point x ∈ O. The hybrid system (2) is said to be forward complete on O if for all x ∈ O, each φ ∈ S(x) is complete. The hybrid system (2) is said to be forward invariant on a set A ⊂ O if, for all x ∈ A, each φ ∈ S(x) is such that φ(t, j) ∈ A for all (t, j) ∈ domφ (this definition of forward invariance does not need existence of trajectories for all points of A). A point x 0 is called an equilibrium of (2) if it is forward invariant for the system.

If G(D)∩D = ∅, then multiple consequent switches or a Zeno-behavior are possible. In many cases (if D ⊂ C, for example) such a behavior does not correspond to physical nature of the considered system, but it is hard to avoid it in the formulation Further, talking about a solution φ of the system (2) we will understand a maximal solution φ ∈ S(x) for some x ∈ O, except the restrictionφ ∈ Ŝ(x) is explicitly stated.

B. Stability of hybrid systems

A continuous function ω : O → R + is proper on O if ω(x i ) → +∞ when x i converges to a boundary of O or |x i | → +∞. For a compact set A ⊂ O, the function ω is a proper indicator for A on O if it is proper and A = {x ∈ O : ω(x) = 0}. Definition 1. [23] Let A ⊂ O be compact.
1) A is prestable for the system (2) if for each ε > 0, there exists δ > 0 such that any solution φ with |φ(0, 0)

| A ≤ δ satisfies |φ(t, j)| A ≤ ε for all (t, j) ∈ dom φ.
2) A is preattractive for the system (2) if there exists δ > 0 such that any solution φ with |φ(0, 0)| A ≤ δ is bounded with respect to O, and if it is complete, then lim t+j→+∞ |φ(t, j)| A = 0.

3) A is preasymptotically stable if it is prestable and preattractive.

4)

A is asymptotically stable if it is preasymptotically stable and there exists δ > 0 such that any maximal solution φ with

|φ(0, 0)| A ≤ δ is complete.
The set of all x ∈ O from which all solutions are bounded with respect to O and the complete ones converge to A is called the basin of preattraction and denoted as O a A . If the conditions of Definition 1 are valid for any δ ∈ R + , then the set A is globally preattractive/preasymptotically or asymptotically stable for (2). Definition 2. A compact set A ⊂ O is called Non-Zeno (NZ) prestable (preattractive/preasymptotically stable/asymptotically stable) for the system (2) if the corresponding properties of Definition 1 are satisfied for all φ ∈ Ŝ only. If for the system (2) the set A is (pre)asymptotically stable, then it is also NZ (pre)asymptotically stable, the converse is in general false.

Definition 3. Let X ⊂ O be open and ω : X → R + be continuous. A function V : X → R + is said to be a smooth Lyapunov function for (X , F, G, C, D, ω) if there exist α 1 , α 2 ∈ K ∞ such that α 1 (ω(x)) ≤ V (x) ≤ α 2 (ω(x)) ∀x ∈ X .
In this work we will use the following designations:

Sup DV (x) = sup f ∈F(x) ∇V (x), f , Inf DV (x) = inf f ∈F(x) ∇V (x), f ; Sup V + (x) = sup g∈G(x)∩X V (g), Inf V + (x) = inf g∈G(x)∩X V (g).
This definition does not impose a restriction on derivative or increment of V since both stability and instability properties will be studied below.

Theorem 1. [23]

Let A be compact and X be open, A ⊂ X ⊂ O, and let ω be a proper indicator for A on X . The following statements are equivalent for (2):

1) There exists a smooth Lyapunov function V for (X , F, G, C, D, ω) such that

Sup DV (x) ≤ -V (x) ∀x ∈ X ∩ C; (3)
Sup V + (x) ≤ e -1 V (x) ∀x ∈ X ∩ D and G(D ∩ X ) ⊂ X .
2) The set A is preasymptotically stable, its basin of preattraction O a A contains X , and X is forward invariant for (2).

Remark 1. Note that the sufficient part of this theorem holds if the condition (3) is replaced with the following one:

Sup DV (x) ≤ -ρ[ω(x)] ∀x ∈ X ∩ C; (4) Sup V + (x) ≤ λV (x) ∀x ∈ X ∩ D,
for some function ρ ∈ K and a constant λ ∈ [0, 1). Indeed, in this case V [φ(t, j)] ≤ σ 1 (V [φ(t j , j)])σ 2 (e tj -t ) for all j ∈ Z + and t ∈ [t j , t j+1 ] such that (t, j) ∈ dom φ and φ(t, j) ∈ C for some functions σ 1 , σ 2 ∈ K; and V [φ(t, j + 1)] ≤ λV [φ(t, j)] for all (t, j + 1) ∈ dom φ and φ(t, j) ∈ D. Therefore, according to Definition 1 the system is prestable and preattractive with respect to the set A (A is preasymptotically stable).

Remark 2. For the sufficient part of Theorem 1 it is enough to have a continuously differentiable or a Lipschitz continuous function V (the same for the results presented below). For a Lipschitz continuous function V the Dini derivatives have to be used.

Following [START_REF] Ye | Stability theory for hybrid dynamical systems[END_REF], a more simple for verification sufficient condition can be proposed for NZ preasymptotical stability of A based on the following assumption.

Assumption 1. For all φ ∈ Ŝ, sup (t,j)∈dom φ t = +∞ iff sup (t,j)∈dom φ j = +∞.

These conditions imply that the hybrid time domains of complete solutions in Ŝ are unbounded with respect to both t and j. This assumption excludes from consideration a case with a finite number of jumps and the Zeno-behavior. The reasons for this condition introduction are the following. If the number of switches is finite for a complete solution φ ∈ S, then the asymptotic transients are predefined by the continuous dynamics only. The Zeno-behavior is rather usual in hybrid models, but frequently it has no physical meaning (like in the compass-gait biped robot or the bouncing ball [START_REF] Cai | Smooth lyapunov functions for hybrid systems part ii: (pre)asymptotically stable compact sets[END_REF], for instance) since multiple impacts are possible in an equilibrium only.

Theorem 2. Let A be compact and X be open, A ⊂ X ⊂ O, and let ω be a proper indicator for A on X . Let Assumption 1 be satisfied and there exist a continuously differentiable Lyapunov function V for (X , F, G, C, D, ω) such that

Sup DV (x) ≤ 0 ∀x ∈ X ∩ C, (5)
G(D ∩ X ) ⊂ X and for any solution φ ∈ Ŝ for all (t j+1 , j + 1) ∈ dom φ such that φ(t j , j) ∈ X and φ(t j+1 , j + 1) ∈ X for

some α ∈ K V [φ(t j+1 , j + 1)] -V [φ(t j , j)] t j+1 -t j ≤ -α(V [φ(t j+1 , j + 1)]). (6)
Then the set A is NZ preasymptotically stable, its basin of preattraction O a A contains X and X is forward invariant for φ ∈ Ŝ of (2).

Proof: Note that V [φ(t, j)] ≤ V [φ(t j , j)] for all j ∈ Z + and t ∈ [t j , t j+1 ] such that (t, j) ∈ dom φ and φ(t, j) ∈ C (this relation follows from [START_REF] Rogers | Power System Oscillations[END_REF]). According to [START_REF] Tass | Phase Resetting in Medicine and Biology: Stochastic Modelling and Data Analysis[END_REF], in addition, V [φ(t j+1 , j + 1)] ≤ V [φ(t j , j)] for all (t j+1 , j + 1) ∈ dom φ.

Combining these facts, since V [φ(t j+1 , j)] ≤ V [φ(t j , j)], we have that V [φ(t, j)] ≤ V [φ(0, 0)] for all (t, j) ∈ dom φ, φ ∈ Ŝ
and the set A is NZ prestable for the system (2) and it is forward complete. Define V j = V [φ(t j , j)] to simplify the notation, then

∆V j = V j+1 -V j ≤ -(t j+1 -t j )α(V j+1 ) ≤ -(t j+1 -t j )α(V N +1 )
for all (t N , N ) ∈ dom φ and j ≤ N (as it was noted above the sequence V j is not increasing). In addition, t j+1 -t j > 0 since φ ∈ Ŝ. Thus

V N +1 -V 0 = N j=0 ∆V j ≤ -(t N +1 -t 0 )α(V N +1 ), and
α(V N +1 ) ≤ V 0 -V N +1 t N +1 -t 0 ≤ V 0 t N +1 -t 0 .
Under Assumption 1 there is an infinite number of switches and t N +1 → +∞, consequently lim t+j→+∞ V [φ(t, j)] = 0 or equivalently (V is a Lyapunov function) lim t+j→+∞ |φ(t, j)| A = 0, therefore A is NZ preattractive for the system (2).

Despite of the conditions of Theorem 2 include the inequality (6) in the time domain for φ ∈ Ŝ, in some cases it is more simple to apply Theorem 2 than Theorem 1 (that we will show later on the examples). Note that without the condition

Sup V + (x) ≤ V (x) ∀x ∈ X ∩ D the function V decreasing may be non-monotone, i.e. V [φ(t j+1 , j + 1)] > V [φ(t j+1 , j)]
for some (t j+1 , j + 1) ∈ domφ and φ(t j+1 , j) ∈ D. In another case V may be decreasing on jumps in D, but not necessarily strictly decreasing for the continuous flow in C.

C. Instability of hybrid systems

This section contains a definition of instability (similarly to [START_REF] Forni | Instability and overshoots of solutions for a class of homogeneous hybrid systems by lyapunov-like analysis[END_REF], [START_REF] Ye | Stability theory for hybrid dynamical systems[END_REF]) and the corresponding sufficient Lyapunov-type conditions.

Definition 4. Let U ⊂ O be compact.

1) U is preunstable for the system (2) if it is not prestable in the sense of Definition 1, i.e. there exists ε > 0 such that for each δ > 0 there is a solution φ with |φ(0, 0)| U ≤ δ that satisfies |φ(t, j)| U > ε for some (t, j) ∈ domφ.

2) U is prerepulsive for the system (2) if there is some δ > 0 such that for any 0 < δ ≤ δ there exists a solution φ with |φ(0, 0)| U ≤ δ for which there are some (t , j ) ∈ dom φ such that the property |φ(t, j)| U > δ holds for all (t, j) ∈ dom φ with t ≥ t , j ≥ j .

3) U is strictly prerepulsive for the system (2) if for some δ > 0 and any solution φ with |φ(0, 0)| U ≤ δ there exist some (t , j ) ∈ dom φ such that the property |φ(t, j)| U > δ holds for all (t, j) ∈ dom φ with t ≥ t , j ≥ j .

As in the continuous case [START_REF] Hahn | Stability of Motion[END_REF], the hybrid system (2) can be preattractive with respect to a set A and preunstable with respect to this set simultaneously. If the set U is prerepulsive for the system (2), then it is preunstable, but this set still can be attractive for a some subset of initial states in a vicinity of U. If it is strictly prerepulsive, then all existent solutions exit a neighborhood of U. The domain of prerepulsion O r U contains x ∈ O such that φ(0, 0) = x and |φ(t, j)| U > δ for all (t, j) ∈ domφ with t ≥ t , j ≥ j ((t , j ) ∈ dom φ). Definition 5. A compact set U ⊂ O is called NZ preunstable (prerepulsive/strictly prerepulsive) for the system (2) if the corresponding properties of Definition 4 are satisfied for φ ∈ Ŝ only.

If U is NZ preunstable/prerepulsive for (2), then it is also preunstable/prerepulsive. If U is strictly prerepulsive, then it is also NZ strictly prerepulsive. Theorem 3. Let U ⊂ X be compact and X ⊂ O be open, G(D ∩ X ) ⊂ X . Let a continuous function ω : X → R + be a proper indicator for U on X . Let also each φ ∈ S with φ(0, 0) ∈ X be nontrivial.

1)

The set U is strictly prerepulsive for the system (2), its domain of prerepulsion O r U contains X , if there exists a continuously differentiable Lyapunov function U : X → R + for (X , F, G, C, D, ω) such that for some function ρ ∈ K and a constant λ ∈ [1, +∞)

Inf DU (x) ≥ ρ[ω(x)] ∀x ∈ X ∩ C; (7)
Inf U + (x) ≥ λU (x) ∀x ∈ X ∩ D.
2) Let Assumption 1 be satisfied. The set U is NZ strictly prerepulsive for the system (2), its domain of prerepulsion O r U contains X for φ ∈ Ŝ, if there exists a continuously differentiable Lyapunov function U : X → R + for (X , F, G, C, D, ω)

such that for any solution φ ∈ Ŝ for all (t j+1 , j + 1) ∈ dom φ such that φ(t j , j) ∈ X and φ(t j+1 , j + 1) ∈ X for some

α ∈ K U [φ(t j+1 , j + 1)] -U [φ(t j , j)] t j+1 -t j ≥ α(U [φ(t j , j)]) (8)
and

Inf DU (x) ≥ 0 ∀x ∈ X ∩ C. (9)
Proof: In the first case U [φ(t, j)] ≥ σ 1 (U [φ(t j , j)])σ 2 (e t-tj ) for all j ∈ Z + and t ∈ [t j , t j+1 ] such that (t, j) ∈ dom φ and φ(t, j) ∈ C for some functions σ 1 , σ 2 ∈ K; and U [φ(t, j + 1)] ≥ λU [φ(t, j)] for all (t, j + 1) ∈ dom φ and φ(t, j) ∈ D.

Therefore, U [φ(t, j)] is a strictly increasing function of (t, j) ∈ dom φ while φ(t, j) ∈ X . By Definition 4 the system (2) is preunstable and prerepulsive with respect to the set U. The domain of prerepulsion O r U contains X . In the second case, from [START_REF] Yakubovich | Frequency oscillations conditions in nonlinear systems with stationary single nonlinearity[END_REF], U [φ(t, j)] ≥ U [φ(t j , j)] for all j ∈ Z + and t ∈ [t j , t j+1 ] such that (t, j) ∈ dom φ and φ(t, j) ∈ C, and U [φ(t j+1 , j + 1)] ≥ U [φ(t j , j)] for all (t j+1 , j + 1) ∈ dom φ in accordance with [START_REF] Fradkov | Introduction to oscillations and chaos[END_REF]. Consequently,

∆U j = U j+1 -U j ≥ [t j+1 -t j ]α(U j ) ≥ [t j+1 -t j ]α(U 0 )
for all (t j+1 , j + 1) ∈ dom φ since the sequence U j = U [φ(t j , j)] is not decreasing (t j+1 > t j since φ ∈ Ŝ). Finally we obtain

U j+1 -U 0 = j i=0 ∆U i ≥ [t j+1 -t 0 ]α(U 0 ) and U j+1 ≥ U 0 + [t j+1 -t 0 ]α(U 0 ).
Since under Assumption 1 for a complete solution φ ∈ Ŝ, we have t j+1 → +∞ with j → +∞, then the sequence U j is unbounded. If the solution φ is not complete, then the trajectory escapes to infinity in a finite time, hence U is again unbounded.

A solution φ may also leave X in a finite time. Then by Definition 4 the system ( 2) is preunstable and prerepulsive with respect to the set U (again X ⊂ O r U ). Note that the condition ( 8) is sufficient to prove that the set U is preunstable for the system (2). Contrarily the condition (7), the conditions ( 8), ( 9) do not restrict the Lyapunov function U behavior after switches (it may decrease on D).

IV. OSCILLATION CONDITIONS

Extending the results from [START_REF] Efimov | Oscillatority of nonlinear systems with static feedback[END_REF], [START_REF] Yakubovich | Frequency oscillations conditions in nonlinear systems with stationary single nonlinearity[END_REF], the definition of Y-oscillations for hybrid systems can be as follows (recall, a function g : R n → R is called monotone if the condition x 1 ≤ x 1 , . . . , x n ≤ x n implies that everywhere either g(x 1 , . . . , x n ) ≤ g(x 1 , . . . , x n ) or g(x 1 , . . . , x n ) ≥ g(x 1 , . . . , x n )). Definition 6. [START_REF] Yakubovich | Frequency oscillations conditions in nonlinear systems with stationary single nonlinearity[END_REF], [START_REF] Efimov | Oscillatority of nonlinear systems with static feedback[END_REF] A complete solution φ ∈ S(x) with x ∈ O of the system (2) is called [π -, π + ]-oscillation with respect to the output ψ = η(x) (where η :

O → R is a continuous monotone function) if lim t+j→+∞ η[φ(t, j)] = π -; lim t+j→+∞ η[φ(t, j)] = π + ; -∞ < π -< π + < +∞.
A complete solution φ ∈ S(x) with x ∈ O of the system (2) is called oscillating, if there exist some output ψ and constants π -, π + such that φ is [π -, π + ]-oscillation with respect to the output ψ. A forward complete system (2) is called oscillatory, if for almost all x ∈ O the solutions φ ∈ S(x) of the system are oscillating. The oscillatory system (2) is called uniformly oscillatory, if for almost all x ∈ O for corresponding solutions φ ∈ S(x) there exist the output ψ and the constants π -, π + not depending on initial conditions.

In other words the solution φ is oscillating if the output ψ(t, j) = η[φ(t, j)] is asymptotically bounded and there is no single limit value of ψ(t, j) for t + j → +∞. The term "almost all solutions" is used to emphasize that generally the system (2) has a nonempty set of equilibrium points, thus there exists a set of initial conditions with zero measure such that the corresponding solutions are not oscillating. The notion of oscillations in the sense of Yakubovich is rather generic including periodical oscillations (limit cycles), quasi-periodical, recurrent and chaotic trajectories. The oscillating trajectories could be repelling being oscillating. The trajectories could also be unbounded, it is required to find a function of the state vector η, that is bounded and admits certain requirements introduced in Definition 6.

Theorem 4. Let A ⊂ X A and U ⊂ X U ⊂ A be compact and X = X A ∪ X U ⊂ O be open. Let continuous functions ω A : X → R + , ω U : X → R + be proper indicators for A and U on X respectively. Let G(D ∩ X ) ⊂ X , each φ ∈ S with φ(0, 0) ∈ X be nontrivial and the system (2) have two continuously differentiable Lyapunov functions V and U such that for all x ∈ X :

υ 1 [ω A (x)] ≤ V (x) ≤ υ 2 [ω A (x)], υ 3 [ω U (x)] ≤ U (x) ≤ υ 4 [ω U (x)], υ 1 , υ 2 , υ 3 , υ 4 ∈ K ∞ and 1
) the conditions ( 4) or ( 5), ( 6) with Assumption 1 are valid for the function V on X A ;

2) the conditions ( 7) or ( 8), ( 9) with Assumption 1 are valid for the function U on X U ;

3) Ω ∩ Ξ = ∅, where Ω = A ∩ X U and Ξ ⊂ R n is a set with zero Lebesgue measure containing all equilibriums of the system.

Then the system (2) has an oscillating solution.

Proof: According to the results of theorems 1, 2 and 3, the set A is NZ preasymptotically stable (its basin of preattraction O a A contains X A , and X A is forward invariant for (2)) and the set U is NZ strictly prerepulsive for the system (2), its domain of prerepulsion O r U contains X U . Thus in X A all solutions φ ∈ Ŝ(x) approach the set A, while inside this set there is a repulsing subset U with the repulsion region X U , then all φ ∈ Ŝ(x) have to approach the intersection Ω = A ∩ X U , a forward invariant compact set for the system (2) for φ ∈ Ŝ(x). The set Ω does not contain the system equilibriums by the theorem conditions. Consider a solution φ ∈ Ŝ(x) with x ∈ Ω of the system (2) (it is complete by the consideration above). Then there exists ψ = η(x) = x i for some 1 ≤ i ≤ n (the i th coordinate of the system) such that the solution φ is oscillating with respect to this output ψ. Assume it is not, it implies that lim t+j→+∞ φ i (t, j) = lim t+j→+∞ φ i (t, j) for all 1 ≤ i ≤ n and that φ converges to an equilibrium point. Since there is no equilibriums in Ω, there is such an index i with -∞ < lim t+j→+∞ φ i (t, j) < lim t+j→+∞ φ i (t, j) < +∞ and the solution φ is oscillating. Remark 3. If the conditions (4) and ( 7) are satisfied, then the system (2) is oscillating.

Remark 4. The obtained in Theorem 4 conditions extend the results of the work [START_REF] Efimov | On periodical oscillations of lurie systems with discontinuous nonlinearity[END_REF], obtained for a linear system with a sliding-mode feedback (in that case the solutions are absolutely continuous), on a more general case.

In the next section we will show several examples of this theorem application.

V. EXAMPLES

In this section we will consider three examples, in all cases D ⊂ C.

A. A planar hybrid system

Consider the following planar hybrid system:

F(x) = (1 -|x|)x 1 + x 2 (1 -|x|)x 2 -x 1 , G(x) = x 2 x 1 , (10)
where

x = [x 1 x 2 ] T , O = R 2 , D = {x ∈ O : x 1 + kx 2 = 0}, C = O and k ∈ R \ {-1, 1} is a parameter.
The system (10) has the single equilibrium at the origin for k = ±1 (Ξ = {0}). Next, it is easy to verify that the Lyapunov function W = |x| 2 for the continuous dynamics F has the derivative

Ẇ = 2(1 - √ W )W for all x ∈ R 2 , then W (t, t 0 , W 0 ) = [1 -(1 -W -0.5 0
)e t0-t ] -2 for any initial condition W (t 0 ) = W 0 , which according to [START_REF] Efimov | Oscillatority of nonlinear systems with static feedback[END_REF], [START_REF] Efimov | Oscillations conditions in homogenous systems[END_REF] implies existence of oscillating trajectories into the set W (x) = 1 (for ẋ = F(x)) and forward completeness of [START_REF] Mallet-Paret | The poincaré-bendixson theorem for monotone cyclic feedback systems with delay[END_REF]. To prove that for any initial condition x(t , j) ∈ O, t ∈ R + , j ∈ Z + for solutions of the system (10) always there exists a time instant t ≥ t such that x(t , j) ∈ D assume that it is not, then from the expression of W (t, t 0 , W 0 ) given above the trajectories of the continuous dynamics converge to the set where W (x) = 1 and the continuous dynamics is reduced to a linear oscillator, thus t exists and Assumption 1 is satisfied.

For the chosen G we have that W (x + ) = W (x) for all x ∈ D, therefore the value of W is not changing for the discrete dynamics. In addition, for k = ±1 the set D \ {0} is not invariant for [START_REF] Mallet-Paret | The poincaré-bendixson theorem for monotone cyclic feedback systems with delay[END_REF], therefore for almost all instants of time the behavior of W is predefined by the continuous dynamics.

To apply Theorem 4 we can choose the origin as the set U, U (x) = W (x) and

X U = {x ∈ O : |x| < 1}, then U (t, j)
is a strictly increasing function on X U and for any solution φ ∈ Ŝ of (10) with φ(0, 0) ∈ X U we have U [φ(t j+1 , j + 1)] -U [φ(t j , j)] ≥ (t j+1 -t j )2(1 -U [φ(t j , j)])U [φ(t j , j)] for all (t j , j) ∈ dom φ while φ rests in X U , i.e. the conditions of Theorem 3 are satisfied.

Let us select the set A = {x ∈ O : |x| ≤ 1}, then X A = O \ A. Using the function W we can prove that this A is NZ preasymptotically stable for [START_REF] Mallet-Paret | The poincaré-bendixson theorem for monotone cyclic feedback systems with delay[END_REF]. Indeed, for any solution φ ∈ Ŝ of (10) with φ(0, 0) ∈ X A we have W [φ(t j+1 , j + 1)] -W [φ(t j , j)] ≤ (t j+1 -t j )2(1 -W [φ(t j+1 , j + 1)])W [φ(t j+1 , j + 1)] for all (t j , j) ∈ dom φ while φ rests in X A . Select

V (x) = max{0, W (x)-1}, then outside of the set A we have V [φ(t j+1 , j+1)]-V [φ(t j , j)] ≤ -(t j+1 -t j )α{V [φ(t j+1 , j+1)]} for α(s) = -(1 - √ 1 + s)(1 + s) ∈ K
and all conditions of Theorem 2 are verified.

The set Ω = {x ∈ O : |x| = 1} does not contain the equilibrium of [START_REF] Mallet-Paret | The poincaré-bendixson theorem for monotone cyclic feedback systems with delay[END_REF]. Therefore, all conditions of Theorem 4 are satisfied and the system [START_REF] Mallet-Paret | The poincaré-bendixson theorem for monotone cyclic feedback systems with delay[END_REF] has an oscillating solution.

The system trajectories for different values of k are shown in Fig. 2. As we can deduce from these results of simulation, the system stay oscillating for all values of k, but the types of the oscillating trajectories differ significantly (there exist two disjoint oscillating subsets in Ω for k = |1.5| and a discontinuous limit cycle for k = |0.5|).

2 - 1 - 0 1 2 2 - 1 - 0 1 2 2 - 1 - 0 1 2 2 - 1 - 0 1 2 2 - 1 - 0 1 2 2 - 1 - 0 1 2 2 - 1 - 0 1 2 2 -

B. A third order hybrid system

Another example is a third order hybrid system: The last example is the motivating two-link compass-gait biped robot. In this case the function G maps the set D to itself, thus formally the system possesses a Zeno-behavior. However, this type of trajectories do not correspond to physical processes in a biped robot, that is why they are neglected and only Ŝ is considered. We will be interested in the passive walking, i.e.

F(x) = A 1 x if |x| < 1 A 2 x if |x| ≥ 1 , G(x) = ϕ(|x|)    x 2 x 1 x 3    , (11)
A 1 =    1 1 0 -1 1 1 0 -1 1    , A 2 =    0 1 0 -1 0 1 0 -1 0    ,
when u = 0. In this case the energy E(x) = 0.5p 1 x 2 2 -p 2 cos(x 1 -x 3 )x 2 x 4 + 0.5p 3 x 2 4 + p 4 (cos(x 1 ) -1) + p 5 (1 -cos(x 3 ))

is the system invariant while x ∈ C (during the continuous transition phase the system is passive). In the discontinuous phase for x ∈ D the energy E can change its value since during the impact the system momentum is preserved only. Moreover, in the vicinity of the stable limit cycle the system can lose or gain the energy during the impact. For the chosen parameters of the system (1) the energy is conserved on the limit cycle.

For this example we will assume that the jump set is defined as

D = {x ∈ R 4 : cos(x 1 + ψ + u ψ ) = cos(x 3 + ψ + u ψ )},
where u ψ is the slope control from the robot side, i.e. the robot has adjustable ends that can be used to locally advance or postpone the impact (that is equivalent to small variation of the slope angle by a foot regulation). This control is acting during the impacts only by changing the jump map properties. The idea is to use the energy E as a Lyapunov function to determine the stability/instability of the robot. Let E ψ corresponds to the system energy on the limit cycle, then ensuring the energy increasing for E ≤ E ψ (the system instability) or its decreasing for E ≥ E ψ (stability) by the control u ψ choice we can guarantee the system oscillations (for E = E ψ = 0 the corresponding trajectories in C are oscillating). It is worth to stress that such a regulation of energy serves as a good example of Theorem 4 application for the system (1), however, the energy stabilization does not imply the robot walking. As we will show in simulation, at least locally this approach leads to the robot locomotion.

By its construction the control u ψ is rather small, thus for its design we can use the first approximation of the system energy after an impact. Let φ(t, j) be a trajectory of the compass-gait biped robot from Ŝ, then E(φ(t j+1 , j + 1)) = E[G(φ(t j+1 , j))]

and since the map G depends on ψ + u f let us define Σ(φ(t j+1 , j), ψ

+ u ψ ) = E[G(φ(t j+1 , j))], then E(φ(t j+1 , j + 1)) = Σ(φ(t j+1 , j), ψ + u ψ ) = Σ(φ(t j+1 , j), ψ) + ∇ ψ Σ(φ(t j+1 , j), ψ)u ψ + o[u 2 ψ ],
where o(•) is the Landau notation (for sufficiently small amplitudes of u ψ the term o[u 2 ψ ] can be neglected). It is assumed that the state vector is available for measurements, then

u ψ (x) = -ϕ 1 [∇ ψ Σ(x, ψ)]ϕ 2 [Σ(x, ψ) -E ψ ] (12)
is an admissible control strategy with ϕ i (s)s ≥ 0, i = 1, 2 for any s ∈ R. Indeed, introducing the error (t, j) = E(φ(t, j)) -E ψ with the substitution of the control [START_REF] Yakubovich | Conditions for self-induced oscillations in nonlinear systems[END_REF] into the first approximation we get at the impact:

(t j+1 , j + 1) = (t j+1 , j) + ∇ ψ Σ(φ(t j+1 , j), ψ)u ψ (φ(t j+1 , j))

= (t j+1 , j) -∇ ψ Σ(φ(t j+1 , j), ψ)ϕ 1 [∇ ψ Σ(φ(t j+1 , j + 1), ψ)]ϕ 2 ( (t j+1 , j)), and
| (t j+1 , j + 1)| < | (t j+1 , j)| for a properly chosen ϕ i , i = 1, 2.
The following fact can be formulated.

Corollary 1. For the system (1) with the control [START_REF] Yakubovich | Conditions for self-induced oscillations in nonlinear systems[END_REF], there exist ε > 0, ψ > 0, E ψ > 0 and ϕ i , i = 1, 2 such that into the forward invariant set X = {x ∈ R 4 : |E(x) -E ψ | ≤ ε} there is an oscillating trajectory.

Proof: Note that (t, j) = (t j , j) for all t ∈ [t j , t j+j ] with (t j+1 , j + 1) ∈ dom φ, while for a properly selected ϕ i , i = 1, 2 the property | (t j+1 , j + 1)| < | (t j+1 , j)| is satisfied for φ(0, 0) ∈ X with a corresponding ε > 0 (the first order approximation arguments). Therefore, X is a forward invariant set for (1), [START_REF] Yakubovich | Conditions for self-induced oscillations in nonlinear systems[END_REF]. Next, there are ψ > 0 and E ψ > 0 such that using Theorem 4 it is possible to show that into the set {x ∈ R 4 : E(x) = E ψ } there is an oscillating trajectory.

The results of the system simulation with ϕ 1 (s) = sign(s) and ϕ 2 (s) = 0.1 tanh(0.1s) are given in Fig. 4. The initial conditions does not belong to the set of attraction of the limit cycle, and the red dash line at the top of the figure represents the system trajectory without the control. The controlled trajectory is shown by the blue solid line, it converges to the limit cycle. The energy E behavior (it converges to the desired level E ψ ) and the control u ψ are given at the bottom of the figure.

The control deviates the slope value during the impact in order to regulate the energy level. During the continuous dynamics phase, the system performs the passive walking (the E is piecewise constant). As we can conclude from the given results of simulation, the proposed control increases the basin of preattraction of the limit cycle and improves the robustness abilities of the robot with respect to the slope variation (the disturbance δ ψ (t)). The paper develops the framework of oscillations in the sense of Yakubovich to the class of hybrid systems (having continuous and discrete dynamics). Some sufficient Lyapunov conditions for stability/instability of hybrid systems are presented. In order to facilitate analysis of physical/mechanical systems described by hybrid models without Zeno-behavior, a part of results is formulated for a subset of maximal solutions of the hybrid system. The utility of the arguments is demonstrated on several examples. For the compass-gait biped robot a control strategy is established, which locally stabilizes an unstable limit cycle and ensures locomotion for a time-varying slope of the walking surface. Extensions to different robotic systems are planned in future works.

Figure 1 .

 1 Figure 1. Schematic of the compass-gait biped on a shallow slope

  where x ∈ O ⊂ R n is the state vector (it may contain continuously changing and discrete or logic-based variables); C ⊂ O and D ⊂ O are the sets where flow or jumps can occur respectively; the (set-valued) mapsF : O ⇒ R n and G : O ⇒ R ndetermine the continuous and discrete dynamics of the system (2). Following[START_REF] Cai | Smooth lyapunov functions for hybrid systems part ii: (pre)asymptotically stable compact sets[END_REF] we impose the Standing Assumptions (SAs) on (2): (SA1) the sets C ⊂ O and D ⊂ O are relatively closed in O (i.e. there exist some closed sets Y C , Y D such that D = Y D ∩O and C = Y C ∩ O); (SA2) the map F : O ⇒ R n is outer semicontinuous 1 and locally bounded, F(x) is nonempty and convex for any x ∈ C; (SA3) the map G : O ⇒ R n is outer semicontinuous and locally bounded, G(x) is nonempty subset of O for any x ∈ D.

( 2 )

 2 since the sets C and D have to be closed. In addition, if D ⊂ C then on the intersection C ∩ D multiple solutions can be initiated, while frequently in practice only a jump is possible in D. Let us define the set of maximal solutions initiated from a point x ∈ O without multiple consequent switches with jumps in D as Ŝ(x), i.e. Ŝ(x) = {φ ∈ S(x) : t j = t j+1 , t j+1 = arg inf t>tj φ(t, j) ∈ D \ Ξ, (t j+1 , j + 1) ∈ dom φ} where the set Ξ ⊂ O contains all equilibriums of (2). This set can be empty or, if G(D) ∩ D = ∅, then Ŝ(x) = S(x).

Figure 2 .

 2 Figure 2. The hybrid system (10) trajectories (x 1 ,x 2 ) for k ∈ {-0.5, 0.5, -1.5, 1.5}

Figure 3 .

 3 Figure 3. The hybrid system (11) trajectories

Figure 4 .

 4 Figure 4. The passive walking with regulated impacts

Figure 5 .

 5 Figure 5. The passive walking with time-varying slope

A set-valued map F : O ⇒ R n is outer semicontinuous at x ∈ O if lim sup x→x F(x) ⊂ F(x), the map F : O ⇒ R n iscalled outer semicontinuous if it has this property for all x ∈ O.

The set D is not invariant for the continuous dynamics for |k| = 1, the system [START_REF] Yakubovich | Oscillations in systems with discontinuous and hysteresis nonlinearities[END_REF] has the single equilibrium (Ξ = {0}).

The matrix A 1 has all eigenvalues with a positive real part and the matrix A 2 has eigenvalues with zero real part (two with non-zero imaginary parts). For the continuous dynamics F consider a Lyapunov function W (x) = |x| 2 , then

Therefore, for any initial condition x(t , j) ∈ O, t ∈ R + , j ∈ Z + for solutions of the system (11) always there exists a time instant t ≥ t such that x(t , j) ∈ D. Indeed, assume that it is not, then from the above expression of Ẇ , for W ≥ 1 the system ẋ = F(x) is purely oscillating and always there exists a crossing of D for |k| = 1, for W < 1 we have

and in a finite time the set with W ≥ 1 is reached. Therefore, t exists and Assumption 1 is valid. In addition, since for W ≥ 1 the system ( 9) between jumps has a linear oscillating dynamics with the frequency √ 2, then there is the maximal time between jumps τ max = √ 2π.

, the continuous dynamics F has the derivative V = 0 for all x ∈ C ∩ X A and

where α(s) = {s + 1 -ϕ √ 1 + s 2 [1 + s]}/τ max , then all conditions of Theorem 2 are satisfied and the set A is NZ preasymptotically stable for [START_REF] Yakubovich | Oscillations in systems with discontinuous and hysteresis nonlinearities[END_REF]. Choose U (x) = W (x), then the condition ( 9) is clearly valid for all x ∈ X U ∩ C since Ẇ = 2W , the time domain condition (8) can also be checked for this system, indeed define W j = W [φ(t j , j)], then (ϕ(|x|) = 1

= φ(t j , j) T e (A T 1 +A1)(tj+1-tj ) φ(t j , j);

Therefore,

where I is the identity matrix and the inequality e η -1 ≥ η for η ≥ 0 was used in the last step. Thus all conditions of Theorem 3 (the second part) are verified and the set U is NZ strictly prerepulsive for the system [START_REF] Yakubovich | Oscillations in systems with discontinuous and hysteresis nonlinearities[END_REF]. Finally, all conditions of Theorem 4 are satisfied and the system [START_REF] Yakubovich | Oscillations in systems with discontinuous and hysteresis nonlinearities[END_REF] has an oscillating solution.

The system trajectories for k = 1.5 are shown in Fig. 3, it has a limit cycle with four jumps.


Riccardo Spica
email: riccardo.spica@irisa.fr

Paolo Robuffo Giordano

Markus Ryll

Heinrich H Bülthoff

Antonio Franchi
email: antonio.franchi@tuebingen.mpg.de

An Open-Source Hardware/Software Architecture for Quadrotor UAVs

In this paper, we illustrate an open-source ready-to-use hardware/software architecture for a quadrotor UAV. The presented platform is price effective, highly customizable, and easily exploitable by other researchers involved in high-level UAV control tasks and for educational purposes as well. The use of object-oriented programming and full support of Robot Operating System (ROS) and Matlab Simulink allows for an efficient customization, code reuse, functionality expansion and rapid prototyping of new algorithms. We provide an extensive illustration of the various UAV components and a thorough description of the main basic algorithms and calibration procedures. Finally, we present some experimental case studies aimed at showing the effectiveness of the proposed architecture.

INTRODUCTION

Quadrotors are Unmanned Aerial Vehicles (UAVs) actuated by four fixed-pitch independent propellers located at the vertexes of a cross-shaped structure [START_REF] Mahony | Multirotor Aerial Vehicles: Modeling, Estimation, and Control of Quadrotor[END_REF]. In addition to being able to take-off and land vertically, quadrotors can reach high angular accelerations thanks to the relatively long lever arm between opposing motors. This makes them more agile than most standard helicopters or similar rotorcraft UAVs. For these reasons, and also because of their affordability and mechanical robustness, this platform has witnessed a considerable growing interest in the robotics and hobbyist community over the last decade, being, in some cases, also commercialized as off-the-shelf products [1, asc].

Nowadays, most of the research efforts are not dealing with more accurate dynamical modeling, or novel flight control design, as these issues have reached an adequate level of maturity and do not need further major improvements, at least concerning standard applications. The challenge is rather on how to exploit these systems as a flexible platform for implementing and validating complex and higherlevel tasks, often involving multiple robots interacting at the same time in partially structured or unstructured environments. Some examples in this sense can be found in, e.g., [START_REF] Schwager | Optimal coverage for multiple hovering robots with downward facing cameras[END_REF][START_REF] Michael | Cooperative manipulation and transportation with aerial robots[END_REF][START_REF] Franchi | Modeling and control of UAV bearing-formations with bilateral high-level steering[END_REF].

Nevertheless, obtaining a reliable and operative quadrotor platform for experimental purposes still requires a considerable investment of time and resources for properly tuning the system, e.g., for taking care of identification and calibration of all the various parameters (offsets, biases, motor curves), gain tuning of the flight controller, implementation of accurate onboard state estimation filters, and, last but not least, software development and debugging. These issues are obviously common to any robot, but, we believe, are particularly sensitive for flying robots as the malfunctioning of any hardware/software component can easily result in a severe crash and loss of the vehicle. Because of these reasons, many ready-to-use platforms have been proposed over the last years, each of them with their pros/cons and spanning a wide range of prices.

The Hummingbird and Pelican quadrotors by Ascending Technologies are probably the most widespread platforms in the research community (see, e.g., [START_REF] Schwager | Optimal coverage for multiple hovering robots with downward facing cameras[END_REF][START_REF] Michael | Cooperative manipulation and transportation with aerial robots[END_REF]). They are characterized by a solid flight performance thanks to their lightweight rigid structure, but their price is still expensive compared to other solutions. An alternative consumeroriented and cheaper quadrotor platform, employed in some research projects, is the AR.Drone Parrot [START_REF] Bristeau | The Navigation and Control technology inside the AR.Drone micro UAV[END_REF] which is able to achieve a stable near hovering flight only using its onboard sensors. However, being the AR.Drone an inexpensive quadrotor designed for the general consumer market, it possesses clearly some drawbacks when used within a research project, e.g., 1) it does not allow for the tuning and customization that are often required in research projects, 2) it is not possible to integrate additional sensors and high performance computational units (due to the extremely low payload and fragile structure), and 3) the usability of the integrated sensors for testing new algorithms is limited by their quality and the latencies at which the corresponding signals can be accessed. Furthermore, both these architectures are closed-source and the low-level controller cannot be freely customized to comply with the users' needs.

As for open-source architectures, an interesting survey of some existing projects can be found in [START_REF] Lim | Build Your Own Quadrotor: Open-Source Projects on Unmanned Aerial Vehicles[END_REF]. Many of the referenced solutions are extremely competitive when compared to their price. Nevertheless, in our opinion none of them possesses enough computational resources to run onboard the typical complex algorithms required for a truly autonomous operation of quadrotor UAVs, such as, e.g., intensive image processing and/or high load communication management. Another project, similar to the one presented in this paper, is described in [START_REF] Sa | System identification, estimation and control for a cost effective open-source quadcopter[END_REF]. The authors exploit the same hardware used in this work, but they do not customize the low-level controller provided by the manufacturer and strongly rely on a ground-based computer for all computations.

With respect to all these solutions, the goal of this paper is to propose a open-source and relatively cheap architecture for a quadrotor UAV characterized by the basic hardware and software components that can be used within an educational and research project to perform autonomous flight and onboard implementation of complex behaviors. For this purpose we are making the complete source code and specs available on a freely-accessible repository for dissemination purposes and easy re-use by other research groups (see https://svn.kyb.mpg.de/kyb-robotics/). From a hardware point of view, the proposed platform is purposely simple, highly customizable and easy to maintain. Every component can be bought off-the-shelf, and the use of standard software architectures and protocol allows for some flexibility in the hardware setup. As from a software point-of-view, the use of object-oriented programming and the support of the Robot Operating System (ROS) and Matlab Simulink environments allow for efficient customization, code reuse, functionality expansion and rapid prototyping of new algorithms. Finally, for the sake of completeness and re-usability, we also give handson details on all the calibration procedures necessary for identifying the dynamic parameters of the quadrotor, and on its basic flight control algorithms.

The paper is organized as follows: Secs. 2-3 describe the general hardware and software architecture of the proposed quadrotor platform. Then, Secs. 4-6 provide details on the implemented trajectory tracking controller, relevant calibration procedures, and state estimation schemes. Finally, Sec. 7 reports some experimental results for some relevant test cases meant to illustrate the potential and performance of the proposed platform, and Sec. 8 draws some conclusions and discusses future directions.

HARDWARE DESIGN

The main goals driving our hardware design can be summarized as follows:

• keeping the cost of the overall components as low as possible, when compared to equivalent solutions, • guaranteeing a wide availability of all the components by avoiding the use of customized hardware, • ensuring an adequate onboard computational power, for the reduced payload of a typical small-size quadrotor platform, • ensuring that both the actuation and sensing capabilities of the platform are made fully interfaced with the ROS environment.

Taking into account the aforementioned specifications, we found a suitable solution in using the mechanical frame, actuators, microcontrollers, and inertial measurement unit (IMU) of the MK-Quadro1 platform, by nevertheless replacing the control layer provided by the manufacturer, introducing new components, and completely redesigning the interface to the platform actuation and sensing resources.

The actuation system of the MK-Quadro consists of four plastic propellers with a diameter of 0.254 m, and a total span and weight of the frame of 0.5 m and 0. In addition to this, the platform is also equipped with a pressure sensor. The whole system is powered by a 2600 mAh LiPo battery which guarantees an endurance of around 10 min of flight in normal regimes. The complete system has a weight of approximately 1.300 kg.

On top of this basic setup the system permits the easy addition of onboard exteroceptive sensor, as done, e.g., in [START_REF] Grabe | On-board velocity estimation and closed-loop control of a quadrotor UAV based on optical flow[END_REF][START_REF] Franchi | Modeling and control of UAV bearing-formations with bilateral high-level steering[END_REF] by using onboard monocular cameras. The use of an onboard RGB-D sensor is currently under development.

We finally note that this hardware solution should be regarded as a suggested configuration, i.e., as a well tested and consolidated, rather than as a mandatory one. Indeed, other hardware choices are possible: in particular, the high-level control software can be run on any standard x86/64 architecture without additional effort. Moreover, since the source code is freely distributed, it is possible (with, probably, some more effort in this case) to compile and run it on different architectures. Finally, concerning the mechanical components, any equivalent configuration can be used by suitably adjusting the control/estimation gains.

SOFTWARE ARCHITECTURE

The software solution developed for our platform is meant to comply with the following specifications:

• possibility of distributing the full open-source code,

• use of an object-oriented high-level programming language (such as C++ or Python), • possibility of providing the user with a working system in terms of flight control, making the platform ready-to-use in laboratory conditions so as to, e.g., test new state estimation and planning algorithms for exploration or multi-robot formation controllers, • possibility to customize, if desired, the single basic functionalities, such as the controllers or the basic state estimator provided with the platfrom,

Base station

Exterocept. • possibility to automatize as much as possible the calibration procedure, • easy integration of additional functionalities in the architecture, and easy interface to Matlab and Simulink.

Sensors

In order to comply with the first two constraints, we opted for Ubuntu Linux Server 12.04 LTS as operating system, as this distribution is currently among the most widespread ones. Furthermore, for facilitating the integration and/or replacement of code parts, we adopted the ROS framework: indeed, this is more and more becoming a standard for robotic applications, with many common functionalities and drivers freely available on the web.

The overall system architecture is represented by the block diagram in Fig. 1. The whole architecture is developed within the the TeleKyb framework (see [START_REF] Grabe | The TeleKyb framework for a modular and extendible ROS-based quadrotor con-trol[END_REF] for more details). TeleKyb manages the generation of the desired trajectories from a pool of behaviors and selects the current state estimator and the motion controller. The control loop frequency depends on the kind of controller used and may vary, e.g., from 60 Hz to 650 Hz. Both the state estimation and the control algorithms are incapsulated in separate plugins. It is therefore easy, if desired, to re-implement them in a customized way.

Behaviors describe the basic operations available for the UAV (e.g., taking-off, landing, following a trajectory, executing velocity commands from a joystick, and so on) at a high level and in a parametrized way. Selection of the behaviors is governed by the behavior controller according to the robot state, to the experimental flow (managed by the experiment node), and to the joystick inputs from the user.

The experiment node is an independent ROS node running in a remote base station and communicating with the Atom board through a standard WiFi connection. From here, the user can control the experiment, modify the parameters online and log data. Finally the base station can also run Matlab Simulink models able to exchange data with the other software components via ROS publishers and subscribers instantiated in s-function blocks.

The MkInterface node is responsible for the communication with the low-level Atmega1284p microcontroller throughout the two serial ports. It receives the IMU measurements, the barometer measurements, the battery status, and the value of other internal variables. On the other hand it allows both to send control commands and change the logical states of the microcontroller. This latter can operate in two different modalities:

• In the near-hovering mode, the microcontroller receives, from the high-level control loop, a desired total thrust force, roll and pitch angles and yaw rate.

It regulates these quantities using the information provided by the IMU and employing a standard linearized hovering controller that is described in [START_REF] Lee | Semi-autonomous haptic teleoperation control architecture of multiple unmanned aerial vehicles[END_REF]; • In the direct mode, instead, the microcontroller is only acting as a communication interface with the sensory and actuation hardware. It collects the IMU readings and then forwards the motor velocity commands (from the high-level control loop) to the brushless (BL) controllers via the I 2 C bus.

Additionally, an emergency mode is also possible: transitions to this mode can be toggled by either the user or by the occurrence of dangerous situations such as, e.g., a very low battery level or an interruption of the communication link with the high-level control. In this case, the microcontroller ignores any further command from the Atom board and tries to land by only relying on the nearhovering controller.

Finally, in all cases the microcontroller runs an onboard state estimation based on a complementary filter as described in [START_REF] Mahony | Complementary filter design on the special orthogonal group SO(3)[END_REF]. This estimation is not strictly necessary in the direct mode but, in order to obtain a smooth transition to the emergency mode at any time, the estimation is nevertheless always kept updated.

MODELING AND CONTROL

In this section we provide details on the implemented control algorithm for trajectory tracking when the microcontroller is operating in direct mode. As typically done, we neglect the dynamics of the propeller actuation and consider the generated thrusts as the actual control inputs for control design. As for the control scheme itself, we chose to employ a geometric tracking technique similar to the one presented in [START_REF] Lee | Geometric tracking control of a quadrotor UAV on SE(3)[END_REF] because of its almost global convergence and excellent tracking capabilities.

Throughout the following, a right subscript will indicate the reference frame to which the quantities belong, and a left superscript the frame in which the quantities are expressed. An omitted left superscript will indicate a quantity expressed in the world inertial frame W. The mapping between the thrusts f i generated by each propeller along its spinning axis and the resulting total thrust u 1 along the vertical axis of the robot, and the three independent torques u 2,3,4 actuating the quadrotor body, is:

u =    u 1 u 2 u 3 u 4    =    1 1 1 1 0 l 0 -l -l 0 l 0 c -c c -c       f 1 f 2 f 3 f 4    = Mũ, (1)
where l represents the distance of the propeller rotational axes from the center of the robot frame and c is a constant related to the propeller characteristics.

The robot state consists of the position and orientation of the body frame B w.r.t. the world inertial frame W, namely r B and W R B , together with the corresponding time-derivatives, i.e. the linear and angular velocity ( ṙB and B ω BW ). Assuming that the robot barycenter is positioned in the geometric center, we obtain

mr B + mgz W = u 1 z B , (2)
J B ωBW + [ B ω BW ] × J B ω BW = (u 2 , u 3 , u 4 ) T , (3)
where m and J are the total mass and the (constant) inertia tensor of the robot respectively, z B is the Z-axis of the body frame and [ B ω BW ] × is the skew-symmetric matrix associated to vector B ω BW .

Given a desired trajectory (r B,t , ψ t ) for the robot position and yaw angle with all necessary derivatives, and having defined the position and velocity errors as e p = r B,tr B and e v = ṙB,t -ṙB , we can introduce a desired force as a PID action plus gravity cancellation

f d = mr B,t +K hl p e p +K hl p t 0 e p dt+K hl v e v +mgz W , (4)
where K hl p , K hl p and K hl v are positive diagonal matrices. The desired force f d can be realized by choosing u 1 = f T d z B and by aligning the body vertical axis with the direction of f d . The desired robot attitude, described by the rotation matrix W R B,d , can then be chosen so that its third column is z B,d = f d f d . The remaining degrees of freedom are then set so as to minimize the yaw error. This is obtained by defining

y C,d = (-sin ψ t cos ψ t 0) T ,
and by choosing the first and second columns of W R B,d as:

x B,d = y C,d × z B,d y C,d × z B,d , y B,d = z B,d × x B,d .
The orientation error can then be defined as

e R = 1 2 W R T B W R B,d -W R T B,d W R B ∨ ,
where ∨ indicates the vee-map that relates a skewsymmetric matrix to its corresponding 3-vector. The angular velocity error is instead

e ω = W R T B W R B,d B,d ω BW,d -B ω BW
where B,d ω BW,d can be computed as a function of the specified trajectory of the robot barycenter and its derivatives up to the third order (see [START_REF] Spica | Planning and control for aerial grasping with a quadrotor UAV[END_REF]). Finally the torque inputs can be chosen as

(u 2 , u 3 , u 4 ) T = K hl R e R + K hl R t 0 e R dt + K hl ω e ω + [ B ω BW ] × J B ω BW + J W R T B W R B,d B,d ωBW,d -[ B ω BW ] × W R T B W R B,d B,d ω BW,d (5)
where K hl R , K hl R and K hl ω are positive diagonal matrices. Note that this control law essentially consists of a PID action, a cancellation of the gyroscopic force and a feedforward term B,d ωBW,d that can be computed as a function of the specified trajectory of the robot barycenter and its derivatives up to the fourth order (see again [START_REF] Spica | Planning and control for aerial grasping with a quadrotor UAV[END_REF]).

At this point, since matrix M in (1) has full rank, it can be inverted to obtain the single motor thrusts f i . As well-known, and shown in Sec. 5, these are related to the spinning velocities of the propellers which are then treated as references for the motor controllers. Nevertheless, directly feeding these references may conflict with any internal saturations due to the limited motor power. The resulting actuation can therefore significantly differ from the desired one. It is then necessary to properly scale down the control inputs in order to make them feasible by the motor controllers. We can achieve this goal by splitting the reference control input u in a component u g = (mg, 0, 0, 0) T , necessary for compensating gravity in hovering, and a second component u res = uu g on which a uniform scaling is applied until all the resulting propeller thrusts lie in the admissible range [f min , f max ]. In practice one has to solve the following optimization problem:

min {λ} 1 2 u -(u g + λu res ) 2 s.t. f min ≤ A -1 (u g + λu res ) ≤ f max , λ ≤ 1.
Note that an admissible solution always exists, provided that u g is a feasible input. Moreover the use of a uniform scaling guarantees that the direction of the control input is maintained.

CALIBRATION AND IDENTIFICATION

A necessary condition for obtaining good flight performance is a proper calibration of all the dynamic parameters of the system. In this section we describe how one can easily obtain an accurate estimation of: i) the accelerometer bias, ii) the gyroscope bias, iii) the motor/propeller characteristics iv) the mass and inertia matrix of the robot. For the sake of simplicity we will assume that the world Z-axis is aligned with the direction of gravity.

Gyroscope bias calibration

The bias of the gyroscope is defined as the output of the sensor when not experiencing any rotation. This bias is influenced by many variable factors like, e.g. the temperature, and it is then affected by considerable drift. As a consequence the calibration of this sensor must be repeated often enough. We therefore decided to completely automate this operation: each time the MkInterface node is restarted, a fixed number N of gyro samples are registered. The mean of the readings is then taken as a bias, provided that their standard deviation is below a certain threshold. Otherwise an error is thrown.

Accelerometer bias calibration

As well-known, an accelerometer measures its acceleration w.r.t. the inertial frame apart from any gravity component. Any accelerometer is affected by a bias, defined as the output of the sensor during free-fall under gravity. In addition, a scaling factor and a roto-translation w.r.t. the UAV body frame can also be present after being mounted. Differently from the gyroscope case, the accelerometer bias is almost constant over time, and it is thus not necessary to repeat the calibration unless the mounting of the sensor is changed. To describe the employed calibration procedure, let rB = W R B T (a mb) + g where b is a bias vector, matrix T takes into account any rotation and scaling factor, and g is the nominal value of the gravity acceleration. When the robot is at rest (r

B = 0), it is Ta m -Tb = -W R T B g.
Assuming that a measure of W R B is available, this equation is linear in the unknowns T and Tb. One can then repeat N measurements of a m for different (constant) robot orientations and stack them in matrix form:

   a T m,1 -1 . . . . . . a T m,N -1    T T b T T T =    -g T W R B,1 . . . -g T W R B,N    .
If the left matrix has rank 4, a simple pseudoinversion solves the system in a least-square sense. We note that this requires at least four available measurements for different orientations of the accelerometer.

As a final remark, we note that a similar rotation and scaling factor are also in general present for the gyroscope. Nevertheless the estimation of these quantities is more difficult because of the absence of a reference ground truth, as it is gravity for the accelerometer case. On the other hand, differently from the case of the accelerometer, any remaining miscalibration in the gyroscope sensor was experimentally found to be negligible.

Motor calibration

We identified the static characteristics of the brushless controller/motor/propeller system by using a Nano17 force/torque sensor 2 . The identification of the static characteristic was obtained by feeding the brush-less controller with the reference signal represented (in green) in Fig. 2a. We recall that this reference is linearly proportional to the desired rotational speed of the propeller. For each step, the motor reference was kept constant for 5 s, allowing for an attenuation of the noise by averaging the corresponding sensor measurements. Between each phase, a resting phase of 10 s was introduced with a low rotational speed for the propeller. This was necessary in order to avoid a propeller overheating which would have distorted the results. The reference was designed in order to vary from 30 to 190 in steps of 10.

As expected, the relationship between the rotational speed and the generated force/torque can be well approximated by a quadratic function. Since we are interested in the inverse relation, necessary to generate the reference for the brush-less controllers, we report the fitting curve in Fig. 2b with the force as the independent variable. We 2 http://www.ati-ia.com decided to exclude from the calibration any reference signal smaller than 30 or greater that 180 (resulting in a force of 0.3131 and 6.03 N respectively). This also allows to obtain a better fitting in the range most likely encountered during normal flight. We also note that, with the chosen limitations on the maximum rotor speeds, and considering the weight of the robot, the motors approximately operate in the center of their dynamical range during hovering.

The ratio between generated forces and torques is also quite constant in the range of interest, as well-known in the literature, and results about 0.0174 m.

Finally, the dynamics between commanded and actual propeller speeds was found to be well approximated by a first order linear system with a time constant of 0.047 s.

Even if present, we decided to neglect this effect in the control design.

Estimation of the mass and inertia tensor

In order to estimate the robot mass, we made use of an algorithm taken from classical adaptive control techniques [START_REF] Slotine | Applied nonlinear control[END_REF]. We start by considering the scalar vertical dynamics of a UAV, obtained by projecting (2) along z W . This can be rearranged as:

z -g = 1 m u 1 z T B z W = 1 m υ = m inv υ. (6)
For convenience, we will estimate the inverse of the mass m inv = 1/m. Equation ( 6) depends on the vertical acceleration z which, in practice, is usually a very noisy quantity.

To obtain a dependence in terms of the linear vertical velocity ż, we multiply both sides of ( 6) by a low-pass filter P (s) = a/(s + a) P (s)(z -g) = m inv P (s)υ.

(7) The term P (s)z in [START_REF] Lim | Build Your Own Quadrotor: Open-Source Projects on Unmanned Aerial Vehicles[END_REF] can be expanded as

a s + a z = as s + a ż = a ż - a s + a ż = a ( ż -P (s) ż) , (8) i
.e., as a function of the sole velocity ż. Now, let g f il = P (s)g, υ f il = P (s)υ, and żfil = P (s) ż. Equation ( 7) can be rewritten as

a( ż -żfil ) -g f il = m inv υ f il . (9)
Assume an estimation of the (inverse of the) mass minv is available, and let the lhs of (9) be y = a( ż -żfil ) -g f il and the rhs of (9) evaluated on the estimated mass be ŷ = minv υ f il . Defining the error function

V (e) = 1 2 k est (y -ŷ) 2 , k est,m > 0,
one can minimize it in terms of m inv by implementing the following gradient update rule

ṁinv = - ∂V (e) ∂ minv = k est,m eυ f il . (10)
Note that the implementation of (10) requires knowledge of i) the vertical (world frame) velocity ż, ii) the direction of the axis z B , iii) the commanded thrust u 1 , and iv) the evaluation of the low-pass filtered signals υ f il , g f il and żfil . In our implementation, the following values were used: a = 0.06π and k est,m = 1.

Finally, we note that the same algorithm can be used for obtaining an estimation of the robot inertia tensor. Indeed, because of the cylindrical symmetry of the robot, the inertia matrix results practically diagonal and with the same X and Y components. Under this assumption the gyroscopic term in (3) becomes negligible for a small Z component of B ω BW . The rotational dynamics can then be split among three independent scalar systems similar to [START_REF] Lee | Geometric tracking control of a quadrotor UAV on SE(3)[END_REF], with the additional simplification of the absence of any gravity term.

STATE ESTIMATION

The control algorithm described in Sec. 4 requires knowledge, at high rate, of the full state of the robot, i.e., of the position r B , the velocity ṙB , the orientation W R B and the angular velocity B ω BW . In this section we describe how to obtain an estimation of these quantities using measurements of linear acceleration and angular velocity available from the onboard IMU, and a low rate measurements of the position and orientation provided by any possible source (e.g., an external or onboard sensor).

Complementary filter in SE(3)

If some (possibly noisy) measurement of the robot orientation is available, the translational part of the robot dynamics χ T = (r B , ṙB ) T has the linear expression

χT = d dt r B ṙB = 0 3×3 I 3×3 0 3×3 0 3×3 χ T + 0 3×1 W R B T (a m -b) + g = Aχ T + υ.
By assuming a low-rate and possibly noisy measurement of the position r B is available, it is straightforward to check that the system is observable and a simple linear state observer can be easily designed as:

χT = A χT + υ + K est,p (r B -rB ) . (11)
We note that one can also estimate the accelerometer bias b. Indeed by considering the extended state χ ′ T = (r B , ṙB , b)

T , one has

χ′ T = d dt r B ṙB b = 0 3×3 I 3×3 0 3×3 0 3×3 0 3×3 -W R B T 0 3×3 0 3×3 0 3×3 χ ′ T + + 0 3×1 W R B Ta m + g 0 3×1 = A ′ χ T + υ ′ .
.

It is easy to prove observability w.r.t. the position r B measurement in this case, too, thus allowing to implement an observer similar to [START_REF] Mahony | Multirotor Aerial Vehicles: Modeling, Estimation, and Control of Quadrotor[END_REF].

As for the orientation, we implemented the complementary filter in SO(3) described in [START_REF] Mahony | A study of non-linear complementary filter design for kinematic systems on the special orthogonal group[END_REF], in particular the socalled passive complementary filter without gyro bias estimation (already compensated for during the calibration of Sec. 5), and using the quaternion-based formulation. Let W q B = ( W q B,s , W q B,v ) be the unit-norm quaternion corresponding to W R B with scalar and vector part W q B,s and W q B,v , respectively. As well-known, the kinematics of the quaternion is given by:

W qB = 1 2 W q B ⊗ 0, B ω BW
where ⊗ indicates the quaternion product defined as:

p ⊗ q = p s q s -p T v q v , p s q v + q s p v + p v × q v .
Assuming that a low-rate and possibly noisy measure of the orientation W q B is available, we can define an orientation error as δq = W q * B ⊗ W q B = (δq s , δq) where the conjugate is defined as q * = (q s , -q v ). The estimation of the orientation can then be updated as follows:

W qB = 1 2 W qB ⊗ 0, B ω BW + 2K est,q δq s δq v . (12)
In order to enforce the unit-norm constraint during the numerical integration, we slightly modified [START_REF] Michael | Cooperative manipulation and transportation with aerial robots[END_REF] by adding a penalty term proportional to the squared norm of the estimated quaternion:

W qB = 1 2 W qB ⊗ 2λ(1 -W qB 2 ), ω m + 2K est,q δq s δq v
where ω m is the angular velocity measured by the gyroscope corrected for the calibrated bias. Finally, we employed the following values: K est,p = 126.3I 3 , K est,v = 3987.3I 3 , and K est,q = 62.8I 3 .

The last part of the state to be estimated is the angular velocity B ω BW . In principle this quantity is already provided by the onboard gyroscope. In practice the output of this sensor during a typical flight is too noisy to be directly fed to the controller. Therefore, we low-pass filtered B ω BW as ω = 2πf ω (ω m -ω) with f ω = 150 Hz.

Extension to bearing-only measurements

The state estimation algorithm described above requires a low-rate and possibly noisy direct measurement of the position and orientation of the robot. Since these measurements are not always available in a laboratory, we here detail a ready-to-use algorithm able to recover these quantities from a set of (at least 4) bearing measurements from some beacons with known identity and position in the world frame. This is, indeed, representative of many situations involving an onboard camera extracting features from the scene.

The reconstruction of the robot pose can be obtained by exploiting standard techniques from geometrical computer vision. The basic algorithm, described in [START_REF] Ma | An invitation to 3D vision, from images to models[END_REF], is here briefly summarized: first, consider a set of N points p belonging to the same 2D plane and assume that two bearing measurements are available for each point, namely x i,1 and x i,2 for i = 1 . . . N . These two bearings are associated to two frames (in our case the world and the body frame respectively). The plane is described by the following equation in the first reference frame:

P = X 1 : n T 1 X 1 = d ⇔ 1 d n T 1 X 1 = 1 (13)
where X 1 contains the coordinates of the point p in the first frame, n 1 is the unit normal vector to the plane and d is the distance of the plane from the origin of the first frame. Note that, since we have assumed a known beacon position in the world frame, all these quantities are also available.

Letting ( W R B , r B ) represent the rotation/translation among the two frames, and exploiting (13), one has

X 2 = W R B + 1 d r B n T 1 X 1 = HX 1
where H is the so-called homography matrix. The homography H can be directly retrieved from a set of at least four bearing pairs (in the reference and current frame) by means of the standard four-point algorithm for planar scene as described in [START_REF] Ma | An invitation to 3D vision, from images to models[END_REF]. Once recovered, H can then be decomposed back into the rotation W R B and the translation r B : let (n ′ 1 , n ′′ 1 ) be a pair of unit-norm vectors such that (n 1 , n ′ 1 , n ′′ 1 ) forms a right hand frame. Note that (n ′ 1 , n ′′ 1 ) are a basis of the null space of n 1 T . By projecting the homography matrix along these vectors, we obtain

H (n ′ 1 n ′′ 1 ) = (m ′ 1 m ′′ 1 ) = W R B (n ′ 1 n ′′ 1 ) . By then taking m 1 = m ′ 1 × m ′′ 1 , W R B is finally recovered as W R B = (m 1 m ′ 1 m ′′ 1 ) (n 1 n ′ 1 n ′′ 1 ) T .
Once W R B has been found, t is simply given by t = d H -W R B where d is assumed known, as explained above. This then allows to recover a measurement of the robot pose ( W R B , r B ) that can be used in the complementary filter described in Sec. 6.1 .

EXPERIMENTAL CASES STUDIES

In this section we present some experimental results obtained with the proposed platform. In the first two experiments we used an external tracking system 3 to provide pose measurements directly to the estimator described in 6.1. For the last experiment, instead, we used the external tracking system only for simulating low-rate and noisy bearing measurements similar to those that would be provided by typical onboard localization systems (e.g., vision-based) available in the literature.

The aim of the first experiment was to prove the performance of the mass estimation algorithm in hovering conditions and during arbitrary motion. To this end, we intentionally initialized the system from an initial guess for the estimated mass of 1 kg, i.e. 0.308 kg less than the nominal value. The result of the experiment is shown in Fig. 3. At the beginning of the experiment, the controller commanded a vertical thrust not sufficient for the robot to take off due to the low value of the estimated mass.

After the take-off, the robot was then commanded to hover at a constant position until the estimated mass had reached its final value (about 1.394 kg). The convergence of the estimation can be detected by setting a threshold on its variance. We then manually controlled the cartesian velocity of the robot in a random way for verifying that the mass estimation was remaining constant despite the erratic motion. Figure 3c shows the actual robot velocity as estimated via the techniques illustrated in Sec. 6.

In the second experiment we show the tracking performance of our control algorithm while following an 8-shaped horizontal trajectory in X and Y (see Fig. 4). The 8 shape has a bounding box of size 1 m×1 m and is completed in 8 s. Both the Z position and the yaw reference were varied during motion by following sinusoidal-like trajectories with an amplitude of 0.2 m and frequency of 1 rad, respectively. The timing law was chosen to ensure the needed smoothness for the reference sent to the controller, i.e., continuous up to the fourth order for the position and up to the second order for the yaw angle. In Fig. 5 some snapshots of this experiment are shown and correspond to the phases of maximum tilt angles for the UAV, i.e. at the four sharp curves of the 8-shape.

As a final experiment, we report in Fig. 6 the results of the state estimation relying on bearing measurements. For this experiment the desired robot trajectory was generated manually with a joystick. We emulated the presence of 4 virtual beacons placed at 10 cm above ground and forming a square of 2 m size. The bearing measurement rate was set to 30 Hz with an additional gaussian angular noise with 1.5 deg STD. The pose of the robot was calculated from the simulated bearings by using the algorithm described in Sec. 6.2 and implemented in a Simulink model. From the comparison with the ground truth in the plots of Fig. 6, one can note that the pose estimation algorithm based on the Homography decomposition was still able to reconstruct a high-quality estimation of the robot state, and consequently to allow for a successful flight.

CONCLUSIONS

In this paper we have presented a completely opensource, ready-to-use and cheap platform for prototyping complex algorithms for single or multiple quadrotor UAVs, and relying on only onboard hardware. The platform can be fully customized from high-level behaviors to low-level motor control. The hardware and software architecture was described in detail with all the needed calibration and tuning procedures. Furthermore, we also freely distribute the whole software package for the interested reader at the SVN repository https://svn.kyb.mpg.de/kyb-robotics/. This allows an easy porting of our solution to any other lab interested in developing quadrotor applications for research purposes.

Fig. 1 .

 1 Fig. 1. Block diagram of the hardware and software architecture. The TeleKyb framework is used for the software part.

Fig. 2 .

 2 Fig. 2. Motor calibration: input/output signals (a) and fitting model (b).

3Fig. 3 .

 3 Fig. 3. Mass estimation during a normal flight. Fig. (a): estimated mass. Fig. (b): UAV position. Fig. (c) and UAV velocity.

Fig. 4 .

 4 Fig. 4. Tracking of an eight shaped trajectory. Fig. (a): UAV cartesian position. Fig. (b): UAV yaw angle. Fig. (c): position error. Fig. (d): yaw angle error.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Snapshots of the quadrotor executing an 8-shaped trajectory

http://www.mikrokopter.de

ACKNOWLEDGEMENTS

The authors want to thank Martin Riedel for his support with the TeleKyb framework.